49 #ifndef __INTREPID2_PROJECTIONTOOLSDEFHDIV_HPP__
50 #define __INTREPID2_PROJECTIONTOOLSDEFHDIV_HPP__
58 namespace Experimental {
61 template<
typename ViewType1,
typename ViewType2,
typename ViewType3>
63 const ViewType1 sideBasisNormalAtBasisEPoints_;
64 const ViewType1 basisAtBasisEPoints_;
65 const ViewType2 basisEWeights_;
66 const ViewType1 wBasisDofAtBasisEPoints_;
67 const ViewType2 targetEWeights_;
68 const ViewType1 basisAtTargetEPoints_;
69 const ViewType1 wBasisDofAtTargetEPoints_;
70 const ViewType3 tagToOrdinal_;
71 const ViewType1 targetAtEPoints_;
72 const ViewType1 targetAtTargetEPoints_;
73 const ViewType1 refSidesNormal_;
74 ordinal_type sideCardinality_;
75 ordinal_type offsetBasis_;
76 ordinal_type offsetTarget_;
77 ordinal_type sideDim_;
82 const ViewType1 basisAtBasisEPoints,
const ViewType2 basisEWeights,
const ViewType1 wBasisDofAtBasisEPoints,
const ViewType2 targetEWeights,
83 const ViewType1 basisAtTargetEPoints,
const ViewType1 wBasisDofAtTargetEvalPoint,
const ViewType3 tagToOrdinal,
84 const ViewType1 targetAtEPoints,
const ViewType1 targetAtTargetEPoints,
85 const ViewType1 refSidesNormal, ordinal_type sideCardinality, ordinal_type offsetBasis,
86 ordinal_type offsetTarget, ordinal_type sideDim,
87 ordinal_type dim, ordinal_type iside) :
88 sideBasisNormalAtBasisEPoints_(sideBasisNormalAtBasisEPoints),
89 basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisDofAtBasisEPoints_(wBasisDofAtBasisEPoints), targetEWeights_(targetEWeights),
90 basisAtTargetEPoints_(basisAtTargetEPoints), wBasisDofAtTargetEPoints_(wBasisDofAtTargetEvalPoint),
91 tagToOrdinal_(tagToOrdinal), targetAtEPoints_(targetAtEPoints),
92 targetAtTargetEPoints_(targetAtTargetEPoints),
93 refSidesNormal_(refSidesNormal), sideCardinality_(sideCardinality), offsetBasis_(offsetBasis),
94 offsetTarget_(offsetTarget), sideDim_(sideDim), dim_(dim), iside_(iside)
98 KOKKOS_INLINE_FUNCTION
99 operator()(
const ordinal_type ic)
const {
102 for(ordinal_type j=0; j <sideCardinality_; ++j) {
103 ordinal_type jdof = tagToOrdinal_(sideDim_, iside_, j);
105 for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq) {
106 for(ordinal_type d=0; d <dim_; ++d)
107 sideBasisNormalAtBasisEPoints_(ic,j,iq) += refSidesNormal_(iside_,d)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
108 wBasisDofAtBasisEPoints_(ic,j,iq) = sideBasisNormalAtBasisEPoints_(ic,j,iq) * basisEWeights_(iq);
110 for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq) {
111 typename ViewType2::value_type sum=0;
112 for(ordinal_type d=0; d <dim_; ++d)
113 sum += refSidesNormal_(iside_,d)*basisAtTargetEPoints_(ic,jdof,offsetTarget_+iq,d);
114 wBasisDofAtTargetEPoints_(ic,j,iq) = sum * targetEWeights_(iq);
118 for(ordinal_type d=0; d <dim_; ++d)
119 for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq)
120 targetAtTargetEPoints_(ic,iq) += refSidesNormal_(iside_,d)*targetAtEPoints_(ic,offsetTarget_+iq,d);
125 template<
typename ViewType1,
typename ViewType2,
typename ViewType3,
126 typename ViewType4,
typename ViewType5>
128 const ViewType1 basisCoeffs_;
129 const ViewType2 negPartialProjAtBasisEPoints_;
130 const ViewType2 nonWeightedBasisAtBasisEPoints_;
131 const ViewType2 basisAtBasisEPoints_;
132 const ViewType3 basisEWeights_;
133 const ViewType2 wBasisAtBasisEPoints_;
134 const ViewType3 targetEWeights_;
135 const ViewType2 basisAtTargetEPoints_;
136 const ViewType2 wBasisAtTargetEPoints_;
137 const ViewType4 computedDofs_;
138 const ViewType5 cellDof_;
139 ordinal_type numCellDofs_;
140 ordinal_type offsetBasis_;
141 ordinal_type offsetTarget_;
142 ordinal_type numSideDofs_;
144 ComputeBasisCoeffsOnCells_HDiv(
const ViewType1 basisCoeffs, ViewType2 negPartialProjAtBasisEPoints,
const ViewType2 nonWeightedBasisAtBasisEPoints,
145 const ViewType2 basisAtBasisEPoints,
const ViewType3 basisEWeights,
const ViewType2 wBasisAtBasisEPoints,
const ViewType3 targetEWeights,
146 const ViewType2 basisAtTargetEPoints,
const ViewType2 wBasisAtTargetEPoints,
const ViewType4 computedDofs,
const ViewType5 cellDof,
147 ordinal_type numCellDofs, ordinal_type offsetBasis, ordinal_type offsetTarget, ordinal_type numSideDofs) :
148 basisCoeffs_(basisCoeffs), negPartialProjAtBasisEPoints_(negPartialProjAtBasisEPoints), nonWeightedBasisAtBasisEPoints_(nonWeightedBasisAtBasisEPoints),
149 basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisAtBasisEPoints_(wBasisAtBasisEPoints), targetEWeights_(targetEWeights),
150 basisAtTargetEPoints_(basisAtTargetEPoints), wBasisAtTargetEPoints_(wBasisAtTargetEPoints),
151 computedDofs_(computedDofs), cellDof_(cellDof),numCellDofs_(numCellDofs), offsetBasis_(offsetBasis),
152 offsetTarget_(offsetTarget), numSideDofs_(numSideDofs) {}
155 KOKKOS_INLINE_FUNCTION
156 operator()(
const ordinal_type ic)
const {
158 for(ordinal_type j=0; j <numCellDofs_; ++j) {
159 ordinal_type idof = cellDof_(j);
160 for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq) {
161 nonWeightedBasisAtBasisEPoints_(ic,j,iq) = basisAtBasisEPoints_(ic,idof,offsetBasis_+iq);
162 wBasisAtBasisEPoints_(ic,j,iq) = nonWeightedBasisAtBasisEPoints_(ic,j,iq) * basisEWeights_(iq);
164 for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq) {
165 wBasisAtTargetEPoints_(ic,j,iq) = basisAtTargetEPoints_(ic,idof,offsetTarget_+iq)* targetEWeights_(iq);
168 for(ordinal_type j=0; j < numSideDofs_; ++j) {
169 ordinal_type jdof = computedDofs_(j);
170 for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq)
171 negPartialProjAtBasisEPoints_(ic,iq) -= basisCoeffs_(ic,jdof)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq);
177 template<
typename ViewType1,
typename ViewType2,
typename ViewType3,
178 typename ViewType4,
typename ViewType5,
typename ViewType6>
180 const ViewType1 basisCoeffs_;
181 const ViewType2 negPartialProjAtBasisEPoints_;
182 const ViewType2 nonWeightedBasisAtBasisEPoints_;
183 const ViewType2 basisAtBasisEPoints_;
184 const ViewType2 hcurlBasisCurlAtBasisEPoints_;
185 const ViewType3 basisEWeights_;
186 const ViewType2 wHCurlBasisAtBasisEPoints_;
187 const ViewType3 targetEWeights_;
188 const ViewType2 hcurlBasisCurlAtTargetEPoints_;
189 const ViewType2 wHCurlBasisAtTargetEPoints_;
190 const ViewType4 tagToOrdinal_;
191 const ViewType5 computedDofs_;
192 const ViewType6 hCurlDof_;
193 ordinal_type numCellDofs_;
194 ordinal_type offsetBasis_;
195 ordinal_type numSideDofs_;
199 const ViewType2 basisAtBasisEPoints,
const ViewType2 hcurlBasisCurlAtBasisEPoints,
const ViewType3 basisEWeights,
const ViewType2 wHCurlBasisAtBasisEPoints,
const ViewType3 targetEWeights,
200 const ViewType2 hcurlBasisCurlAtTargetEPoints,
const ViewType2 wHCurlBasisAtTargetEPoints,
const ViewType4 tagToOrdinal,
const ViewType5 computedDofs,
const ViewType6 hCurlDof,
201 ordinal_type numCellDofs, ordinal_type offsetBasis, ordinal_type numSideDofs, ordinal_type dim) :
202 basisCoeffs_(basisCoeffs), negPartialProjAtBasisEPoints_(negPartialProjAtBasisEPoints), nonWeightedBasisAtBasisEPoints_(nonWeightedBasisAtBasisEPoints),
203 basisAtBasisEPoints_(basisAtBasisEPoints), hcurlBasisCurlAtBasisEPoints_(hcurlBasisCurlAtBasisEPoints), basisEWeights_(basisEWeights), wHCurlBasisAtBasisEPoints_(wHCurlBasisAtBasisEPoints), targetEWeights_(targetEWeights),
204 hcurlBasisCurlAtTargetEPoints_(hcurlBasisCurlAtTargetEPoints), wHCurlBasisAtTargetEPoints_(wHCurlBasisAtTargetEPoints),
205 tagToOrdinal_(tagToOrdinal), computedDofs_(computedDofs), hCurlDof_(hCurlDof),numCellDofs_(numCellDofs), offsetBasis_(offsetBasis),
206 numSideDofs_(numSideDofs), dim_(dim) {}
209 KOKKOS_INLINE_FUNCTION
210 operator()(
const ordinal_type ic)
const {
212 ordinal_type numBasisEPoints = basisEWeights_.extent(0);
214 for(ordinal_type i=0; i<numSideDofs_; ++i) {
215 ordinal_type idof = computedDofs_(i);
216 for(ordinal_type iq=0; iq<numBasisEPoints; ++iq){
217 for(ordinal_type d=0; d <dim_; ++d)
218 negPartialProjAtBasisEPoints_(ic,iq,d) -= basisCoeffs_(ic,idof)*basisAtBasisEPoints_(ic,idof,offsetBasis_+iq,d);
222 for(ordinal_type i=0; i<numCellDofs_; ++i) {
223 ordinal_type idof = tagToOrdinal_(dim_, 0, i);
224 for(ordinal_type iq=0; iq<numBasisEPoints; ++iq) {
225 for(ordinal_type d=0; d<dim_; ++d)
226 nonWeightedBasisAtBasisEPoints_(ic,i,iq,d) = basisAtBasisEPoints_(ic,idof,offsetBasis_+iq,d);
230 for(ordinal_type i=0; i<static_cast<ordinal_type>(hCurlDof_.extent(0)); ++i) {
231 ordinal_type idof = hCurlDof_(i);
232 for(ordinal_type d=0; d<dim_; ++d) {
233 for(ordinal_type iq=0; iq<numBasisEPoints; ++iq) {
234 wHCurlBasisAtBasisEPoints_(ic,i,iq,d) = hcurlBasisCurlAtBasisEPoints_(idof,iq,d)*basisEWeights_(iq);
236 for(ordinal_type iq=0; iq<static_cast<ordinal_type>(targetEWeights_.extent(0)); ++iq) {
237 wHCurlBasisAtTargetEPoints_(ic,i,iq,d) = hcurlBasisCurlAtTargetEPoints_(idof,iq,d)*targetEWeights_(iq);
245 template<
typename SpT>
246 template<
typename BasisType,
247 typename ortValueType,
class ...ortProperties>
250 typename BasisType::ScalarViewType targetDivEPoints,
251 const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
252 const BasisType* cellBasis,
254 const EvalPointsType evalPointType){
255 typedef typename BasisType::scalarType scalarType;
256 typedef Kokkos::DynRankView<scalarType,SpT> ScalarViewType;
257 typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
258 auto refTopologyKey = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->
getTopologyKey());
260 ordinal_type dim = cellBasis->getBaseCellTopology().getDimension();
261 ordinal_type sideDim = dim-1;
263 ordinal_type numSides = cellBasis->getBaseCellTopology().getSideCount();
265 ordinal_type numCells = orts.extent(0);
268 typename CellTools<SpT>::subcellParamViewType subcellParamSide;
272 auto evalPointsRange = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->
getPointsRange(evalPointType));
274 for(ordinal_type is=0; is<numSides; ++is) {
275 ScalarViewType sideWorkview(
"sideWorkview", numCells, projStruct->
getMaxNumEvalPoints(evalPointType), sideDim);
277 const auto topoKey = refTopologyKey(sideDim,is);
278 auto sideBasisEPoints = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->
getEvalPoints(sideDim,is,evalPointType));
281 (
"Evaluate Points Sides ",
282 Kokkos::RangePolicy<SpT, int> (0, numCells),
283 KOKKOS_LAMBDA (
const size_t ic) {
284 auto sidePointsRange = evalPointsRange(sideDim, is);
285 auto sideRefPointsRange = range_type(0, range_size(sidePointsRange));
286 auto orientedBasisEPoints = Kokkos::subview(sideWorkview, ic, sideRefPointsRange, range_type(0,sideDim));
288 ordinal_type sOrt[6];
290 orts(ic).getFaceOrientation(sOrt, numSides);
292 orts(ic).getEdgeOrientation(sOrt, numSides);
293 ordinal_type ort = sOrt[is];
300 if(cellBasis->getDofCount(dim,0) <= 0)
303 auto evalDivPointsRange = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->
getDerivPointsRange(evalPointType));
304 auto cellDivPointsRange = evalDivPointsRange(dim, 0);
305 auto cellBasisDivEPoints = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->
getDerivEvalPoints(dim,0,evalPointType));
307 RealSpaceTools<SpT>::clone(Kokkos::subview(targetDivEPoints, Kokkos::ALL(), cellDivPointsRange, Kokkos::ALL()), cellBasisDivEPoints);
312 auto cellPointsRange = evalPointsRange(dim, 0);
313 auto cellBasisEPoints = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->
getEvalPoints(dim,0,evalPointType));
319 template<
typename SpT>
320 template<
typename basisCoeffsValueType,
class ...basisCoeffsProperties,
321 typename funValsValueType,
class ...funValsProperties,
323 typename ortValueType,
class ...ortProperties>
326 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEPoints,
327 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetDivAtDivEPoints,
328 const typename BasisType::ScalarViewType targetEPoints,
329 const typename BasisType::ScalarViewType targetDivEPoints,
330 const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
331 const BasisType* cellBasis,
332 ProjectionStruct<SpT, typename BasisType::scalarType> * projStruct){
333 typedef typename BasisType::scalarType scalarType;
334 typedef Kokkos::DynRankView<scalarType,SpT> ScalarViewType;
335 typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
336 const auto cellTopo = cellBasis->getBaseCellTopology();
337 ordinal_type dim = cellTopo.getDimension();
338 ordinal_type numTotalEvaluationPoints(targetAtEPoints.extent(1)),
339 numTotalDivEvaluationPoints(targetDivAtDivEPoints.extent(1));
340 ordinal_type basisCardinality = cellBasis->getCardinality();
341 ordinal_type numCells = targetAtEPoints.extent(0);
342 const ordinal_type sideDim = dim-1;
344 const std::string& name = cellBasis->getName();
346 ordinal_type numSides = cellBasis->getBaseCellTopology().getSideCount();
347 ScalarViewType refSideNormal(
"refSideNormal", dim);
349 ordinal_type numSideDofs(0);
350 for(ordinal_type is=0; is<numSides; ++is)
351 numSideDofs += cellBasis->getDofCount(sideDim,is);
353 Kokkos::View<ordinal_type*, SpT> computedDofs(
"computedDofs",numSideDofs);
355 const Kokkos::RangePolicy<SpT> policy(0, numCells);
357 auto targetEPointsRange = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->getTargetPointsRange());
358 auto basisEPointsRange = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->getBasisPointsRange());
360 auto tagToOrdinal = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(), cellBasis->getAllDofOrdinal());
362 ordinal_type numTotalBasisEPoints = projStruct->getNumBasisEvalPoints(), numTotalBasisDivEPoints = projStruct->getNumBasisDerivEvalPoints();
363 ScalarViewType basisEPoints_(
"basisEPoints",numCells,numTotalBasisEPoints, dim);
364 ScalarViewType basisDivEPoints(
"basisDivEPoints",numCells,numTotalBasisDivEPoints, dim);
365 getHDivEvaluationPoints(basisEPoints_, basisDivEPoints, orts, cellBasis, projStruct, EvalPointsType::BASIS);
367 ScalarViewType basisAtBasisEPoints(
"basisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, dim);
368 ScalarViewType basisAtTargetEPoints(
"basisAtTargetEPoints",numCells,basisCardinality, numTotalEvaluationPoints, dim);
370 ScalarViewType nonOrientedBasisAtBasisEPoints(
"nonOrientedBasisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, dim);
371 ScalarViewType nonOrientedBasisAtTargetEPoints(
"nonOrientedBasisAtTargetEPoints",numCells,basisCardinality, numTotalEvaluationPoints, dim);
372 for(ordinal_type ic=0; ic<numCells; ++ic) {
373 cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtTargetEPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(targetEPoints, ic, Kokkos::ALL(), Kokkos::ALL()));
374 cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtBasisEPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(basisEPoints_, ic, Kokkos::ALL(), Kokkos::ALL()));
381 ScalarViewType basisDivAtBasisDivEPoints;
382 ScalarViewType basisDivAtTargetDivEPoints;
383 if(numTotalDivEvaluationPoints>0) {
384 ScalarViewType nonOrientedBasisDivAtTargetDivEPoints, nonOrientedBasisDivAtBasisDivEPoints;
385 basisDivAtBasisDivEPoints = ScalarViewType (
"basisDivAtBasisDivEPoints",numCells,basisCardinality, numTotalBasisDivEPoints);
386 nonOrientedBasisDivAtBasisDivEPoints = ScalarViewType (
"nonOrientedBasisDivAtBasisDivEPoints",numCells,basisCardinality, numTotalBasisDivEPoints);
387 basisDivAtTargetDivEPoints = ScalarViewType(
"basisDivAtTargetDivEPoints",numCells,basisCardinality, numTotalDivEvaluationPoints);
388 nonOrientedBasisDivAtTargetDivEPoints = ScalarViewType(
"nonOrientedBasisDivAtTargetDivEPoints",numCells,basisCardinality, numTotalDivEvaluationPoints);
390 for(ordinal_type ic=0; ic<numCells; ++ic) {
391 cellBasis->getValues(Kokkos::subview(nonOrientedBasisDivAtBasisDivEPoints,ic,Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(basisDivEPoints, ic, Kokkos::ALL(), Kokkos::ALL()),OPERATOR_DIV);
392 cellBasis->getValues(Kokkos::subview(nonOrientedBasisDivAtTargetDivEPoints,ic,Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(targetDivEPoints, ic, Kokkos::ALL(), Kokkos::ALL()),OPERATOR_DIV);
398 ScalarViewType refSidesNormal(
"refSidesNormal", numSides, dim);
399 ordinal_type computedDofsCount = 0;
400 for(ordinal_type is=0; is<numSides; ++is) {
402 ordinal_type sideCardinality = cellBasis->getDofCount(sideDim,is);
404 ordinal_type numTargetEPoints = range_size(targetEPointsRange(sideDim,is));
405 ordinal_type numBasisEPoints = range_size(basisEPointsRange(sideDim,is));
407 for(ordinal_type i=0; i<sideCardinality; ++i)
408 computedDofs(computedDofsCount++) = cellBasis->getDofOrdinal(sideDim, is, i);
410 auto sideNormal = Kokkos::subview(refSidesNormal,is,Kokkos::ALL());
411 auto sideNormalHost = Kokkos::create_mirror_view(sideNormal);
413 Kokkos::deep_copy(sideNormal, sideNormalHost);
415 ScalarViewType basisNormalAtBasisEPoints(
"normalBasisAtBasisEPoints",numCells,sideCardinality, numBasisEPoints);
416 ScalarViewType wBasisNormalAtBasisEPoints(
"wBasisNormalAtBasisEPoints",numCells,sideCardinality, numBasisEPoints);
417 ScalarViewType wBasisNormalAtTargetEPoints(
"wBasisNormalAtTargetEPoints",numCells,sideCardinality, numTargetEPoints);
418 ScalarViewType targetNormalAtTargetEPoints(
"targetNormalAtTargetEPoints",numCells, numTargetEPoints);
420 ordinal_type offsetBasis = basisEPointsRange(sideDim,is).first;
421 ordinal_type offsetTarget = targetEPointsRange(sideDim,is).first;
422 auto targetEWeights = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->getTargetEvalWeights(sideDim,is));
423 auto basisEWeights = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->getBasisEvalWeights(sideDim,is));
426 typedef ComputeBasisCoeffsOnSides_HDiv<ScalarViewType, decltype(basisEWeights), decltype(tagToOrdinal)> functorTypeSide;
427 Kokkos::parallel_for(policy, functorTypeSide(basisNormalAtBasisEPoints, basisAtBasisEPoints,
428 basisEWeights, wBasisNormalAtBasisEPoints, targetEWeights,
429 basisAtTargetEPoints, wBasisNormalAtTargetEPoints, tagToOrdinal,
430 targetAtEPoints, targetNormalAtTargetEPoints,
431 refSidesNormal, sideCardinality, offsetBasis,
432 offsetTarget, sideDim,
435 ScalarViewType sideMassMat_(
"sideMassMat_", numCells, sideCardinality+1, sideCardinality+1),
436 sideRhsMat_(
"rhsMat_", numCells, sideCardinality+1);
438 ScalarViewType targetEWeights_(
"targetEWeights", numCells, 1, targetEWeights.extent(0));
441 range_type range_H(0, sideCardinality);
442 range_type range_B(sideCardinality, sideCardinality+1);
443 ScalarViewType ones(
"ones",numCells,1,numBasisEPoints);
444 Kokkos::deep_copy(ones,1);
452 typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, SpT> WorkArrayViewType;
453 ScalarViewType t_(
"t",numCells, sideCardinality+1);
454 WorkArrayViewType w_(
"w",numCells, sideCardinality+1);
456 auto sideDof = Kokkos::subview(tagToOrdinal, sideDim, is, Kokkos::ALL());
458 ElemSystem sideSystem(
"sideSystem",
false);
459 sideSystem.solve(basisCoeffs, sideMassMat_, sideRhsMat_, t_, w_, sideDof, sideCardinality, 1);
464 ordinal_type numCellDofs = cellBasis->getDofCount(dim,0);
468 Basis<SpT,scalarType,scalarType> *hcurlBasis = NULL;
469 if(cellTopo.getKey() == shards::getCellTopologyData<shards::Hexahedron<8> >()->key)
470 hcurlBasis =
new Basis_HCURL_HEX_In_FEM<SpT,scalarType,scalarType>(cellBasis->getDegree());
471 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Tetrahedron<4> >()->key)
472 hcurlBasis =
new Basis_HCURL_TET_In_FEM<SpT,scalarType,scalarType>(cellBasis->getDegree());
473 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Quadrilateral<4> >()->key)
474 hcurlBasis =
new Basis_HGRAD_QUAD_Cn_FEM<SpT,scalarType,scalarType>(cellBasis->getDegree());
475 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Triangle<3> >()->key)
476 hcurlBasis =
new Basis_HGRAD_TRI_Cn_FEM<SpT,scalarType,scalarType>(cellBasis->getDegree());
478 std::stringstream ss;
479 ss <<
">>> ERROR (Intrepid2::ProjectionTools::getHDivEvaluationPoints): "
480 <<
"Method not implemented for basis " << name;
481 INTREPID2_TEST_FOR_EXCEPTION(
true, std::runtime_error, ss.str().c_str() );
483 if(hcurlBasis == NULL)
return;
486 auto targetDivEPointsRange = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->getTargetDerivPointsRange());
487 auto basisDivEPointsRange = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->getBasisDerivPointsRange());
489 ordinal_type numTargetDivEPoints = range_size(targetDivEPointsRange(dim,0));
490 ordinal_type numBasisDivEPoints = range_size(basisDivEPointsRange(dim,0));
492 auto targetDivEWeights = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->getTargetDerivEvalWeights(dim,0));
493 auto divEWeights = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->getBasisDerivEvalWeights(dim,0));
495 ordinal_type offsetBasisDiv = basisDivEPointsRange(dim,0).first;
496 ordinal_type offsetTargetDiv = targetDivEPointsRange(dim,0).first;
498 ScalarViewType weightedBasisDivAtBasisEPoints(
"weightedBasisDivAtBasisEPoints",numCells,numCellDofs, numBasisDivEPoints);
499 ScalarViewType weightedBasisDivAtTargetEPoints(
"weightedBasisDivAtTargetEPoints",numCells, numCellDofs, numTargetDivEPoints);
500 ScalarViewType basisDivAtBasisEPoints(
"basisDivAtBasisEPoints",numCells,numCellDofs, numBasisDivEPoints);
501 ScalarViewType targetSideDivAtBasisEPoints(
"targetSideDivAtBasisEPoints",numCells, numBasisDivEPoints);
503 auto cellDofs = Kokkos::subview(tagToOrdinal, dim, 0, Kokkos::ALL());
504 typedef ComputeBasisCoeffsOnCells_HDiv<decltype(basisCoeffs), ScalarViewType, decltype(divEWeights), decltype(computedDofs), decltype(cellDofs)> functorType;
505 Kokkos::parallel_for(policy, functorType( basisCoeffs, targetSideDivAtBasisEPoints, basisDivAtBasisEPoints,
506 basisDivAtBasisDivEPoints, divEWeights, weightedBasisDivAtBasisEPoints, targetDivEWeights, basisDivAtTargetDivEPoints, weightedBasisDivAtTargetEPoints,
507 computedDofs, cellDofs, numCellDofs, offsetBasisDiv, offsetTargetDiv, numSideDofs));
510 ordinal_type hcurlBasisCardinality = hcurlBasis->getCardinality();
511 ordinal_type numCurlInteriorDOFs = hcurlBasis->getDofCount(dim,0);
513 range_type range_H(0, numCellDofs);
514 range_type range_B(numCellDofs, numCellDofs+numCurlInteriorDOFs);
517 ScalarViewType massMat_(
"massMat_",numCells,numCellDofs+numCurlInteriorDOFs,numCellDofs+numCurlInteriorDOFs),
518 rhsMatTrans(
"rhsMatTrans",numCells,numCellDofs+numCurlInteriorDOFs);
524 if(numCurlInteriorDOFs>0){
525 ordinal_type numTargetEPoints = range_size(targetEPointsRange(dim,0));
526 ordinal_type numBasisEPoints = range_size(basisEPointsRange(dim,0));
528 auto targetEWeights = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->getTargetEvalWeights(dim,0));
529 auto basisEWeights = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(),projStruct->getBasisEvalWeights(dim,0));
531 ordinal_type offsetBasis = basisEPointsRange(dim,0).first;
533 auto basisEPoints = Kokkos::subview(basisEPoints_, 0, basisEPointsRange(dim,0), Kokkos::ALL());
535 ScalarViewType negPartialProjAtBasisEPoints(
"targetSideAtBasisEPoints",numCells, numBasisEPoints, dim);
536 ScalarViewType nonWeightedBasisAtBasisEPoints(
"basisAtBasisEPoints",numCells,numCellDofs, numBasisEPoints, dim);
537 ScalarViewType hcurlBasisCurlAtBasisEPoints(
"hcurlBasisCurlAtBasisEPoints",hcurlBasisCardinality, numBasisEPoints,dim);
538 ScalarViewType hcurlBasisCurlAtTargetEPoints(
"hcurlBasisCurlAtTargetEPoints", hcurlBasisCardinality,numTargetEPoints, dim);
539 ScalarViewType wHcurlBasisCurlAtBasisEPoints(
"wHcurlBasisHcurlAtBasisEPoints", numCells, numCurlInteriorDOFs, numBasisEPoints,dim);
540 ScalarViewType wHcurlBasisCurlAtTargetEPoints(
"wHcurlBasisHcurlAtTargetEPoints",numCells, numCurlInteriorDOFs, numTargetEPoints,dim);
542 hcurlBasis->getValues(hcurlBasisCurlAtBasisEPoints, basisEPoints, OPERATOR_CURL);
543 hcurlBasis->getValues(hcurlBasisCurlAtTargetEPoints, Kokkos::subview(targetEPoints,0,targetEPointsRange(dim,0),Kokkos::ALL()), OPERATOR_CURL);
545 auto hCurlTagToOrdinal = Kokkos::create_mirror_view_and_copy(
typename SpT::memory_space(), hcurlBasis->getAllDofOrdinal());
546 auto cellHCurlDof = Kokkos::subview(hCurlTagToOrdinal, dim, 0, range_type(0, numCurlInteriorDOFs));
548 typedef ComputeHCurlBasisCoeffsOnCells_HDiv<decltype(basisCoeffs), ScalarViewType, decltype(divEWeights),
549 decltype(tagToOrdinal), decltype(computedDofs), decltype(cellDofs)> functorTypeHCurlCells;
550 Kokkos::parallel_for(policy, functorTypeHCurlCells(basisCoeffs, negPartialProjAtBasisEPoints, nonWeightedBasisAtBasisEPoints,
551 basisAtBasisEPoints, hcurlBasisCurlAtBasisEPoints, basisEWeights, wHcurlBasisCurlAtBasisEPoints, targetEWeights,
552 hcurlBasisCurlAtTargetEPoints, wHcurlBasisCurlAtTargetEPoints, tagToOrdinal, computedDofs, cellHCurlDof,
553 numCellDofs, offsetBasis, numSideDofs, dim));
556 FunctionSpaceTools<SpT >::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_B), Kokkos::subview(targetAtEPoints, Kokkos::ALL(), targetEPointsRange(dim,0), Kokkos::ALL()), wHcurlBasisCurlAtTargetEPoints);
561 typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, SpT> WorkArrayViewType;
562 ScalarViewType t_(
"t",numCells, numCellDofs+numCurlInteriorDOFs);
563 WorkArrayViewType w_(
"w",numCells, numCellDofs+numCurlInteriorDOFs);
565 ElemSystem cellSystem(
"cellSystem",
true);
566 cellSystem.solve(basisCoeffs, massMat_, rhsMatTrans, t_, w_, cellDofs, numCellDofs, numCurlInteriorDOFs);
Header file for the abstract base class Intrepid2::DefaultCubatureFactory.
const range_tag getPointsRange(const EvalPointsType type) const
Returns the range tag of the basis/target evaluation points in subcells.
const range_tag getDerivPointsRange(const EvalPointsType type) const
Returns the range tag of the basis/target derivative evaluation points on subcells.
ordinal_type getMaxNumEvalPoints(const EvalPointsType type) const
Returns the maximum number of evaluation points across all the subcells.
view_type getDerivEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId, EvalPointsType type) const
Returns the evaluation points for basis/target derivatives on a subcell.
view_type getTargetEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the points where to evaluate the target function on a subcell.
An helper class to compute the evaluation points and weights needed for performing projections...
const key_tag getTopologyKey() const
Returns the key tag for subcells.
view_type getEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId, EvalPointsType type) const
Returns the basis/target evaluation points on a subcell.