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 ViewType4 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 ViewType4 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);
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);
318 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEPoints,
319 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetDivAtDivEPoints,
320 const typename BasisType::ScalarViewType targetEPoints,
321 const typename BasisType::ScalarViewType targetDivEPoints,
322 const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
323 const BasisType* cellBasis,
325 typedef typename BasisType::scalarType scalarType;
326 typedef Kokkos::DynRankView<scalarType,DeviceType> ScalarViewType;
327 typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
328 const auto cellTopo = cellBasis->getBaseCellTopology();
329 ordinal_type dim = cellTopo.getDimension();
330 ordinal_type numTotalEvaluationPoints(targetAtEPoints.extent(1)),
331 numTotalDivEvaluationPoints(targetDivAtDivEPoints.extent(1));
332 ordinal_type basisCardinality = cellBasis->getCardinality();
333 ordinal_type numCells = targetAtEPoints.extent(0);
334 const ordinal_type sideDim = dim-1;
336 const std::string& name = cellBasis->getName();
338 ordinal_type numSides = cellBasis->getBaseCellTopology().getSideCount();
339 ScalarViewType refSideNormal(
"refSideNormal", dim);
341 ordinal_type numSideDofs(0);
342 for(ordinal_type is=0; is<numSides; ++is)
343 numSideDofs += cellBasis->getDofCount(sideDim,is);
345 Kokkos::View<ordinal_type*, DeviceType> computedDofs(
"computedDofs",numSideDofs);
347 const Kokkos::RangePolicy<ExecSpaceType> policy(0, numCells);
352 auto tagToOrdinal = Kokkos::create_mirror_view_and_copy(MemSpaceType(), cellBasis->getAllDofOrdinal());
355 ScalarViewType basisEPoints_(
"basisEPoints",numCells,numTotalBasisEPoints, dim);
356 ScalarViewType basisDivEPoints(
"basisDivEPoints",numCells,numTotalBasisDivEPoints, dim);
357 getHDivEvaluationPoints(basisEPoints_, basisDivEPoints, orts, cellBasis, projStruct, EvalPointsType::BASIS);
359 ScalarViewType basisAtBasisEPoints(
"basisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, dim);
360 ScalarViewType basisAtTargetEPoints(
"basisAtTargetEPoints",numCells,basisCardinality, numTotalEvaluationPoints, dim);
362 ScalarViewType nonOrientedBasisAtBasisEPoints(
"nonOrientedBasisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, dim);
363 ScalarViewType nonOrientedBasisAtTargetEPoints(
"nonOrientedBasisAtTargetEPoints",numCells,basisCardinality, numTotalEvaluationPoints, dim);
364 for(ordinal_type ic=0; ic<numCells; ++ic) {
365 cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtTargetEPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(targetEPoints, ic, Kokkos::ALL(), Kokkos::ALL()));
366 cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtBasisEPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(basisEPoints_, ic, Kokkos::ALL(), Kokkos::ALL()));
373 ScalarViewType basisDivAtBasisDivEPoints;
374 ScalarViewType basisDivAtTargetDivEPoints;
375 if(numTotalDivEvaluationPoints>0) {
376 ScalarViewType nonOrientedBasisDivAtTargetDivEPoints, nonOrientedBasisDivAtBasisDivEPoints;
377 basisDivAtBasisDivEPoints = ScalarViewType (
"basisDivAtBasisDivEPoints",numCells,basisCardinality, numTotalBasisDivEPoints);
378 nonOrientedBasisDivAtBasisDivEPoints = ScalarViewType (
"nonOrientedBasisDivAtBasisDivEPoints",numCells,basisCardinality, numTotalBasisDivEPoints);
379 basisDivAtTargetDivEPoints = ScalarViewType(
"basisDivAtTargetDivEPoints",numCells,basisCardinality, numTotalDivEvaluationPoints);
380 nonOrientedBasisDivAtTargetDivEPoints = ScalarViewType(
"nonOrientedBasisDivAtTargetDivEPoints",numCells,basisCardinality, numTotalDivEvaluationPoints);
382 for(ordinal_type ic=0; ic<numCells; ++ic) {
383 cellBasis->getValues(Kokkos::subview(nonOrientedBasisDivAtBasisDivEPoints,ic,Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(basisDivEPoints, ic, Kokkos::ALL(), Kokkos::ALL()),OPERATOR_DIV);
384 cellBasis->getValues(Kokkos::subview(nonOrientedBasisDivAtTargetDivEPoints,ic,Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(targetDivEPoints, ic, Kokkos::ALL(), Kokkos::ALL()),OPERATOR_DIV);
390 ScalarViewType refSidesNormal(
"refSidesNormal", numSides, dim);
391 ordinal_type computedDofsCount = 0;
392 for(ordinal_type is=0; is<numSides; ++is) {
394 ordinal_type sideCardinality = cellBasis->getDofCount(sideDim,is);
396 ordinal_type numTargetEPoints = range_size(targetEPointsRange(sideDim,is));
397 ordinal_type numBasisEPoints = range_size(basisEPointsRange(sideDim,is));
399 auto sideDofs = Kokkos::subview(tagToOrdinal, sideDim, is, range_type(0,sideCardinality));
400 auto computedSideDofs = Kokkos::subview(computedDofs, range_type(computedDofsCount,computedDofsCount+sideCardinality));
401 deep_copy(computedSideDofs, sideDofs);
402 computedDofsCount += sideCardinality;
404 auto sideNormal = Kokkos::subview(refSidesNormal,is,Kokkos::ALL());
407 ScalarViewType basisNormalAtBasisEPoints(
"normalBasisAtBasisEPoints",numCells,sideCardinality, numBasisEPoints);
408 ScalarViewType wBasisNormalAtBasisEPoints(
"wBasisNormalAtBasisEPoints",numCells,sideCardinality, numBasisEPoints);
409 ScalarViewType wBasisNormalAtTargetEPoints(
"wBasisNormalAtTargetEPoints",numCells,sideCardinality, numTargetEPoints);
410 ScalarViewType targetNormalAtTargetEPoints(
"targetNormalAtTargetEPoints",numCells, numTargetEPoints);
412 ordinal_type offsetBasis = basisEPointsRange(sideDim,is).first;
413 ordinal_type offsetTarget = targetEPointsRange(sideDim,is).first;
414 auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->
getTargetEvalWeights(sideDim,is));
415 auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->
getBasisEvalWeights(sideDim,is));
418 typedef ComputeBasisCoeffsOnSides_HDiv<ScalarViewType,
decltype(basisEWeights),
decltype(tagToOrdinal),
decltype(targetAtEPoints)> functorTypeSide;
419 Kokkos::parallel_for(policy, functorTypeSide(basisNormalAtBasisEPoints, basisAtBasisEPoints,
420 basisEWeights, wBasisNormalAtBasisEPoints, targetEWeights,
421 basisAtTargetEPoints, wBasisNormalAtTargetEPoints, tagToOrdinal,
422 targetAtEPoints, targetNormalAtTargetEPoints,
423 refSidesNormal, sideCardinality, offsetBasis,
424 offsetTarget, sideDim,
427 ScalarViewType sideMassMat_(
"sideMassMat_", numCells, sideCardinality+1, sideCardinality+1),
428 sideRhsMat_(
"rhsMat_", numCells, sideCardinality+1);
430 ScalarViewType targetEWeights_(
"targetEWeights", numCells, 1, targetEWeights.extent(0));
433 range_type range_H(0, sideCardinality);
434 range_type range_B(sideCardinality, sideCardinality+1);
435 ScalarViewType ones(
"ones",numCells,1,numBasisEPoints);
436 Kokkos::deep_copy(ones,1);
444 typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
445 ScalarViewType t_(
"t",numCells, sideCardinality+1);
446 WorkArrayViewType w_(
"w",numCells, sideCardinality+1);
448 auto sideDof = Kokkos::subview(tagToOrdinal, sideDim, is, Kokkos::ALL());
451 sideSystem.
solve(basisCoeffs, sideMassMat_, sideRhsMat_, t_, w_, sideDof, sideCardinality, 1);
456 ordinal_type numCellDofs = cellBasis->getDofCount(dim,0);
461 if(cellTopo.getKey() == shards::getCellTopologyData<shards::Hexahedron<8> >()->key)
463 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Tetrahedron<4> >()->key)
465 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Quadrilateral<4> >()->key)
467 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Triangle<3> >()->key)
470 std::stringstream ss;
471 ss <<
">>> ERROR (Intrepid2::ProjectionTools::getHDivEvaluationPoints): "
472 <<
"Method not implemented for basis " << name;
473 INTREPID2_TEST_FOR_EXCEPTION(
true, std::runtime_error, ss.str().c_str() );
475 if(hcurlBasis == NULL)
return;
481 ordinal_type numTargetDivEPoints = range_size(targetDivEPointsRange(dim,0));
482 ordinal_type numBasisDivEPoints = range_size(basisDivEPointsRange(dim,0));
484 auto targetDivEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->
getTargetDerivEvalWeights(dim,0));
487 ordinal_type offsetBasisDiv = basisDivEPointsRange(dim,0).first;
488 ordinal_type offsetTargetDiv = targetDivEPointsRange(dim,0).first;
490 ScalarViewType weightedBasisDivAtBasisEPoints(
"weightedBasisDivAtBasisEPoints",numCells,numCellDofs, numBasisDivEPoints);
491 ScalarViewType weightedBasisDivAtTargetEPoints(
"weightedBasisDivAtTargetEPoints",numCells, numCellDofs, numTargetDivEPoints);
492 ScalarViewType basisDivAtBasisEPoints(
"basisDivAtBasisEPoints",numCells,numCellDofs, numBasisDivEPoints);
493 ScalarViewType targetSideDivAtBasisEPoints(
"targetSideDivAtBasisEPoints",numCells, numBasisDivEPoints);
495 auto cellDofs = Kokkos::subview(tagToOrdinal, dim, 0, Kokkos::ALL());
496 typedef ComputeBasisCoeffsOnCells_HDiv<
decltype(basisCoeffs), ScalarViewType,
decltype(divEWeights),
decltype(computedDofs),
decltype(cellDofs)> functorType;
497 Kokkos::parallel_for(policy, functorType( basisCoeffs, targetSideDivAtBasisEPoints, basisDivAtBasisEPoints,
498 basisDivAtBasisDivEPoints, divEWeights, weightedBasisDivAtBasisEPoints, targetDivEWeights, basisDivAtTargetDivEPoints, weightedBasisDivAtTargetEPoints,
499 computedDofs, cellDofs, numCellDofs, offsetBasisDiv, offsetTargetDiv, numSideDofs));
502 ordinal_type hcurlBasisCardinality = hcurlBasis->
getCardinality();
503 ordinal_type numCurlInteriorDOFs = hcurlBasis->
getDofCount(dim,0);
505 range_type range_H(0, numCellDofs);
506 range_type range_B(numCellDofs, numCellDofs+numCurlInteriorDOFs);
509 ScalarViewType massMat_(
"massMat_",numCells,numCellDofs+numCurlInteriorDOFs,numCellDofs+numCurlInteriorDOFs),
510 rhsMatTrans(
"rhsMatTrans",numCells,numCellDofs+numCurlInteriorDOFs);
516 if(numCurlInteriorDOFs>0){
517 ordinal_type numTargetEPoints = range_size(targetEPointsRange(dim,0));
518 ordinal_type numBasisEPoints = range_size(basisEPointsRange(dim,0));
520 auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->
getTargetEvalWeights(dim,0));
521 auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->
getBasisEvalWeights(dim,0));
523 ordinal_type offsetBasis = basisEPointsRange(dim,0).first;
525 auto basisEPoints = Kokkos::subview(basisEPoints_, 0, basisEPointsRange(dim,0), Kokkos::ALL());
527 ScalarViewType negPartialProjAtBasisEPoints(
"targetSideAtBasisEPoints",numCells, numBasisEPoints, dim);
528 ScalarViewType nonWeightedBasisAtBasisEPoints(
"basisAtBasisEPoints",numCells,numCellDofs, numBasisEPoints, dim);
529 ScalarViewType hcurlBasisCurlAtBasisEPoints(
"hcurlBasisCurlAtBasisEPoints",hcurlBasisCardinality, numBasisEPoints,dim);
530 ScalarViewType hcurlBasisCurlAtTargetEPoints(
"hcurlBasisCurlAtTargetEPoints", hcurlBasisCardinality,numTargetEPoints, dim);
531 ScalarViewType wHcurlBasisCurlAtBasisEPoints(
"wHcurlBasisHcurlAtBasisEPoints", numCells, numCurlInteriorDOFs, numBasisEPoints,dim);
532 ScalarViewType wHcurlBasisCurlAtTargetEPoints(
"wHcurlBasisHcurlAtTargetEPoints",numCells, numCurlInteriorDOFs, numTargetEPoints,dim);
534 hcurlBasis->
getValues(hcurlBasisCurlAtBasisEPoints, basisEPoints, OPERATOR_CURL);
535 hcurlBasis->
getValues(hcurlBasisCurlAtTargetEPoints, Kokkos::subview(targetEPoints,0,targetEPointsRange(dim,0),Kokkos::ALL()), OPERATOR_CURL);
537 auto hCurlTagToOrdinal = Kokkos::create_mirror_view_and_copy(MemSpaceType(), hcurlBasis->
getAllDofOrdinal());
538 auto cellHCurlDof = Kokkos::subview(hCurlTagToOrdinal, dim, 0, range_type(0, numCurlInteriorDOFs));
541 decltype(tagToOrdinal),
decltype(computedDofs),
decltype(cellDofs)> functorTypeHCurlCells;
542 Kokkos::parallel_for(policy, functorTypeHCurlCells(basisCoeffs, negPartialProjAtBasisEPoints, nonWeightedBasisAtBasisEPoints,
543 basisAtBasisEPoints, hcurlBasisCurlAtBasisEPoints, basisEWeights, wHcurlBasisCurlAtBasisEPoints, targetEWeights,
544 hcurlBasisCurlAtTargetEPoints, wHcurlBasisCurlAtTargetEPoints, tagToOrdinal, computedDofs, cellHCurlDof,
545 numCellDofs, offsetBasis, numSideDofs, dim));
548 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_B), Kokkos::subview(targetAtEPoints, Kokkos::ALL(), targetEPointsRange(dim,0), Kokkos::ALL()), wHcurlBasisCurlAtTargetEPoints);
553 typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
554 ScalarViewType t_(
"t",numCells, numCellDofs+numCurlInteriorDOFs);
555 WorkArrayViewType w_(
"w",numCells, numCellDofs+numCurlInteriorDOFs);
558 cellSystem.
solve(basisCoeffs, massMat_, rhsMatTrans, t_, w_, cellDofs, numCellDofs, numCurlInteriorDOFs);