From 75a8b216afa33c5db43cf5794e98ef1dfc22acc7 Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Wed, 6 Sep 2023 18:22:34 -0600 Subject: [PATCH] Intrepid2: Further improvements to Projections method Now we are fully working in the reference element and map back to the oriented frame at the end of the projection. This allows significant savings. --- .../Projection/Intrepid2_ProjectionTools.hpp | 2 +- .../Intrepid2_ProjectionToolsDefHCURL.hpp | 346 ++++++++--------- .../Intrepid2_ProjectionToolsDefHDIV.hpp | 335 ++++++++--------- .../Intrepid2_ProjectionToolsDefHGRAD.hpp | 255 ++++++------- .../Intrepid2_ProjectionToolsDefL2.hpp | 355 +++++++----------- 5 files changed, 569 insertions(+), 724 deletions(-) diff --git a/packages/intrepid2/src/Projection/Intrepid2_ProjectionTools.hpp b/packages/intrepid2/src/Projection/Intrepid2_ProjectionTools.hpp index 45856acc9edf..26f5d9b1d58f 100644 --- a/packages/intrepid2/src/Projection/Intrepid2_ProjectionTools.hpp +++ b/packages/intrepid2/src/Projection/Intrepid2_ProjectionTools.hpp @@ -369,7 +369,7 @@ class ProjectionTools { static void getHVolBasisCoeffs(Kokkos::DynRankView basisCoeffs, const Kokkos::DynRankView targetAtEvalPoints, - const Kokkos::DynRankView cellOrientations, + [[maybe_unused]] const Kokkos::DynRankView cellOrientations, const BasisType* cellBasis, ProjectionStruct * projStruct); diff --git a/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHCURL.hpp b/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHCURL.hpp index fd79d832d146..6d38c8139401 100644 --- a/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHCURL.hpp +++ b/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHCURL.hpp @@ -61,37 +61,31 @@ namespace FunctorsProjectionTools { template struct ComputeBasisCoeffsOnEdges_HCurl { - const ViewType1 basisTanAtBasisEPoints_; - const ViewType1 basisAtBasisEPoints_; - const ViewType2 basisEWeights_; - const ViewType1 wTanBasisAtBasisEPoints_; const ViewType2 targetEWeights_; const ViewType1 basisAtTargetEPoints_; const ViewType1 wTanBasisAtTargetEPoints_; const ViewType3 tagToOrdinal_; const ViewType4 targetAtTargetEPoints_; const ViewType1 targetTanAtTargetEPoints_; - const ViewType1 refEdgesTangent_; + const ViewType1 refEdgeTangent_; ordinal_type edgeCardinality_; - ordinal_type offsetBasis_; ordinal_type offsetTarget_; ordinal_type edgeDim_; ordinal_type dim_; ordinal_type iedge_; - ComputeBasisCoeffsOnEdges_HCurl(const ViewType1 basisTanAtBasisEPoints, - const ViewType1 basisAtBasisEPoints, const ViewType2 basisEWeights, const ViewType1 wTanBasisAtBasisEPoints, const ViewType2 targetEWeights, + ComputeBasisCoeffsOnEdges_HCurl( + const ViewType2 targetEWeights, const ViewType1 basisAtTargetEPoints, const ViewType1 wTanBasisAtTargetEPoints, const ViewType3 tagToOrdinal, const ViewType4 targetAtTargetEPoints, const ViewType1 targetTanAtTargetEPoints, - const ViewType1 refEdgesTangent, ordinal_type edgeCardinality, ordinal_type offsetBasis, + const ViewType1 refEdgeTangent, ordinal_type edgeCardinality, ordinal_type offsetTarget, ordinal_type edgeDim, ordinal_type dim, ordinal_type iedge) : - basisTanAtBasisEPoints_(basisTanAtBasisEPoints), - basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wTanBasisAtBasisEPoints_(wTanBasisAtBasisEPoints), targetEWeights_(targetEWeights), + targetEWeights_(targetEWeights), basisAtTargetEPoints_(basisAtTargetEPoints), wTanBasisAtTargetEPoints_(wTanBasisAtTargetEPoints), tagToOrdinal_(tagToOrdinal), targetAtTargetEPoints_(targetAtTargetEPoints), targetTanAtTargetEPoints_(targetTanAtTargetEPoints), - refEdgesTangent_(refEdgesTangent), edgeCardinality_(edgeCardinality), offsetBasis_(offsetBasis), + refEdgeTangent_(refEdgeTangent), edgeCardinality_(edgeCardinality), offsetTarget_(offsetTarget), edgeDim_(edgeDim), dim_(dim), iedge_(iedge) {} @@ -99,61 +93,53 @@ struct ComputeBasisCoeffsOnEdges_HCurl { KOKKOS_INLINE_FUNCTION operator()(const ordinal_type ic) const { - ordinal_type numBasisEPoints = basisEWeights_.extent(0); ordinal_type numTargetEPoints = targetEWeights_.extent(0); + typename ViewType1::value_type tmp = 0; for(ordinal_type j=0; j +typename ViewType5> struct ComputeBasisCoeffsOnFaces_HCurl { const ViewType1 basisCoeffs_; - const ViewType2 orts_; - const ViewType3 negPartialProjTan_; - const ViewType3 negPartialProjCurlNormal_; - const ViewType3 hgradBasisGradAtBasisEPoints_; - const ViewType3 wHgradBasisGradAtBasisEPoints_; - const ViewType3 basisCurlAtBasisCurlEPoints_; - const ViewType3 basisCurlNormalAtBasisCurlEPoints_; - const ViewType3 basisAtBasisEPoints_; - const ViewType3 normalTargetCurlAtTargetEPoints_; - const ViewType3 basisTanAtBasisEPoints_; - const ViewType3 hgradBasisGradAtTargetEPoints_; - const ViewType3 wHgradBasisGradAtTargetEPoints_; - const ViewType3 wNormalBasisCurlAtBasisCurlEPoints_; - const ViewType3 basisCurlAtTargetCurlEPoints_; - const ViewType3 wNormalBasisCurlBasisAtTargetCurlEPoints_; - const ViewType4 targetAtTargetEPoints_; - const ViewType3 targetTanAtTargetEPoints_; - const ViewType4 targetCurlAtTargetCurlEPoints_; - const ViewType5 basisEWeights_; - const ViewType5 targetEWeights_; - const ViewType5 basisCurlEWeights_; - const ViewType5 targetCurlEWeights_; - const ViewType6 tagToOrdinal_; - const ViewType6 hGradTagToOrdinal_; - const ViewType7 faceParametrization_; - const ViewType8 computedDofs_; - unsigned refTopologyKey_; + const ViewType1 negPartialProjTan_; + const ViewType1 negPartialProjCurlNormal_; + const ViewType1 hgradBasisGradAtBasisEPoints_; + const ViewType1 wHgradBasisGradAtBasisEPoints_; + const ViewType1 basisCurlAtBasisCurlEPoints_; + const ViewType1 basisCurlNormalAtBasisCurlEPoints_; + const ViewType1 basisAtBasisEPoints_; + const ViewType1 normalTargetCurlAtTargetEPoints_; + const ViewType1 basisTanAtBasisEPoints_; + const ViewType1 hgradBasisGradAtTargetEPoints_; + const ViewType1 wHgradBasisGradAtTargetEPoints_; + const ViewType1 wNormalBasisCurlAtBasisCurlEPoints_; + const ViewType1 basisCurlAtTargetCurlEPoints_; + const ViewType1 wNormalBasisCurlBasisAtTargetCurlEPoints_; + const ViewType2 targetAtTargetEPoints_; + const ViewType1 targetTanAtTargetEPoints_; + const ViewType2 targetCurlAtTargetCurlEPoints_; + const ViewType3 basisEWeights_; + const ViewType3 targetEWeights_; + const ViewType3 basisCurlEWeights_; + const ViewType3 targetCurlEWeights_; + const ViewType4 tagToOrdinal_; + const ViewType4 hGradTagToOrdinal_; + const ViewType5 computedDofs_; + const ViewType1 refFaceNormal_; ordinal_type offsetBasis_; ordinal_type offsetBasisCurl_; ordinal_type offsetTarget_; @@ -170,27 +156,27 @@ struct ComputeBasisCoeffsOnFaces_HCurl { ComputeBasisCoeffsOnFaces_HCurl(const ViewType1 basisCoeffs, - const ViewType2 orts, const ViewType3 negPartialProjTan, const ViewType3 negPartialProjCurlNormal, - const ViewType3 hgradBasisGradAtBasisEPoints, const ViewType3 wHgradBasisGradAtBasisEPoints, - const ViewType3 basisCurlAtBasisCurlEPoints, const ViewType3 basisCurlNormalAtBasisCurlEPoints, - const ViewType3 basisAtBasisEPoints, - const ViewType3 normalTargetCurlAtTargetEPoints, - const ViewType3 basisTanAtBasisEPoints, - const ViewType3 hgradBasisGradAtTargetEPoints, const ViewType3 wHgradBasisGradAtTargetEPoints, - const ViewType3 wNormalBasisCurlAtBasisCurlEPoints, const ViewType3 basisCurlAtTargetCurlEPoints, - const ViewType3 wNormalBasisCurlBasisAtTargetCurlEPoints, const ViewType4 targetAtTargetEPoints, - const ViewType3 targetTanAtTargetEPoints, const ViewType4 targetCurlAtTargetCurlEPoints, - const ViewType5 basisEWeights, const ViewType5 targetEWeights, - const ViewType5 basisCurlEWeights, const ViewType5 targetCurlEWeights, const ViewType6 tagToOrdinal, - const ViewType6 hGradTagToOrdinal, const ViewType7 faceParametrization, - const ViewType8 computedDofs, unsigned refTopologyKey, ordinal_type offsetBasis, + const ViewType1 negPartialProjTan, const ViewType1 negPartialProjCurlNormal, + const ViewType1 hgradBasisGradAtBasisEPoints, const ViewType1 wHgradBasisGradAtBasisEPoints, + const ViewType1 basisCurlAtBasisCurlEPoints, const ViewType1 basisCurlNormalAtBasisCurlEPoints, + const ViewType1 basisAtBasisEPoints, + const ViewType1 normalTargetCurlAtTargetEPoints, + const ViewType1 basisTanAtBasisEPoints, + const ViewType1 hgradBasisGradAtTargetEPoints, const ViewType1 wHgradBasisGradAtTargetEPoints, + const ViewType1 wNormalBasisCurlAtBasisCurlEPoints, const ViewType1 basisCurlAtTargetCurlEPoints, + const ViewType1 wNormalBasisCurlBasisAtTargetCurlEPoints, const ViewType2 targetAtTargetEPoints, + const ViewType1 targetTanAtTargetEPoints, const ViewType2 targetCurlAtTargetCurlEPoints, + const ViewType3 basisEWeights, const ViewType3 targetEWeights, + const ViewType3 basisCurlEWeights, const ViewType3 targetCurlEWeights, const ViewType4 tagToOrdinal, + const ViewType4 hGradTagToOrdinal, const ViewType5 computedDofs, + const ViewType1 refFaceNormal , ordinal_type offsetBasis, ordinal_type offsetBasisCurl, ordinal_type offsetTarget, ordinal_type offsetTargetCurl, ordinal_type iface, ordinal_type hgradCardinality, ordinal_type numFaces, ordinal_type numFaceDofs, ordinal_type numEdgeDofs, ordinal_type faceDim, ordinal_type dim): basisCoeffs_(basisCoeffs), - orts_(orts), negPartialProjTan_(negPartialProjTan), negPartialProjCurlNormal_(negPartialProjCurlNormal), + negPartialProjTan_(negPartialProjTan), negPartialProjCurlNormal_(negPartialProjCurlNormal), hgradBasisGradAtBasisEPoints_(hgradBasisGradAtBasisEPoints), wHgradBasisGradAtBasisEPoints_(wHgradBasisGradAtBasisEPoints), basisCurlAtBasisCurlEPoints_(basisCurlAtBasisCurlEPoints), basisCurlNormalAtBasisCurlEPoints_(basisCurlNormalAtBasisCurlEPoints), basisAtBasisEPoints_(basisAtBasisEPoints), @@ -201,8 +187,8 @@ struct ComputeBasisCoeffsOnFaces_HCurl { targetTanAtTargetEPoints_(targetTanAtTargetEPoints), targetCurlAtTargetCurlEPoints_(targetCurlAtTargetCurlEPoints), basisEWeights_(basisEWeights), targetEWeights_(targetEWeights), basisCurlEWeights_(basisCurlEWeights), targetCurlEWeights_(targetCurlEWeights), tagToOrdinal_(tagToOrdinal), - hGradTagToOrdinal_(hGradTagToOrdinal), faceParametrization_(faceParametrization), - computedDofs_(computedDofs), refTopologyKey_(refTopologyKey), offsetBasis_(offsetBasis), + hGradTagToOrdinal_(hGradTagToOrdinal), computedDofs_(computedDofs), + refFaceNormal_(refFaceNormal), offsetBasis_(offsetBasis), offsetBasisCurl_(offsetBasisCurl), offsetTarget_(offsetTarget), offsetTargetCurl_(offsetTargetCurl), iface_(iface), hgradCardinality_(hgradCardinality), numFaces_(numFaces), @@ -213,15 +199,8 @@ struct ComputeBasisCoeffsOnFaces_HCurl { KOKKOS_INLINE_FUNCTION operator()(const ordinal_type ic) const { - ordinal_type fOrt[6]; - orts_(ic).getFaceOrientation(fOrt, numFaces_); - - ordinal_type ort = fOrt[iface_]; - - typename ViewType3::value_type data[3*3]; - auto tangentsAndNormal = ViewType3(data, dim_, dim_); - Impl::OrientationTools::getRefSideTangentsAndNormal(tangentsAndNormal, faceParametrization_,refTopologyKey_, iface_, ort); - typename ViewType3::value_type n[3] = {tangentsAndNormal(2,0), tangentsAndNormal(2,1), tangentsAndNormal(2,2)}; + typename ViewType1::value_type n[3] = {refFaceNormal_(0), refFaceNormal_(1), refFaceNormal_(2)}; + typename ViewType1::value_type tmp=0; ordinal_type numBasisEPoints = basisEWeights_.extent(0); ordinal_type numTargetEPoints = targetEWeights_.extent(0); @@ -253,22 +232,24 @@ struct ComputeBasisCoeffsOnFaces_HCurl { ordinal_type dp1 = (d+1) % dim_; ordinal_type dp2 = (d+2) % dim_; // basis \times n - basisTanAtBasisEPoints_(ic,j,iq,d) = basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,dp1)*n[dp2] - basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,dp2)*n[dp1]; + basisTanAtBasisEPoints_(0,j,iq,d) = basisAtBasisEPoints_(jdof,offsetBasis_+iq,dp1)*n[dp2] - basisAtBasisEPoints_(jdof,offsetBasis_+iq,dp2)*n[dp1]; } } //Note: we are not considering the jacobian of the orientation map for normals since it is simply a scalar term for the integrals and it does not affect the projection for(ordinal_type iq=0; iq +typename ViewType4> struct ComputeBasisCoeffsOnCell_HCurl { const ViewType1 basisCoeffs_; - const ViewType2 negPartialProj_; - const ViewType2 negPartialProjCurl_; - const ViewType2 cellBasisAtBasisEPoints_; - const ViewType2 cellBasisCurlAtBasisCurlEPoints_; - const ViewType2 basisAtBasisEPoints_; - const ViewType2 hgradBasisGradAtBasisEPoints_; - const ViewType2 basisCurlAtBasisCurlEPoints_; - const ViewType2 hgradBasisGradAtTargetEPoints_; - const ViewType2 basisCurlAtTargetCurlEPoints_; - const ViewType3 basisEWeights_; - const ViewType3 basisCurlEWeights_; - const ViewType2 wHgradBasisGradAtBasisEPoints_; - const ViewType2 wBasisCurlAtBasisCurlEPoints_; - const ViewType3 targetEWeights_; - const ViewType3 targetCurlEWeights_; - const ViewType2 wHgradBasisGradAtTargetEPoints_; - const ViewType2 wBasisCurlAtTargetCurlEPoints_; - const ViewType4 computedDofs_; - const ViewType5 tagToOrdinal_; - const ViewType5 hGradTagToOrdinal_; + const ViewType1 negPartialProj_; + const ViewType1 negPartialProjCurl_; + const ViewType1 cellBasisAtBasisEPoints_; + const ViewType1 cellBasisCurlAtBasisCurlEPoints_; + const ViewType1 basisAtBasisEPoints_; + const ViewType1 hgradBasisGradAtBasisEPoints_; + const ViewType1 basisCurlAtBasisCurlEPoints_; + const ViewType1 hgradBasisGradAtTargetEPoints_; + const ViewType1 basisCurlAtTargetCurlEPoints_; + const ViewType2 basisEWeights_; + const ViewType2 basisCurlEWeights_; + const ViewType1 wHgradBasisGradAtBasisEPoints_; + const ViewType1 wBasisCurlAtBasisCurlEPoints_; + const ViewType2 targetEWeights_; + const ViewType2 targetCurlEWeights_; + const ViewType1 wHgradBasisGradAtTargetEPoints_; + const ViewType1 wBasisCurlAtTargetCurlEPoints_; + const ViewType3 computedDofs_; + const ViewType4 tagToOrdinal_; + const ViewType4 hGradTagToOrdinal_; ordinal_type numCellDofs_; ordinal_type hgradCardinality_; ordinal_type offsetBasis_; @@ -336,16 +317,16 @@ struct ComputeBasisCoeffsOnCell_HCurl { ordinal_type dim_; ordinal_type derDim_; - ComputeBasisCoeffsOnCell_HCurl(const ViewType1 basisCoeffs, ViewType2 negPartialProj, ViewType2 negPartialProjCurl, - const ViewType2 cellBasisAtBasisEPoints, const ViewType2 cellBasisCurlAtBasisCurlEPoints, - const ViewType2 basisAtBasisEPoints, const ViewType2 hgradBasisGradAtBasisEPoints, const ViewType2 basisCurlAtBasisCurlEPoints, - const ViewType2 hgradBasisGradAtTargetEPoints, const ViewType2 basisCurlAtTargetCurlEPoints, - const ViewType3 basisEWeights, const ViewType3 basisCurlEWeights, - const ViewType2 wHgradBasisGradAtBasisEPoints, const ViewType2 wBasisCurlAtBasisCurlEPoints, - const ViewType3 targetEWeights, const ViewType3 targetCurlEWeights, - const ViewType2 wHgradBasisGradAtTargetEPoints, - const ViewType2 wBasisCurlAtTargetCurlEPoints, const ViewType4 computedDofs, - const ViewType5 tagToOrdinal, const ViewType5 hGradTagToOrdinal, + ComputeBasisCoeffsOnCell_HCurl(const ViewType1 basisCoeffs, ViewType1 negPartialProj, ViewType1 negPartialProjCurl, + const ViewType1 cellBasisAtBasisEPoints, const ViewType1 cellBasisCurlAtBasisCurlEPoints, + const ViewType1 basisAtBasisEPoints, const ViewType1 hgradBasisGradAtBasisEPoints, const ViewType1 basisCurlAtBasisCurlEPoints, + const ViewType1 hgradBasisGradAtTargetEPoints, const ViewType1 basisCurlAtTargetCurlEPoints, + const ViewType2 basisEWeights, const ViewType2 basisCurlEWeights, + const ViewType1 wHgradBasisGradAtBasisEPoints, const ViewType1 wBasisCurlAtBasisCurlEPoints, + const ViewType2 targetEWeights, const ViewType2 targetCurlEWeights, + const ViewType1 wHgradBasisGradAtTargetEPoints, + const ViewType1 wBasisCurlAtTargetCurlEPoints, const ViewType3 computedDofs, + const ViewType4 tagToOrdinal, const ViewType4 hGradTagToOrdinal, ordinal_type numCellDofs, ordinal_type hgradCardinality, ordinal_type offsetBasis, ordinal_type offsetBasisCurl, ordinal_type offsetTargetCurl, ordinal_type numEdgeFaceDofs, ordinal_type dim, ordinal_type derDim) : @@ -387,25 +368,25 @@ struct ComputeBasisCoeffsOnCell_HCurl { ordinal_type idof = tagToOrdinal_(dim_, 0, j); for(ordinal_type d=0; d ::getHCurlBasisCoeffs(Kokkos::DynRankViewgetDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0; ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount() : 0; - ScalarViewType refEdgesTangent("refEdgesTangent", numEdges, dim); - ScalarViewType refFacesTangents("refFaceTangents", numFaces, dim, 2); - ScalarViewType refFacesNormal("refFaceNormal", numFaces, dim); + ScalarViewType refEdgeTangent("refEdgeTangent", dim); + ScalarViewType refFaceNormal("refFaceNormal", dim); ordinal_type numEdgeDofs(0); for(ordinal_type ie=0; ie::getHCurlBasisCoeffs(Kokkos::DynRankViewgetBasisPointsRange(); auto basisCurlEPointsRange = projStruct->getBasisDerivPointsRange(); - auto refTopologyKey = projStruct->getTopologyKey(); - ordinal_type numTotalBasisEPoints = projStruct->getNumBasisEvalPoints(), numTotalBasisCurlEPoints = projStruct->getNumBasisDerivEvalPoints(); auto basisEPoints = projStruct->getAllEvalPoints(EvalPointsType::BASIS); auto basisCurlEPoints = projStruct->getAllDerivEvalPoints(EvalPointsType::BASIS); @@ -476,40 +454,29 @@ ProjectionTools::getHCurlBasisCoeffs(Kokkos::DynRankViewgetAllEvalPoints(EvalPointsType::TARGET); auto targetCurlEPoints = projStruct->getAllDerivEvalPoints(EvalPointsType::TARGET); - ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, dim); - ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints, dim); - { - ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtEPoints",basisCardinality, numTotalBasisEPoints, dim); - ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",basisCardinality, numTotalTargetEPoints, dim); - cellBasis->getValues(nonOrientedBasisAtTargetEPoints, targetEPoints); - cellBasis->getValues(nonOrientedBasisAtBasisEPoints, basisEPoints); - - OrientationTools::modifyBasisByOrientation(basisAtBasisEPoints, nonOrientedBasisAtBasisEPoints, orts, cellBasis); - OrientationTools::modifyBasisByOrientation(basisAtTargetEPoints, nonOrientedBasisAtTargetEPoints, orts, cellBasis); - } + ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",basisCardinality, numTotalBasisEPoints, dim); + ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",basisCardinality, numTotalTargetEPoints, dim); + cellBasis->getValues(basisAtTargetEPoints, targetEPoints); + cellBasis->getValues(basisAtBasisEPoints, basisEPoints); ScalarViewType basisCurlAtBasisCurlEPoints; ScalarViewType basisCurlAtTargetCurlEPoints; if(numTotalBasisCurlEPoints>0) { ScalarViewType nonOrientedBasisCurlAtTargetCurlEPoints, nonOrientedBasisCurlAtBasisCurlEPoints; if (dim == 3) { - basisCurlAtBasisCurlEPoints = ScalarViewType ("basisCurlAtBasisCurlEPoints",numCells,basisCardinality, numTotalBasisCurlEPoints, dim); - nonOrientedBasisCurlAtBasisCurlEPoints = ScalarViewType ("nonOrientedBasisCurlAtBasisCurlEPoints", basisCardinality, numTotalBasisCurlEPoints, dim); - basisCurlAtTargetCurlEPoints = ScalarViewType("basisCurlAtTargetCurlEPoints",numCells,basisCardinality, numTotalTargetCurlEPoints, dim); - nonOrientedBasisCurlAtTargetCurlEPoints = ScalarViewType("nonOrientedBasisCurlAtTargetCurlEPoints",basisCardinality, numTotalTargetCurlEPoints, dim); + basisCurlAtBasisCurlEPoints = ScalarViewType ("basisCurlAtBasisCurlEPoints",basisCardinality, numTotalBasisCurlEPoints, dim); + basisCurlAtTargetCurlEPoints = ScalarViewType("basisCurlAtTargetCurlEPoints",basisCardinality, numTotalTargetCurlEPoints, dim); } else { - basisCurlAtBasisCurlEPoints = ScalarViewType ("basisCurlAtBasisCurlEPoints",numCells,basisCardinality, numTotalBasisCurlEPoints); - nonOrientedBasisCurlAtBasisCurlEPoints = ScalarViewType ("nonOrientedBasisCurlAtBasisCurlEPoints",basisCardinality, numTotalBasisCurlEPoints); - basisCurlAtTargetCurlEPoints = ScalarViewType("basisCurlAtTargetCurlEPoints",numCells,basisCardinality, numTotalTargetCurlEPoints); - nonOrientedBasisCurlAtTargetCurlEPoints = ScalarViewType("nonOrientedBasisCurlAtTargetCurlEPoints",basisCardinality, numTotalTargetCurlEPoints); + basisCurlAtBasisCurlEPoints = ScalarViewType ("basisCurlAtBasisCurlEPoints",basisCardinality, numTotalBasisCurlEPoints); + basisCurlAtTargetCurlEPoints = ScalarViewType("basisCurlAtTargetCurlEPoints",basisCardinality, numTotalTargetCurlEPoints); } - cellBasis->getValues(nonOrientedBasisCurlAtBasisCurlEPoints, basisCurlEPoints,OPERATOR_CURL); - cellBasis->getValues(nonOrientedBasisCurlAtTargetCurlEPoints, targetCurlEPoints,OPERATOR_CURL); - OrientationTools::modifyBasisByOrientation(basisCurlAtBasisCurlEPoints, nonOrientedBasisCurlAtBasisCurlEPoints, orts, cellBasis); - OrientationTools::modifyBasisByOrientation(basisCurlAtTargetCurlEPoints, nonOrientedBasisCurlAtTargetCurlEPoints, orts, cellBasis); + cellBasis->getValues(basisCurlAtBasisCurlEPoints, basisCurlEPoints,OPERATOR_CURL); + cellBasis->getValues(basisCurlAtTargetCurlEPoints, targetCurlEPoints,OPERATOR_CURL); } + ScalarViewType refBasisCoeffs("refBasisCoeffs", basisCoeffs.extent(0), basisCoeffs.extent(1)); + ordinal_type computedDofsCount = 0; for(ordinal_type ie=0; ie::getHCurlBasisCoeffs(Kokkos::DynRankView::getReferenceEdgeTangent(refEdgeTan, ie, cellTopo); - } + CellTools::getReferenceEdgeTangent(refEdgeTangent, ie, cellTopo); - ScalarViewType basisTanAtBasisEPoints("basisTanAtBasisEPoints",numCells,edgeCardinality, numBasisEPoints); - ScalarViewType weightedTanBasisAtBasisEPoints("weightedTanBasisAtBasisEPoints",numCells,edgeCardinality, numBasisEPoints); + ScalarViewType basisTanAtBasisEPoints("basisTanAtBasisEPoints",1,edgeCardinality, numBasisEPoints); + ScalarViewType weightedTanBasisAtBasisEPoints("weightedTanBasisAtBasisEPoints",1,edgeCardinality, numBasisEPoints); ScalarViewType weightedTanBasisAtTargetEPoints("weightedTanBasisAtTargetEPoints",numCells,edgeCardinality, numTargetEPoints); ScalarViewType targetTanAtTargetEPoints("normalTargetAtTargetEPoints",numCells, numTargetEPoints); auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(edgeDim,ie)); auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(edgeDim,ie)); - //Note: we are not considering the jacobian of the orientation map since it is simply a scalar term for the integrals and it does not affect the projection ordinal_type offsetBasis = basisEPointsRange(edgeDim, ie).first; ordinal_type offsetTarget = targetEPointsRange(edgeDim, ie).first; - using functorTypeEdge = FunctorsProjectionTools::ComputeBasisCoeffsOnEdges_HCurl; - Kokkos::parallel_for(policy, functorTypeEdge(basisTanAtBasisEPoints,basisAtBasisEPoints,basisEWeights, - weightedTanBasisAtBasisEPoints, targetEWeights, + Kokkos::parallel_for(Kokkos::MDRangePolicy >({0,0}, {edgeCardinality,numBasisEPoints}), + KOKKOS_LAMBDA (const ordinal_type &j, const ordinal_type &iq) { + ordinal_type jdof = tagToOrdinal(edgeDim, ie, j); + for(ordinal_type d=0; d ; + Kokkos::parallel_for(policy, functorTypeEdge(targetEWeights, basisAtTargetEPoints, weightedTanBasisAtTargetEPoints, tagToOrdinal, targetAtTargetEPoints, targetTanAtTargetEPoints, - refEdgesTangent, edgeCardinality, offsetBasis, + refEdgeTangent, edgeCardinality, offsetTarget, edgeDim, dim, ie)); - ScalarViewType edgeMassMat_("edgeMassMat_", numCells, edgeCardinality+1, edgeCardinality+1), + ScalarViewType edgeMassMat_("edgeMassMat_", 1, edgeCardinality+1, edgeCardinality+1), edgeRhsMat_("rhsMat_", numCells, edgeCardinality+1); - ScalarViewType eWeights_("eWeights_", numCells, 1, basisEWeights.extent(0)), targetEWeights_("targetEWeights", numCells, 1, targetEWeights.extent(0)); + ScalarViewType eWeights_("eWeights_", 1, 1, basisEWeights.extent(0)), targetEWeights_("targetEWeights", numCells, 1, targetEWeights.extent(0)); RealSpaceTools::clone(eWeights_, basisEWeights); RealSpaceTools::clone(targetEWeights_, targetEWeights); @@ -562,8 +532,8 @@ ProjectionTools::getHCurlBasisCoeffs(Kokkos::DynRankView::getHCurlBasisCoeffs(Kokkos::DynRankViewgetAllDofOrdinal()); - ScalarViewType basisTanAtBasisEPoints("basisTanAtBasisEPoints",numCells,numFaceDofs, numBasisEPoints,dim); - ScalarViewType basisCurlNormalAtBasisCurlEPoints("normaBasisCurlAtBasisEPoints",numCells,numFaceDofs, numBasisCurlEPoints); + ScalarViewType basisTanAtBasisEPoints("basisTanAtBasisEPoints",1,numFaceDofs, numBasisEPoints,dim); + ScalarViewType basisCurlNormalAtBasisCurlEPoints("normaBasisCurlAtBasisEPoints",1,numFaceDofs, numBasisCurlEPoints); ScalarViewType wNormalBasisCurlAtBasisCurlEPoints("weightedNormalBasisCurlAtBasisEPoints",numCells,numFaceDofs, numBasisCurlEPoints); ScalarViewType targetTanAtTargetEPoints("targetTanAtTargetEPoints",numCells, numTargetEPoints, dim); @@ -633,11 +603,13 @@ ProjectionTools::getHCurlBasisCoeffs(Kokkos::DynRankViewgetTargetEvalWeights(faceDim,iface)); auto targetCurlEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetDerivEvalWeights(faceDim,iface)); auto basisCurlEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisDerivEvalWeights(faceDim,iface)); - const auto topoKey = refTopologyKey(faceDim, iface); - using functorTypeFaces = FunctorsProjectionTools::ComputeBasisCoeffsOnFaces_HCurl; - Kokkos::parallel_for(policy, functorTypeFaces(basisCoeffs, - orts, negPartialProjTan, negPartialProjCurlNormal, + + CellTools::getReferenceFaceNormal(refFaceNormal, iface, cellTopo); + + using functorTypeFaces = FunctorsProjectionTools::ComputeBasisCoeffsOnFaces_HCurl; + Kokkos::parallel_for(policy, functorTypeFaces(refBasisCoeffs, + negPartialProjTan, negPartialProjCurlNormal, hgradBasisGradAtBasisEPoints, wHgradBasisGradAtBasisEPoints, basisCurlAtBasisCurlEPoints, basisCurlNormalAtBasisCurlEPoints, basisAtBasisEPoints, @@ -648,8 +620,8 @@ ProjectionTools::getHCurlBasisCoeffs(Kokkos::DynRankView::getHCurlBasisCoeffs(Kokkos::DynRankView::integrate(Kokkos::subview(faceMassMat_,Kokkos::ALL(),range_H,range_H), basisCurlNormalAtBasisCurlEPoints, wNormalBasisCurlAtBasisCurlEPoints); - FunctionSpaceTools::integrate(Kokkos::subview(faceMassMat_,Kokkos::ALL(),range_H,range_B), basisTanAtBasisEPoints, wHgradBasisGradAtBasisEPoints); + FunctionSpaceTools::integrate(Kokkos::subview(faceMassMat_,Kokkos::ALL(),range_H,range_H), basisCurlNormalAtBasisCurlEPoints, Kokkos::subview(wNormalBasisCurlAtBasisCurlEPoints, std::make_pair(0,1), Kokkos::ALL(), Kokkos::ALL()) ); + FunctionSpaceTools::integrate(Kokkos::subview(faceMassMat_,Kokkos::ALL(),range_H,range_B), basisTanAtBasisEPoints, Kokkos::subview(wHgradBasisGradAtBasisEPoints, std::make_pair(0,1), Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()) ); FunctionSpaceTools::integrate(Kokkos::subview(faceRhsMat_,Kokkos::ALL(),range_H), normalTargetCurlAtTargetEPoints, wNormalBasisCurlBasisAtTargetCurlEPoints); FunctionSpaceTools::integrate(Kokkos::subview(faceRhsMat_,Kokkos::ALL(),range_H), negPartialProjCurlNormal, wNormalBasisCurlAtBasisCurlEPoints,true); @@ -676,8 +648,8 @@ ProjectionTools::getHCurlBasisCoeffs(Kokkos::DynRankView::getHCurlBasisCoeffs(Kokkos::DynRankViewgetValues(hgradBasisGradAtBasisEPoints,Kokkos::subview(basisEPoints, basisEPointsRange(dim, 0), Kokkos::ALL()), OPERATOR_GRAD); hgradBasis->getValues(hgradBasisGradAtTargetEPoints,Kokkos::subview(targetEPoints, targetEPointsRange(dim, 0), Kokkos::ALL()),OPERATOR_GRAD); - ScalarViewType cellBasisAtBasisEPoints("basisCellAtEPoints",numCells,numCellDofs, numBasisEPoints, dim); - ScalarViewType cellBasisCurlAtCurlEPoints("cellBasisCurlAtCurlEPoints",numCells,numCellDofs, numBasisCurlEPoints, derDim); + ScalarViewType cellBasisAtBasisEPoints("basisCellAtEPoints",1,numCellDofs, numBasisEPoints, dim); + ScalarViewType cellBasisCurlAtCurlEPoints("cellBasisCurlAtCurlEPoints",1,numCellDofs, numBasisCurlEPoints, derDim); ScalarViewType negPartialProjCurl("negPartialProjCurl", numCells, numBasisEPoints, derDim); ScalarViewType negPartialProj("negPartialProj", numCells, numBasisEPoints, dim); ScalarViewType wBasisCurlAtCurlEPoints("weightedBasisCurlAtBasisEPoints",numCells,numCellDofs, numBasisCurlEPoints,derDim); @@ -741,9 +713,9 @@ ProjectionTools::getHCurlBasisCoeffs(Kokkos::DynRankViewgetAllDofOrdinal()); - using functorTypeCell = FunctorsProjectionTools::ComputeBasisCoeffsOnCell_HCurl; - Kokkos::parallel_for(policy, functorTypeCell(basisCoeffs, negPartialProj, negPartialProjCurl, + Kokkos::parallel_for(policy, functorTypeCell(refBasisCoeffs, negPartialProj, negPartialProjCurl, cellBasisAtBasisEPoints, cellBasisCurlAtCurlEPoints, basisAtBasisEPoints, hgradBasisGradAtBasisEPoints, basisCurlAtBasisCurlEPoints, hgradBasisGradAtTargetEPoints, basisCurlAtTargetCurlEPoints, @@ -756,13 +728,13 @@ ProjectionTools::getHCurlBasisCoeffs(Kokkos::DynRankView::integrate(Kokkos::subview(cellMassMat_,Kokkos::ALL(),range_H,range_H), cellBasisCurlAtCurlEPoints, wBasisCurlAtCurlEPoints); - FunctionSpaceTools::integrate(Kokkos::subview(cellMassMat_,Kokkos::ALL(),range_H,range_B), cellBasisAtBasisEPoints, wHgradBasisGradAtBasisEPoints); + FunctionSpaceTools::integrate(Kokkos::subview(cellMassMat_,Kokkos::ALL(),range_H,range_H), cellBasisCurlAtCurlEPoints, Kokkos::subview(wBasisCurlAtCurlEPoints, std::make_pair(0,1), Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()) ); + FunctionSpaceTools::integrate(Kokkos::subview(cellMassMat_,Kokkos::ALL(),range_H,range_B), cellBasisAtBasisEPoints, Kokkos::subview(wHgradBasisGradAtBasisEPoints, std::make_pair(0,1), Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()) ); if(dim==3) FunctionSpaceTools::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_H), Kokkos::subview(targetCurlAtTargetCurlEPoints,Kokkos::ALL(),cellCurlPointsRange,Kokkos::ALL()), wBasisCurlBasisAtTargetCurlEPoints); else @@ -778,10 +750,12 @@ ProjectionTools::getHCurlBasisCoeffs(Kokkos::DynRankView::modifyBasisByOrientationInverse(basisCoeffs, refBasisCoeffs, orts, cellBasis, true); } } // Intrepid2 namespace diff --git a/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHDIV.hpp b/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHDIV.hpp index 30cff3ffbaaf..f3d07dc245e9 100644 --- a/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHDIV.hpp +++ b/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHDIV.hpp @@ -60,37 +60,30 @@ namespace FunctorsProjectionTools { template struct ComputeBasisCoeffsOnSides_HDiv { - const ViewType1 sideBasisNormalAtBasisEPoints_; - const ViewType1 basisAtBasisEPoints_; - const ViewType2 basisEWeights_; - const ViewType1 wBasisDofAtBasisEPoints_; const ViewType2 targetEWeights_; const ViewType1 basisAtTargetEPoints_; const ViewType1 wBasisDofAtTargetEPoints_; const ViewType3 tagToOrdinal_; const ViewType4 targetAtEPoints_; const ViewType1 targetAtTargetEPoints_; - const ViewType1 refSidesNormal_; + const ViewType1 refSideNormal_; ordinal_type sideCardinality_; - ordinal_type offsetBasis_; ordinal_type offsetTarget_; ordinal_type sideDim_; ordinal_type dim_; ordinal_type iside_; - ComputeBasisCoeffsOnSides_HDiv(const ViewType1 sideBasisNormalAtBasisEPoints, - const ViewType1 basisAtBasisEPoints, const ViewType2 basisEWeights, const ViewType1 wBasisDofAtBasisEPoints, const ViewType2 targetEWeights, + ComputeBasisCoeffsOnSides_HDiv(const ViewType2 targetEWeights, const ViewType1 basisAtTargetEPoints, const ViewType1 wBasisDofAtTargetEvalPoint, const ViewType3 tagToOrdinal, const ViewType4 targetAtEPoints, const ViewType1 targetAtTargetEPoints, - const ViewType1 refSidesNormal, ordinal_type sideCardinality, ordinal_type offsetBasis, + const ViewType1 refSideNormal, ordinal_type sideCardinality, ordinal_type offsetTarget, ordinal_type sideDim, ordinal_type dim, ordinal_type iside) : - sideBasisNormalAtBasisEPoints_(sideBasisNormalAtBasisEPoints), - basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisDofAtBasisEPoints_(wBasisDofAtBasisEPoints), targetEWeights_(targetEWeights), + targetEWeights_(targetEWeights), basisAtTargetEPoints_(basisAtTargetEPoints), wBasisDofAtTargetEPoints_(wBasisDofAtTargetEvalPoint), tagToOrdinal_(tagToOrdinal), targetAtEPoints_(targetAtEPoints), targetAtTargetEPoints_(targetAtTargetEPoints), - refSidesNormal_(refSidesNormal), sideCardinality_(sideCardinality), offsetBasis_(offsetBasis), + refSideNormal_(refSideNormal), sideCardinality_(sideCardinality), offsetTarget_(offsetTarget), sideDim_(sideDim), dim_(dim), iside_(iside) {} @@ -98,52 +91,46 @@ struct ComputeBasisCoeffsOnSides_HDiv { KOKKOS_INLINE_FUNCTION operator()(const ordinal_type ic) const { - //Note: we are not considering the jacobian of the orientation map since it is simply a scalar term for the integrals and it does not affect the projection + typename ViewType1::value_type tmp =0; for(ordinal_type j=0; j +typename ViewType4> struct ComputeBasisCoeffsOnCells_HDiv { const ViewType1 basisCoeffs_; - const ViewType2 negPartialProjAtBasisEPoints_; - const ViewType2 nonWeightedBasisAtBasisEPoints_; - const ViewType2 basisAtBasisEPoints_; - const ViewType3 basisEWeights_; - const ViewType2 wBasisAtBasisEPoints_; - const ViewType3 targetEWeights_; - const ViewType2 basisAtTargetEPoints_; - const ViewType2 wBasisAtTargetEPoints_; - const ViewType4 computedDofs_; - const ViewType5 cellDof_; + const ViewType1 negPartialProjAtBasisEPoints_; + const ViewType1 nonWeightedBasisAtBasisEPoints_; + const ViewType1 basisAtBasisEPoints_; + const ViewType2 basisEWeights_; + const ViewType1 wBasisAtBasisEPoints_; + const ViewType2 targetEWeights_; + const ViewType1 basisAtTargetEPoints_; + const ViewType1 wBasisAtTargetEPoints_; + const ViewType3 computedDofs_; + const ViewType4 cellDof_; ordinal_type numCellDofs_; ordinal_type offsetBasis_; ordinal_type offsetTarget_; ordinal_type numSideDofs_; - ComputeBasisCoeffsOnCells_HDiv(const ViewType1 basisCoeffs, ViewType2 negPartialProjAtBasisEPoints, const ViewType2 nonWeightedBasisAtBasisEPoints, - const ViewType2 basisAtBasisEPoints, const ViewType3 basisEWeights, const ViewType2 wBasisAtBasisEPoints, const ViewType3 targetEWeights, - const ViewType2 basisAtTargetEPoints, const ViewType2 wBasisAtTargetEPoints, const ViewType4 computedDofs, const ViewType5 cellDof, + ComputeBasisCoeffsOnCells_HDiv(const ViewType1 basisCoeffs, ViewType1 negPartialProjAtBasisEPoints, const ViewType1 nonWeightedBasisAtBasisEPoints, + const ViewType1 basisAtBasisEPoints, const ViewType2 basisEWeights, const ViewType1 wBasisAtBasisEPoints, const ViewType2 targetEWeights, + const ViewType1 basisAtTargetEPoints, const ViewType1 wBasisAtTargetEPoints, const ViewType3 computedDofs, const ViewType4 cellDof, ordinal_type numCellDofs, ordinal_type offsetBasis, ordinal_type offsetTarget, ordinal_type numSideDofs) : basisCoeffs_(basisCoeffs), negPartialProjAtBasisEPoints_(negPartialProjAtBasisEPoints), nonWeightedBasisAtBasisEPoints_(nonWeightedBasisAtBasisEPoints), basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisAtBasisEPoints_(wBasisAtBasisEPoints), targetEWeights_(targetEWeights), @@ -158,46 +145,46 @@ struct ComputeBasisCoeffsOnCells_HDiv { for(ordinal_type j=0; j +typename ViewType4, typename ViewType5> struct ComputeHCurlBasisCoeffsOnCells_HDiv { const ViewType1 basisCoeffs_; - const ViewType2 negPartialProjAtBasisEPoints_; - const ViewType2 nonWeightedBasisAtBasisEPoints_; - const ViewType2 basisAtBasisEPoints_; - const ViewType2 hcurlBasisCurlAtBasisEPoints_; - const ViewType3 basisEWeights_; - const ViewType2 wHCurlBasisAtBasisEPoints_; - const ViewType3 targetEWeights_; - const ViewType2 hcurlBasisCurlAtTargetEPoints_; - const ViewType2 wHCurlBasisAtTargetEPoints_; - const ViewType4 tagToOrdinal_; - const ViewType5 computedDofs_; - const ViewType6 hCurlDof_; + const ViewType1 negPartialProjAtBasisEPoints_; + const ViewType1 nonWeightedBasisAtBasisEPoints_; + const ViewType1 basisAtBasisEPoints_; + const ViewType1 hcurlBasisCurlAtBasisEPoints_; + const ViewType2 basisEWeights_; + const ViewType1 wHCurlBasisAtBasisEPoints_; + const ViewType2 targetEWeights_; + const ViewType1 hcurlBasisCurlAtTargetEPoints_; + const ViewType1 wHCurlBasisAtTargetEPoints_; + const ViewType3 tagToOrdinal_; + const ViewType4 computedDofs_; + const ViewType5 hCurlDof_; ordinal_type numCellDofs_; ordinal_type offsetBasis_; ordinal_type numSideDofs_; ordinal_type dim_; - ComputeHCurlBasisCoeffsOnCells_HDiv(const ViewType1 basisCoeffs, ViewType2 negPartialProjAtBasisEPoints, const ViewType2 nonWeightedBasisAtBasisEPoints, - const ViewType2 basisAtBasisEPoints, const ViewType2 hcurlBasisCurlAtBasisEPoints, const ViewType3 basisEWeights, const ViewType2 wHCurlBasisAtBasisEPoints, const ViewType3 targetEWeights, - const ViewType2 hcurlBasisCurlAtTargetEPoints, const ViewType2 wHCurlBasisAtTargetEPoints, const ViewType4 tagToOrdinal, const ViewType5 computedDofs, const ViewType6 hCurlDof, + ComputeHCurlBasisCoeffsOnCells_HDiv(const ViewType1 basisCoeffs, ViewType1 negPartialProjAtBasisEPoints, const ViewType1 nonWeightedBasisAtBasisEPoints, + const ViewType1 basisAtBasisEPoints, const ViewType1 hcurlBasisCurlAtBasisEPoints, const ViewType2 basisEWeights, const ViewType1 wHCurlBasisAtBasisEPoints, const ViewType2 targetEWeights, + const ViewType1 hcurlBasisCurlAtTargetEPoints, const ViewType1 wHCurlBasisAtTargetEPoints, const ViewType3 tagToOrdinal, const ViewType4 computedDofs, const ViewType5 hCurlDof, ordinal_type numCellDofs, ordinal_type offsetBasis, ordinal_type numSideDofs, ordinal_type dim) : basisCoeffs_(basisCoeffs), negPartialProjAtBasisEPoints_(negPartialProjAtBasisEPoints), nonWeightedBasisAtBasisEPoints_(nonWeightedBasisAtBasisEPoints), basisAtBasisEPoints_(basisAtBasisEPoints), hcurlBasisCurlAtBasisEPoints_(hcurlBasisCurlAtBasisEPoints), basisEWeights_(basisEWeights), wHCurlBasisAtBasisEPoints_(wHCurlBasisAtBasisEPoints), targetEWeights_(targetEWeights), @@ -215,7 +202,7 @@ struct ComputeHCurlBasisCoeffsOnCells_HDiv { ordinal_type idof = computedDofs_(i); for(ordinal_type iq=0; iq::getHDivBasisCoeffs(Kokkos::DynRankViewgetName(); ordinal_type numSides = cellBasis->getBaseCellTopology().getSideCount(); - ScalarViewType refSideNormal("refSideNormal", dim); ordinal_type numSideDofs(0); for(ordinal_type is=0; is::getHDivBasisCoeffs(Kokkos::DynRankViewgetAllEvalPoints(EvalPointsType::TARGET); auto targetDivEPoints = projStruct->getAllDerivEvalPoints(EvalPointsType::TARGET); - ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, dim); - ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints, dim); - { - ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",basisCardinality, numTotalBasisEPoints, dim); - ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",basisCardinality, numTotalTargetEPoints, dim); - cellBasis->getValues(nonOrientedBasisAtTargetEPoints, targetEPoints); - cellBasis->getValues(nonOrientedBasisAtBasisEPoints, basisEPoints); - - OrientationTools::modifyBasisByOrientation(basisAtBasisEPoints, nonOrientedBasisAtBasisEPoints, orts, cellBasis); - OrientationTools::modifyBasisByOrientation(basisAtTargetEPoints, nonOrientedBasisAtTargetEPoints, orts, cellBasis); - } + ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",basisCardinality, numTotalBasisEPoints, dim); + ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",basisCardinality, numTotalTargetEPoints, dim); + cellBasis->getValues(basisAtTargetEPoints, targetEPoints); + cellBasis->getValues(basisAtBasisEPoints, basisEPoints); ScalarViewType basisDivAtBasisDivEPoints; ScalarViewType basisDivAtTargetDivEPoints; if(numTotalTargetDivEPoints>0) { - ScalarViewType nonOrientedBasisDivAtTargetDivEPoints, nonOrientedBasisDivAtBasisDivEPoints; - basisDivAtBasisDivEPoints = ScalarViewType ("basisDivAtBasisDivEPoints",numCells,basisCardinality, numTotalBasisDivEPoints); - nonOrientedBasisDivAtBasisDivEPoints = ScalarViewType ("nonOrientedBasisDivAtBasisDivEPoints",basisCardinality, numTotalBasisDivEPoints); - basisDivAtTargetDivEPoints = ScalarViewType("basisDivAtTargetDivEPoints",numCells,basisCardinality, numTotalTargetDivEPoints); - nonOrientedBasisDivAtTargetDivEPoints = ScalarViewType("nonOrientedBasisDivAtTargetDivEPoints",basisCardinality, numTotalTargetDivEPoints); - cellBasis->getValues(nonOrientedBasisDivAtBasisDivEPoints, basisDivEPoints, OPERATOR_DIV); - cellBasis->getValues(nonOrientedBasisDivAtTargetDivEPoints, targetDivEPoints, OPERATOR_DIV); - - OrientationTools::modifyBasisByOrientation(basisDivAtBasisDivEPoints, nonOrientedBasisDivAtBasisDivEPoints, orts, cellBasis); - OrientationTools::modifyBasisByOrientation(basisDivAtTargetDivEPoints, nonOrientedBasisDivAtTargetDivEPoints, orts, cellBasis); + basisDivAtBasisDivEPoints = ScalarViewType ("basisDivAtBasisDivEPoints",basisCardinality, numTotalBasisDivEPoints); + basisDivAtTargetDivEPoints = ScalarViewType("basisDivAtTargetDivEPoints",basisCardinality, numTotalTargetDivEPoints); + cellBasis->getValues(basisDivAtBasisDivEPoints, basisDivEPoints, OPERATOR_DIV); + cellBasis->getValues(basisDivAtTargetDivEPoints, targetDivEPoints, OPERATOR_DIV); } - ScalarViewType refSidesNormal("refSidesNormal", numSides, dim); + ScalarViewType refBasisCoeffs("refBasisCoeffs", basisCoeffs.extent(0), basisCoeffs.extent(1)); + ordinal_type computedDofsCount = 0; for(ordinal_type is=0; is::getHDivBasisCoeffs(Kokkos::DynRankView::getReferenceSideNormal(sideNormal, is, cellTopo); + ScalarViewType refSideNormal("refSideNormal", dim); + CellTools::getReferenceSideNormal(refSideNormal, is, cellTopo); - ScalarViewType basisNormalAtBasisEPoints("normalBasisAtBasisEPoints",numCells,sideCardinality, numBasisEPoints); - ScalarViewType wBasisNormalAtBasisEPoints("wBasisNormalAtBasisEPoints",numCells,sideCardinality, numBasisEPoints); + ScalarViewType basisNormalAtBasisEPoints("normalBasisAtBasisEPoints",1,sideCardinality, numBasisEPoints); + ScalarViewType wBasisNormalAtBasisEPoints("wBasisNormalAtBasisEPoints",1,sideCardinality, numBasisEPoints); ScalarViewType wBasisNormalAtTargetEPoints("wBasisNormalAtTargetEPoints",numCells,sideCardinality, numTargetEPoints); ScalarViewType targetNormalAtTargetEPoints("targetNormalAtTargetEPoints",numCells, numTargetEPoints); @@ -344,17 +318,24 @@ ProjectionTools::getHDivBasisCoeffs(Kokkos::DynRankViewgetTargetEvalWeights(sideDim,is)); auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(sideDim,is)); + Kokkos::parallel_for(Kokkos::MDRangePolicy >({0,0}, {sideCardinality,numBasisEPoints}), + KOKKOS_LAMBDA (const ordinal_type &j, const ordinal_type &iq) { + ordinal_type jdof = tagToOrdinal(sideDim, is, j); + for(ordinal_type d=0; d ; - Kokkos::parallel_for(policy, functorTypeSide(basisNormalAtBasisEPoints, basisAtBasisEPoints, - basisEWeights, wBasisNormalAtBasisEPoints, targetEWeights, + + using functorTypeSide = FunctorsProjectionTools::ComputeBasisCoeffsOnSides_HDiv; + Kokkos::parallel_for(policy, functorTypeSide(targetEWeights, basisAtTargetEPoints, wBasisNormalAtTargetEPoints, tagToOrdinal, targetAtEPoints, targetNormalAtTargetEPoints, - refSidesNormal, sideCardinality, offsetBasis, + refSideNormal, sideCardinality, offsetTarget, sideDim, dim, is)); - ScalarViewType sideMassMat_("sideMassMat_", numCells, sideCardinality+1, sideCardinality+1), + ScalarViewType sideMassMat_("sideMassMat_", 1, sideCardinality+1, sideCardinality+1), sideRhsMat_("rhsMat_", numCells, sideCardinality+1); ScalarViewType targetEWeights_("targetEWeights", numCells, 1, targetEWeights.extent(0)); @@ -362,7 +343,7 @@ ProjectionTools::getHDivBasisCoeffs(Kokkos::DynRankView::integrate(Kokkos::subview(sideMassMat_, Kokkos::ALL(), range_H, range_H), basisNormalAtBasisEPoints, wBasisNormalAtBasisEPoints); @@ -377,115 +358,117 @@ ProjectionTools::getHDivBasisCoeffs(Kokkos::DynRankViewgetDofCount(dim,0); - if(numCellDofs==0) - return; - - Basis *hcurlBasis = NULL; - if(cellTopo.getKey() == shards::getCellTopologyData >()->key) - hcurlBasis = new Basis_HCURL_HEX_In_FEM(cellBasis->getDegree()); - else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) - hcurlBasis = new Basis_HCURL_TET_In_FEM(cellBasis->getDegree()); - else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) - hcurlBasis = new typename DerivedNodalBasisFamily::HCURL_WEDGE(cellBasis->getDegree()); - else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) - hcurlBasis = new Basis_HGRAD_QUAD_Cn_FEM(cellBasis->getDegree()); - else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) - hcurlBasis = new Basis_HGRAD_TRI_Cn_FEM(cellBasis->getDegree()); - else { - std::stringstream ss; - ss << ">>> ERROR (Intrepid2::ProjectionTools::getHDivBasisCoeffs): " - << "Method not implemented for basis " << name; - INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() ); - } - if(hcurlBasis == NULL) return; + if(numCellDofs!=0) { + Basis *hcurlBasis = NULL; + if(cellTopo.getKey() == shards::getCellTopologyData >()->key) + hcurlBasis = new Basis_HCURL_HEX_In_FEM(cellBasis->getDegree()); + else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) + hcurlBasis = new Basis_HCURL_TET_In_FEM(cellBasis->getDegree()); + else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) + hcurlBasis = new typename DerivedNodalBasisFamily::HCURL_WEDGE(cellBasis->getDegree()); + else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) + hcurlBasis = new Basis_HGRAD_QUAD_Cn_FEM(cellBasis->getDegree()); + else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) + hcurlBasis = new Basis_HGRAD_TRI_Cn_FEM(cellBasis->getDegree()); + else { + std::stringstream ss; + ss << ">>> ERROR (Intrepid2::ProjectionTools::getHDivBasisCoeffs): " + << "Method not implemented for basis " << name; + INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() ); + } + if(hcurlBasis == NULL) return; + + auto targetDivEPointsRange = projStruct->getTargetDerivPointsRange(); + auto basisDivEPointsRange = projStruct->getBasisDerivPointsRange(); - auto targetDivEPointsRange = projStruct->getTargetDerivPointsRange(); - auto basisDivEPointsRange = projStruct->getBasisDerivPointsRange(); + ordinal_type numTargetDivEPoints = range_size(targetDivEPointsRange(dim,0)); + ordinal_type numBasisDivEPoints = range_size(basisDivEPointsRange(dim,0)); - ordinal_type numTargetDivEPoints = range_size(targetDivEPointsRange(dim,0)); - ordinal_type numBasisDivEPoints = range_size(basisDivEPointsRange(dim,0)); + auto targetDivEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetDerivEvalWeights(dim,0)); + auto divEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisDerivEvalWeights(dim,0)); - auto targetDivEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetDerivEvalWeights(dim,0)); - auto divEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisDerivEvalWeights(dim,0)); + ordinal_type offsetBasisDiv = basisDivEPointsRange(dim,0).first; + ordinal_type offsetTargetDiv = targetDivEPointsRange(dim,0).first; - ordinal_type offsetBasisDiv = basisDivEPointsRange(dim,0).first; - ordinal_type offsetTargetDiv = targetDivEPointsRange(dim,0).first; + ScalarViewType weightedBasisDivAtBasisEPoints("weightedBasisDivAtBasisEPoints",numCells,numCellDofs, numBasisDivEPoints); + ScalarViewType weightedBasisDivAtTargetEPoints("weightedBasisDivAtTargetEPoints",numCells, numCellDofs, numTargetDivEPoints); + ScalarViewType basisDivAtBasisEPoints("basisDivAtBasisEPoints",1,numCellDofs, numBasisDivEPoints); + ScalarViewType targetSideDivAtBasisEPoints("targetSideDivAtBasisEPoints",numCells, numBasisDivEPoints); - ScalarViewType weightedBasisDivAtBasisEPoints("weightedBasisDivAtBasisEPoints",numCells,numCellDofs, numBasisDivEPoints); - ScalarViewType weightedBasisDivAtTargetEPoints("weightedBasisDivAtTargetEPoints",numCells, numCellDofs, numTargetDivEPoints); - ScalarViewType basisDivAtBasisEPoints("basisDivAtBasisEPoints",numCells,numCellDofs, numBasisDivEPoints); - ScalarViewType targetSideDivAtBasisEPoints("targetSideDivAtBasisEPoints",numCells, numBasisDivEPoints); + auto cellDofs = Kokkos::subview(tagToOrdinal, dim, 0, Kokkos::ALL()); + using functorType = FunctorsProjectionTools::ComputeBasisCoeffsOnCells_HDiv; + Kokkos::parallel_for(policy, functorType( refBasisCoeffs, targetSideDivAtBasisEPoints, basisDivAtBasisEPoints, + basisDivAtBasisDivEPoints, divEWeights, weightedBasisDivAtBasisEPoints, targetDivEWeights, basisDivAtTargetDivEPoints, weightedBasisDivAtTargetEPoints, + computedDofs, cellDofs, numCellDofs, offsetBasisDiv, offsetTargetDiv, numSideDofs)); - auto cellDofs = Kokkos::subview(tagToOrdinal, dim, 0, Kokkos::ALL()); - using functorType = FunctorsProjectionTools::ComputeBasisCoeffsOnCells_HDiv; - Kokkos::parallel_for(policy, functorType( basisCoeffs, targetSideDivAtBasisEPoints, basisDivAtBasisEPoints, - basisDivAtBasisDivEPoints, divEWeights, weightedBasisDivAtBasisEPoints, targetDivEWeights, basisDivAtTargetDivEPoints, weightedBasisDivAtTargetEPoints, - computedDofs, cellDofs, numCellDofs, offsetBasisDiv, offsetTargetDiv, numSideDofs)); + ordinal_type hcurlBasisCardinality = hcurlBasis->getCardinality(); + ordinal_type numCurlInteriorDOFs = hcurlBasis->getDofCount(dim,0); - ordinal_type hcurlBasisCardinality = hcurlBasis->getCardinality(); - ordinal_type numCurlInteriorDOFs = hcurlBasis->getDofCount(dim,0); + range_type range_H(0, numCellDofs); + range_type range_B(numCellDofs, numCellDofs+numCurlInteriorDOFs); - range_type range_H(0, numCellDofs); - range_type range_B(numCellDofs, numCellDofs+numCurlInteriorDOFs); + ScalarViewType massMat_("massMat_",1,numCellDofs+numCurlInteriorDOFs,numCellDofs+numCurlInteriorDOFs), + rhsMatTrans("rhsMatTrans",numCells,numCellDofs+numCurlInteriorDOFs); - ScalarViewType massMat_("massMat_",numCells,numCellDofs+numCurlInteriorDOFs,numCellDofs+numCurlInteriorDOFs), - rhsMatTrans("rhsMatTrans",numCells,numCellDofs+numCurlInteriorDOFs); + FunctionSpaceTools::integrate(Kokkos::subview(massMat_, Kokkos::ALL(), range_H,range_H), basisDivAtBasisEPoints, Kokkos::subview(weightedBasisDivAtBasisEPoints, std::make_pair(0,1), Kokkos::ALL(), Kokkos::ALL()) ); + FunctionSpaceTools::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_H), targetDivAtDivEPoints, weightedBasisDivAtTargetEPoints); + FunctionSpaceTools::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_H), targetSideDivAtBasisEPoints, weightedBasisDivAtBasisEPoints,true); - FunctionSpaceTools::integrate(Kokkos::subview(massMat_, Kokkos::ALL(), range_H,range_H), basisDivAtBasisEPoints, weightedBasisDivAtBasisEPoints); - FunctionSpaceTools::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_H), targetDivAtDivEPoints, weightedBasisDivAtTargetEPoints); - FunctionSpaceTools::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_H), targetSideDivAtBasisEPoints, weightedBasisDivAtBasisEPoints,true); + if(numCurlInteriorDOFs>0){ + ordinal_type numTargetEPoints = range_size(targetEPointsRange(dim,0)); + ordinal_type numBasisEPoints = range_size(basisEPointsRange(dim,0)); - if(numCurlInteriorDOFs>0){ - ordinal_type numTargetEPoints = range_size(targetEPointsRange(dim,0)); - ordinal_type numBasisEPoints = range_size(basisEPointsRange(dim,0)); + auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(dim,0)); + auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(dim,0)); - auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(dim,0)); - auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(dim,0)); + ordinal_type offsetBasis = basisEPointsRange(dim,0).first; - ordinal_type offsetBasis = basisEPointsRange(dim,0).first; + ScalarViewType negPartialProjAtBasisEPoints("targetSideAtBasisEPoints",numCells, numBasisEPoints, dim); + ScalarViewType nonWeightedBasisAtBasisEPoints("basisAtBasisEPoints",1,numCellDofs, numBasisEPoints, dim); + ScalarViewType hcurlBasisCurlAtBasisEPoints("hcurlBasisCurlAtBasisEPoints",hcurlBasisCardinality, numBasisEPoints,dim); + ScalarViewType hcurlBasisCurlAtTargetEPoints("hcurlBasisCurlAtTargetEPoints", hcurlBasisCardinality,numTargetEPoints, dim); + ScalarViewType wHcurlBasisCurlAtBasisEPoints("wHcurlBasisHcurlAtBasisEPoints", numCells, numCurlInteriorDOFs, numBasisEPoints,dim); + ScalarViewType wHcurlBasisCurlAtTargetEPoints("wHcurlBasisHcurlAtTargetEPoints",numCells, numCurlInteriorDOFs, numTargetEPoints,dim); - ScalarViewType negPartialProjAtBasisEPoints("targetSideAtBasisEPoints",numCells, numBasisEPoints, dim); - ScalarViewType nonWeightedBasisAtBasisEPoints("basisAtBasisEPoints",numCells,numCellDofs, numBasisEPoints, dim); - ScalarViewType hcurlBasisCurlAtBasisEPoints("hcurlBasisCurlAtBasisEPoints",hcurlBasisCardinality, numBasisEPoints,dim); - ScalarViewType hcurlBasisCurlAtTargetEPoints("hcurlBasisCurlAtTargetEPoints", hcurlBasisCardinality,numTargetEPoints, dim); - ScalarViewType wHcurlBasisCurlAtBasisEPoints("wHcurlBasisHcurlAtBasisEPoints", numCells, numCurlInteriorDOFs, numBasisEPoints,dim); - ScalarViewType wHcurlBasisCurlAtTargetEPoints("wHcurlBasisHcurlAtTargetEPoints",numCells, numCurlInteriorDOFs, numTargetEPoints,dim); + hcurlBasis->getValues(hcurlBasisCurlAtBasisEPoints, Kokkos::subview(basisEPoints, basisEPointsRange(dim,0), Kokkos::ALL()), OPERATOR_CURL); + hcurlBasis->getValues(hcurlBasisCurlAtTargetEPoints, Kokkos::subview(targetEPoints,targetEPointsRange(dim,0),Kokkos::ALL()), OPERATOR_CURL); - hcurlBasis->getValues(hcurlBasisCurlAtBasisEPoints, Kokkos::subview(basisEPoints, basisEPointsRange(dim,0), Kokkos::ALL()), OPERATOR_CURL); - hcurlBasis->getValues(hcurlBasisCurlAtTargetEPoints, Kokkos::subview(targetEPoints,targetEPointsRange(dim,0),Kokkos::ALL()), OPERATOR_CURL); + auto hCurlTagToOrdinal = Kokkos::create_mirror_view_and_copy(MemSpaceType(), hcurlBasis->getAllDofOrdinal()); + auto cellHCurlDof = Kokkos::subview(hCurlTagToOrdinal, dim, 0, range_type(0, numCurlInteriorDOFs)); - auto hCurlTagToOrdinal = Kokkos::create_mirror_view_and_copy(MemSpaceType(), hcurlBasis->getAllDofOrdinal()); - auto cellHCurlDof = Kokkos::subview(hCurlTagToOrdinal, dim, 0, range_type(0, numCurlInteriorDOFs)); + using functorTypeHCurlCells = FunctorsProjectionTools::ComputeHCurlBasisCoeffsOnCells_HDiv; + Kokkos::parallel_for(policy, functorTypeHCurlCells(refBasisCoeffs, negPartialProjAtBasisEPoints, nonWeightedBasisAtBasisEPoints, + basisAtBasisEPoints, hcurlBasisCurlAtBasisEPoints, basisEWeights, wHcurlBasisCurlAtBasisEPoints, targetEWeights, + hcurlBasisCurlAtTargetEPoints, wHcurlBasisCurlAtTargetEPoints, tagToOrdinal, computedDofs, cellHCurlDof, + numCellDofs, offsetBasis, numSideDofs, dim)); - using functorTypeHCurlCells = FunctorsProjectionTools::ComputeHCurlBasisCoeffsOnCells_HDiv; - Kokkos::parallel_for(policy, functorTypeHCurlCells(basisCoeffs, negPartialProjAtBasisEPoints, nonWeightedBasisAtBasisEPoints, - basisAtBasisEPoints, hcurlBasisCurlAtBasisEPoints, basisEWeights, wHcurlBasisCurlAtBasisEPoints, targetEWeights, - hcurlBasisCurlAtTargetEPoints, wHcurlBasisCurlAtTargetEPoints, tagToOrdinal, computedDofs, cellHCurlDof, - numCellDofs, offsetBasis, numSideDofs, dim)); + FunctionSpaceTools::integrate(Kokkos::subview(massMat_, Kokkos::ALL(), range_H,range_B), nonWeightedBasisAtBasisEPoints, Kokkos::subview(wHcurlBasisCurlAtBasisEPoints, std::make_pair(0,1), Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()) ); + FunctionSpaceTools::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_B), Kokkos::subview(targetAtEPoints, Kokkos::ALL(), targetEPointsRange(dim,0), Kokkos::ALL()), wHcurlBasisCurlAtTargetEPoints); + FunctionSpaceTools::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_B), negPartialProjAtBasisEPoints, wHcurlBasisCurlAtBasisEPoints,true); + } + delete hcurlBasis; + + typedef Kokkos::DynRankView WorkArrayViewType; + ScalarViewType t_("t",numCells, numCellDofs+numCurlInteriorDOFs); + WorkArrayViewType w_("w",numCells, numCellDofs+numCurlInteriorDOFs); - FunctionSpaceTools::integrate(Kokkos::subview(massMat_, Kokkos::ALL(), range_H,range_B), nonWeightedBasisAtBasisEPoints, wHcurlBasisCurlAtBasisEPoints); - FunctionSpaceTools::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_B), Kokkos::subview(targetAtEPoints, Kokkos::ALL(), targetEPointsRange(dim,0), Kokkos::ALL()), wHcurlBasisCurlAtTargetEPoints); - FunctionSpaceTools::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_B), negPartialProjAtBasisEPoints, wHcurlBasisCurlAtBasisEPoints,true); + ElemSystem cellSystem("cellSystem", true); + cellSystem.solve(refBasisCoeffs, massMat_, rhsMatTrans, t_, w_, cellDofs, numCellDofs, numCurlInteriorDOFs); } - delete hcurlBasis; - typedef Kokkos::DynRankView WorkArrayViewType; - ScalarViewType t_("t",numCells, numCellDofs+numCurlInteriorDOFs); - WorkArrayViewType w_("w",numCells, numCellDofs+numCurlInteriorDOFs); + OrientationTools::modifyBasisByOrientationInverse(basisCoeffs, refBasisCoeffs, orts, cellBasis, true); - ElemSystem cellSystem("cellSystem", true); - cellSystem.solve(basisCoeffs, massMat_, rhsMatTrans, t_, w_, cellDofs, numCellDofs, numCurlInteriorDOFs); } } // Intrepid2 namespace diff --git a/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHGRAD.hpp b/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHGRAD.hpp index 992988935b43..215a49c3da35 100644 --- a/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHGRAD.hpp +++ b/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHGRAD.hpp @@ -84,28 +84,28 @@ struct ComputeBasisCoeffsOnVertices_HGRAD { ordinal_type idof = tagToOrdinal_(0, iv, 0); ordinal_type pt = targetEPointsRange_(0,iv).first; //the value of the basis at the vertex might be different than 1; HGrad basis, so the function is scalar - basisCoeffs_(ic,idof) = targetAtTargetEPoints_(ic,pt)/basisAtTargetEPoints_(ic,idof,pt,0); + basisCoeffs_(ic,idof) = targetAtTargetEPoints_(ic,pt)/basisAtTargetEPoints_(idof,pt,0); } } }; template +typename ViewType4, typename ViewType5> struct ComputeBasisCoeffsOnEdges_HGRAD { const ViewType1 basisCoeffs_; - const ViewType2 negPartialProjGrad_; - const ViewType2 basisTanAtEPoints_; - const ViewType2 basisGradAtBasisGradEPoints_; - const ViewType3 basisGradEWeights_; - const ViewType2 wBasisAtBasisGradEPoints_; - const ViewType3 targetGradEWeights_; - const ViewType2 basisGradAtTargetGradEPoints_; - const ViewType2 wBasisAtTargetGradEPoints_; - const ViewType4 computedDofs_; - const ViewType5 tagToOrdinal_; - const ViewType2 targetGradTanAtTargetGradEPoints_; - const ViewType6 targetGradAtTargetGradEPoints_; - const ViewType2 refEdgesTan_; + const ViewType1 negPartialProjGrad_; + const ViewType1 basisTanAtEPoints_; + const ViewType1 basisGradAtBasisGradEPoints_; + const ViewType2 basisGradEWeights_; + const ViewType1 wBasisAtBasisGradEPoints_; + const ViewType2 targetGradEWeights_; + const ViewType1 basisGradAtTargetGradEPoints_; + const ViewType1 wBasisAtTargetGradEPoints_; + const ViewType3 computedDofs_; + const ViewType4 tagToOrdinal_; + const ViewType1 targetGradTanAtTargetGradEPoints_; + const ViewType5 targetGradAtTargetGradEPoints_; + const ViewType1 refEdgeTan_; ordinal_type edgeCardinality_; ordinal_type offsetBasis_; ordinal_type offsetTarget_; @@ -114,17 +114,17 @@ struct ComputeBasisCoeffsOnEdges_HGRAD { ordinal_type dim_; ordinal_type iedge_; - ComputeBasisCoeffsOnEdges_HGRAD(const ViewType1 basisCoeffs, ViewType2 negPartialProjGrad, const ViewType2 basisTanAtEPoints, - const ViewType2 basisGradAtBasisGradEPoints, const ViewType3 basisGradEWeights, const ViewType2 wBasisAtBasisGradEPoints, const ViewType3 targetGradEWeights, - const ViewType2 basisGradAtTargetGradEPoints, const ViewType2 wBasisAtTargetGradEPoints, const ViewType4 computedDofs, const ViewType5 tagToOrdinal, - const ViewType2 targetGradTanAtTargetGradEPoints, const ViewType6 targetGradAtTargetGradEPoints, const ViewType2 refEdgesTan, + ComputeBasisCoeffsOnEdges_HGRAD(const ViewType1 basisCoeffs, ViewType1 negPartialProjGrad, const ViewType1 basisTanAtEPoints, + const ViewType1 basisGradAtBasisGradEPoints, const ViewType2 basisGradEWeights, const ViewType1 wBasisAtBasisGradEPoints, const ViewType2 targetGradEWeights, + const ViewType1 basisGradAtTargetGradEPoints, const ViewType1 wBasisAtTargetGradEPoints, const ViewType3 computedDofs, const ViewType4 tagToOrdinal, + const ViewType1 targetGradTanAtTargetGradEPoints, const ViewType5 targetGradAtTargetGradEPoints, const ViewType1 refEdgeTan, ordinal_type edgeCardinality, ordinal_type offsetBasis, ordinal_type offsetTarget, ordinal_type numVertexDofs, ordinal_type edgeDim, ordinal_type dim, ordinal_type iedge) : basisCoeffs_(basisCoeffs), negPartialProjGrad_(negPartialProjGrad), basisTanAtEPoints_(basisTanAtEPoints), basisGradAtBasisGradEPoints_(basisGradAtBasisGradEPoints), basisGradEWeights_(basisGradEWeights), wBasisAtBasisGradEPoints_(wBasisAtBasisGradEPoints), targetGradEWeights_(targetGradEWeights), basisGradAtTargetGradEPoints_(basisGradAtTargetGradEPoints), wBasisAtTargetGradEPoints_(wBasisAtTargetGradEPoints), computedDofs_(computedDofs), tagToOrdinal_(tagToOrdinal), targetGradTanAtTargetGradEPoints_(targetGradTanAtTargetGradEPoints), - targetGradAtTargetGradEPoints_(targetGradAtTargetGradEPoints), refEdgesTan_(refEdgesTan), + targetGradAtTargetGradEPoints_(targetGradAtTargetGradEPoints), refEdgeTan_(refEdgeTan), edgeCardinality_(edgeCardinality), offsetBasis_(offsetBasis), offsetTarget_(offsetTarget), numVertexDofs_(numVertexDofs), edgeDim_(edgeDim), dim_(dim), iedge_(iedge) {} @@ -135,53 +135,54 @@ struct ComputeBasisCoeffsOnEdges_HGRAD { ordinal_type numBasisGradEPoints = basisGradEWeights_.extent(0); ordinal_type numTargetGradEPoints = targetGradEWeights_.extent(0); + + //this loop could be moved outside of cell loop for(ordinal_type j=0; j +typename ViewType4, typename ViewType5> struct ComputeBasisCoeffsOnFaces_HGRAD { const ViewType1 basisCoeffs_; - const ViewType2 negPartialProjGrad_; - const ViewType2 faceBasisGradAtGradEPoints_; - const ViewType2 basisGradAtBasisGradEPoints_; - const ViewType3 basisGradEWeights_; - const ViewType2 wBasisGradAtGradEPoints_; - const ViewType3 targetGradEWeights_; - const ViewType2 basisGradAtTargetGradEPoints_; - const ViewType2 wBasisGradBasisAtTargetGradEPoints_; - const ViewType4 computedDofs_; - const ViewType5 tagToOrdinal_; - const ViewType6 orts_; - const ViewType2 targetGradTanAtTargetGradEPoints_; - const ViewType7 targetGradAtTargetGradEPoints_; - const ViewType2 faceCoeff_; - const ViewType8 faceParametrization_; + const ViewType1 negPartialProjGrad_; + const ViewType1 faceBasisGradAtGradEPoints_; + const ViewType1 basisGradAtBasisGradEPoints_; + const ViewType2 basisGradEWeights_; + const ViewType1 wBasisGradAtGradEPoints_; + const ViewType2 targetGradEWeights_; + const ViewType1 basisGradAtTargetGradEPoints_; + const ViewType1 wBasisGradBasisAtTargetGradEPoints_; + const ViewType3 computedDofs_; + const ViewType4 tagToOrdinal_; + const ViewType1 targetGradTanAtTargetGradEPoints_; + const ViewType5 targetGradAtTargetGradEPoints_; + const ViewType1 refFaceNormal_; ordinal_type faceCardinality_; ordinal_type offsetBasisGrad_; ordinal_type offsetTargetGrad_; @@ -190,42 +191,34 @@ struct ComputeBasisCoeffsOnFaces_HGRAD { ordinal_type faceDim_; ordinal_type dim_; ordinal_type iface_; - unsigned topoKey_; - ComputeBasisCoeffsOnFaces_HGRAD(const ViewType1 basisCoeffs, ViewType2 negPartialProjGrad, const ViewType2 faceBasisGradAtGradEPoints, - const ViewType2 basisGradAtBasisGradEPoints, const ViewType3 basisGradEWeights, const ViewType2 wBasisGradAtGradEPoints, const ViewType3 targetGradEWeights, - const ViewType2 basisGradAtTargetGradEPoints, const ViewType2 wBasisGradBasisAtTargetGradEPoints, const ViewType4 computedDofs, const ViewType5 tagToOrdinal, - const ViewType6 orts, const ViewType2 targetGradTanAtTargetGradEPoints, const ViewType7 targetGradAtTargetGradEPoints, - const ViewType8 faceParametrization, ordinal_type faceCardinality, ordinal_type offsetBasisGrad, + ComputeBasisCoeffsOnFaces_HGRAD(const ViewType1 basisCoeffs, ViewType1 negPartialProjGrad, const ViewType1 faceBasisGradAtGradEPoints, + const ViewType1 basisGradAtBasisGradEPoints, const ViewType2 basisGradEWeights, const ViewType1 wBasisGradAtGradEPoints, const ViewType2 targetGradEWeights, + const ViewType1 basisGradAtTargetGradEPoints, const ViewType1 wBasisGradBasisAtTargetGradEPoints, const ViewType3 computedDofs, const ViewType4 tagToOrdinal, + const ViewType1 targetGradTanAtTargetGradEPoints, const ViewType5 targetGradAtTargetGradEPoints, + const ViewType1 refFaceNormal, ordinal_type faceCardinality, ordinal_type offsetBasisGrad, ordinal_type offsetTargetGrad, ordinal_type numVertexEdgeDofs, ordinal_type numFaces, ordinal_type faceDim, - ordinal_type dim, ordinal_type iface, unsigned topoKey) : + ordinal_type dim, ordinal_type iface) : basisCoeffs_(basisCoeffs), negPartialProjGrad_(negPartialProjGrad), faceBasisGradAtGradEPoints_(faceBasisGradAtGradEPoints), basisGradAtBasisGradEPoints_(basisGradAtBasisGradEPoints), basisGradEWeights_(basisGradEWeights), wBasisGradAtGradEPoints_(wBasisGradAtGradEPoints), targetGradEWeights_(targetGradEWeights), basisGradAtTargetGradEPoints_(basisGradAtTargetGradEPoints), wBasisGradBasisAtTargetGradEPoints_(wBasisGradBasisAtTargetGradEPoints), - computedDofs_(computedDofs), tagToOrdinal_(tagToOrdinal), orts_(orts), targetGradTanAtTargetGradEPoints_(targetGradTanAtTargetGradEPoints), - targetGradAtTargetGradEPoints_(targetGradAtTargetGradEPoints), faceParametrization_(faceParametrization), + computedDofs_(computedDofs), tagToOrdinal_(tagToOrdinal), targetGradTanAtTargetGradEPoints_(targetGradTanAtTargetGradEPoints), + targetGradAtTargetGradEPoints_(targetGradAtTargetGradEPoints), refFaceNormal_(refFaceNormal), faceCardinality_(faceCardinality), offsetBasisGrad_(offsetBasisGrad), offsetTargetGrad_(offsetTargetGrad), numVertexEdgeDofs_(numVertexEdgeDofs), numFaces_(numFaces), - faceDim_(faceDim), dim_(dim), iface_(iface), topoKey_(topoKey) + faceDim_(faceDim), dim_(dim), iface_(iface) {} void KOKKOS_INLINE_FUNCTION operator()(const ordinal_type ic) const { - typename ViewType2::value_type data[3*3]; - auto tangentsAndNormal = ViewType2(data, 3, 3); - - ordinal_type fOrt[6]; - orts_(ic).getFaceOrientation(fOrt, numFaces_); - ordinal_type ort = fOrt[iface_]; - ordinal_type numBasisGradEPoints = basisGradEWeights_.extent(0); ordinal_type numTargetGradEPoints = targetGradEWeights_.extent(0); - Impl::OrientationTools::getRefSideTangentsAndNormal(tangentsAndNormal, faceParametrization_,topoKey_, iface_, ort); + //normal - typename ViewType2::value_type n[3] = {tangentsAndNormal(2,0), tangentsAndNormal(2,1), tangentsAndNormal(2,2)}; + typename ViewType2::value_type n[3] = {refFaceNormal_(0), refFaceNormal_(1), refFaceNormal_(2)}; for(ordinal_type j=0; j +typename ViewType4> struct ComputeBasisCoeffsOnCells_HGRAD { const ViewType1 basisCoeffs_; - const ViewType2 negPartialProjGrad_; - const ViewType2 cellBasisGradAtGradEPoints_; - const ViewType2 basisGradAtBasisGradEPoints_; - const ViewType3 basisGradEWeights_; - const ViewType2 wBasisGradAtGradEPoints_; - const ViewType3 targetGradEWeights_; - const ViewType2 basisGradAtTargetGradEPoints_; - const ViewType2 wBasisGradBasisAtTargetGradEPoints_; - const ViewType4 computedDofs_; - const ViewType5 elemDof_; + const ViewType1 negPartialProjGrad_; + const ViewType1 cellBasisGradAtGradEPoints_; + const ViewType1 basisGradAtBasisGradEPoints_; + const ViewType2 basisGradEWeights_; + const ViewType1 wBasisGradAtGradEPoints_; + const ViewType2 targetGradEWeights_; + const ViewType1 basisGradAtTargetGradEPoints_; + const ViewType1 wBasisGradBasisAtTargetGradEPoints_; + const ViewType3 computedDofs_; + const ViewType4 elemDof_; ordinal_type dim_; ordinal_type numElemDofs_; ordinal_type offsetBasisGrad_; ordinal_type offsetTargetGrad_; ordinal_type numVertexEdgeFaceDofs_; - ComputeBasisCoeffsOnCells_HGRAD(const ViewType1 basisCoeffs, ViewType2 negPartialProjGrad, const ViewType2 cellBasisGradAtGradEPoints, - const ViewType2 basisGradAtBasisGradEPoints, const ViewType3 basisGradEWeights, const ViewType2 wBasisGradAtGradEPoints, const ViewType3 targetGradEWeights, - const ViewType2 basisGradAtTargetGradEPoints, const ViewType2 wBasisGradBasisAtTargetGradEPoints, const ViewType4 computedDofs, const ViewType5 elemDof, + ComputeBasisCoeffsOnCells_HGRAD(const ViewType1 basisCoeffs, ViewType1 negPartialProjGrad, const ViewType1 cellBasisGradAtGradEPoints, + const ViewType1 basisGradAtBasisGradEPoints, const ViewType2 basisGradEWeights, const ViewType1 wBasisGradAtGradEPoints, const ViewType2 targetGradEWeights, + const ViewType1 basisGradAtTargetGradEPoints, const ViewType1 wBasisGradBasisAtTargetGradEPoints, const ViewType3 computedDofs, const ViewType4 elemDof, ordinal_type dim, ordinal_type numElemDofs, ordinal_type offsetBasisGrad, ordinal_type offsetTargetGrad, ordinal_type numVertexEdgeFaceDofs) : basisCoeffs_(basisCoeffs), negPartialProjGrad_(negPartialProjGrad), cellBasisGradAtGradEPoints_(cellBasisGradAtGradEPoints), basisGradAtBasisGradEPoints_(basisGradAtBasisGradEPoints), basisGradEWeights_(basisGradEWeights), wBasisGradAtGradEPoints_(wBasisGradAtGradEPoints), targetGradEWeights_(targetGradEWeights), @@ -301,11 +294,11 @@ struct ComputeBasisCoeffsOnCells_HGRAD { ordinal_type idof = elemDof_(j); for(ordinal_type d=0; d ::getHGradBasisCoeffs(BasisCoeffsViewType basisCoeffs ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0; ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount() : 0; - ScalarViewType refEdgesTan("refEdgesTan", numEdges, dim); + ScalarViewType refEdgeTan("refTan", dim); + ScalarViewType refFaceNormal("refNormal", dim); ordinal_type numVertexDofs = numVertices; @@ -372,35 +366,23 @@ ProjectionTools::getHGradBasisCoeffs(BasisCoeffsViewType basisCoeffs auto targetGradEPoints = projStruct->getAllDerivEvalPoints(EvalPointsType::TARGET); - ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints); - { - ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",basisCardinality, numTotalTargetEPoints); - cellBasis->getValues(nonOrientedBasisAtTargetEPoints, targetEPoints); - OrientationTools::modifyBasisByOrientation(basisAtTargetEPoints, nonOrientedBasisAtTargetEPoints, orts, cellBasis); - } + ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",basisCardinality, numTotalTargetEPoints); + cellBasis->getValues(basisAtTargetEPoints, targetEPoints); ScalarViewType basisGradAtBasisGradEPoints; ScalarViewType basisGradAtTargetGradEPoints; if(numTotalBasisGradEPoints>0) { - ScalarViewType nonOrientedBasisGradAtTargetGradEPoints, nonOrientedBasisGradAtBasisGradEPoints; - - basisGradAtBasisGradEPoints = ScalarViewType ("basisGradAtBasisGradEPoints",numCells,basisCardinality, numTotalBasisGradEPoints, dim); - nonOrientedBasisGradAtBasisGradEPoints = ScalarViewType ("nonOrientedBasisGradAtBasisGradEPoints",basisCardinality, numTotalBasisGradEPoints, dim); - basisGradAtTargetGradEPoints = ScalarViewType("basisGradAtTargetGradEPoints",numCells,basisCardinality, numTotalTargetGradEPoints, dim); - nonOrientedBasisGradAtTargetGradEPoints = ScalarViewType("nonOrientedBasisGradAtTargetGradEPoints",basisCardinality, numTotalTargetGradEPoints, dim); + basisGradAtBasisGradEPoints = ScalarViewType ("basisGradAtBasisGradEPoints",basisCardinality, numTotalBasisGradEPoints, dim); + basisGradAtTargetGradEPoints = ScalarViewType("basisGradAtTargetGradEPoints",basisCardinality, numTotalTargetGradEPoints, dim); - cellBasis->getValues(nonOrientedBasisGradAtBasisGradEPoints, basisGradEPoints, OPERATOR_GRAD); - cellBasis->getValues(nonOrientedBasisGradAtTargetGradEPoints, targetGradEPoints, OPERATOR_GRAD); - OrientationTools::modifyBasisByOrientation(basisGradAtBasisGradEPoints, nonOrientedBasisGradAtBasisGradEPoints, orts, cellBasis); - OrientationTools::modifyBasisByOrientation(basisGradAtTargetGradEPoints, nonOrientedBasisGradAtTargetGradEPoints, orts, cellBasis); + cellBasis->getValues(basisGradAtBasisGradEPoints, basisGradEPoints, OPERATOR_GRAD); + cellBasis->getValues(basisGradAtTargetGradEPoints, targetGradEPoints, OPERATOR_GRAD); } auto targetEPointsRange = Kokkos::create_mirror_view_and_copy(MemSpaceType(), projStruct->getTargetPointsRange()); auto targetGradEPointsRange = projStruct->getTargetDerivPointsRange(); auto basisGradEPointsRange = projStruct->getBasisDerivPointsRange(); - auto refTopologyKey = projStruct->getTopologyKey(); - auto tagToOrdinal = Kokkos::create_mirror_view_and_copy(MemSpaceType(), cellBasis->getAllDofOrdinal()); typename RefSubcellParametrization::ConstViewType subcellParamFace; @@ -414,13 +396,13 @@ ProjectionTools::getHGradBasisCoeffs(BasisCoeffsViewType basisCoeffs }); ordinal_type computedDofsCount = numVertices; + ScalarViewType refBasisCoeffs("refBasisCoeffs", basisCoeffs.extent(0), basisCoeffs.extent(1)); + const Kokkos::RangePolicy policy(0, numCells); - using functorType = FunctorsProjectionTools::ComputeBasisCoeffsOnVertices_HGRAD; - Kokkos::parallel_for(policy, functorType(basisCoeffs, tagToOrdinal, targetEPointsRange, + Kokkos::parallel_for(policy, functorType(refBasisCoeffs, tagToOrdinal, targetEPointsRange, targetAtTargetEPoints, basisAtTargetEPoints, numVertices)); - -// printView(basisCoeffs, std::cout, "after vertex evaluation, basisCoeffs"); for(ordinal_type ie=0; ie::getHGradBasisCoeffs(BasisCoeffsViewType basisCoeffs auto basisGradEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisDerivEvalWeights(edgeDim,ie)); auto targetGradEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetDerivEvalWeights(edgeDim,ie)); - auto edgeTan = Kokkos::subview(refEdgesTan, ie, Kokkos::ALL()); - CellTools::getReferenceEdgeTangent(edgeTan, ie, cellTopo); + CellTools::getReferenceEdgeTangent(refEdgeTan, ie, cellTopo); - ScalarViewType basisTanAtEPoints("basisTanAtEPoints",numCells,edgeCardinality, numBasisGradEPoints); + ScalarViewType basisTanAtEPoints("basisTanAtEPoints",1,edgeCardinality, numBasisGradEPoints); ScalarViewType targetGradTanAtTargetGradEPoints("tanBasisAtTargetGradEPoints",numCells, numTargetGradEPoints); ScalarViewType wBasisAtBasisGradEPoints("wTanBasisAtBasisGradEPoints",numCells,edgeCardinality, numBasisGradEPoints); ScalarViewType wBasisAtTargetGradEPoints("wTanBasisAtTargetGradEPoints",numCells,edgeCardinality, numTargetGradEPoints); ScalarViewType negPartialProjGrad("negPartialProjGrad", numCells, numBasisGradEPoints); - using functorTypeEdge = FunctorsProjectionTools::ComputeBasisCoeffsOnEdges_HGRAD; - Kokkos::parallel_for(policy, functorTypeEdge(basisCoeffs, negPartialProjGrad, basisTanAtEPoints, + Kokkos::parallel_for(policy, functorTypeEdge(refBasisCoeffs, negPartialProjGrad, basisTanAtEPoints, basisGradAtBasisGradEPoints, basisGradEWeights, wBasisAtBasisGradEPoints, targetGradEWeights, basisGradAtTargetGradEPoints, wBasisAtTargetGradEPoints, computedDofs, tagToOrdinal, - targetGradTanAtTargetGradEPoints, targetGradAtTargetGradEPoints, refEdgesTan, + targetGradTanAtTargetGradEPoints, targetGradAtTargetGradEPoints, refEdgeTan, edgeCardinality, offsetBasis, offsetTarget, numVertexDofs, edgeDim, dim, ie)); - ScalarViewType edgeMassMat_("edgeMassMat_", numCells, edgeCardinality, edgeCardinality), + ScalarViewType edgeMassMat_("edgeMassMat_", 1, edgeCardinality, edgeCardinality), edgeRhsMat_("rhsMat_", numCells, edgeCardinality); - FunctionSpaceTools::integrate(edgeMassMat_, basisTanAtEPoints, wBasisAtBasisGradEPoints); + FunctionSpaceTools::integrate(edgeMassMat_, basisTanAtEPoints, Kokkos::subview(wBasisAtBasisGradEPoints, std::make_pair(0,1), Kokkos::ALL(), Kokkos::ALL()) ); FunctionSpaceTools::integrate(edgeRhsMat_, targetGradTanAtTargetGradEPoints, wBasisAtTargetGradEPoints); FunctionSpaceTools::integrate(edgeRhsMat_, negPartialProjGrad, wBasisAtBasisGradEPoints,true); @@ -464,8 +445,8 @@ ProjectionTools::getHGradBasisCoeffs(BasisCoeffsViewType basisCoeffs auto edgeDofs = Kokkos::subview(tagToOrdinal, edgeDim, ie, Kokkos::ALL()); - ElemSystem edgeSystem("edgeSystem", false); - edgeSystem.solve(basisCoeffs, edgeMassMat_, edgeRhsMat_, t_, w_, edgeDofs, edgeCardinality); + ElemSystem edgeSystem("edgeSystem", true); + edgeSystem.solve(refBasisCoeffs, edgeMassMat_, edgeRhsMat_, t_, w_, edgeDofs, edgeCardinality); Kokkos::parallel_for("Compute Dofs ", Kokkos::RangePolicy (0, edgeCardinality), KOKKOS_LAMBDA (const size_t i) { @@ -477,8 +458,6 @@ ProjectionTools::getHGradBasisCoeffs(BasisCoeffsViewType basisCoeffs for(ordinal_type iface=0; ifacegetDofCount(faceDim,iface); ordinal_type numGradEPoints = range_size(basisGradEPointsRange(faceDim, iface)); @@ -488,28 +467,29 @@ ProjectionTools::getHGradBasisCoeffs(BasisCoeffsViewType basisCoeffs auto basisGradEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisDerivEvalWeights(faceDim,iface)); auto targetGradEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetDerivEvalWeights(faceDim,iface)); - ScalarViewType faceBasisGradAtGradEPoints("faceBasisGradAtGradEPoints",numCells,faceCardinality, numGradEPoints,dim); + CellTools::getReferenceFaceNormal(refFaceNormal, iface, cellTopo); + + ScalarViewType faceBasisGradAtGradEPoints("faceBasisGradAtGradEPoints",1,faceCardinality, numGradEPoints,dim); ScalarViewType wBasisGradAtGradEPoints("wNormalBasisGradAtGradEPoints",numCells,faceCardinality, numGradEPoints,dim); ScalarViewType wBasisGradBasisAtTargetGradEPoints("wNormalBasisGradAtTargetGradEPoints",numCells,faceCardinality, numTargetGradEPoints,dim); ScalarViewType targetGradTanAtTargetGradEPoints("targetGradTanAtTargetGradEPoints",numCells, numTargetGradEPoints,dim); ScalarViewType negPartialProjGrad("mNormalComputedProjection", numCells,numGradEPoints,dim); - using functorTypeFace_HGRAD = FunctorsProjectionTools::ComputeBasisCoeffsOnFaces_HGRAD; + using functorTypeFace_HGRAD = FunctorsProjectionTools::ComputeBasisCoeffsOnFaces_HGRAD; - Kokkos::parallel_for(policy, functorTypeFace_HGRAD(basisCoeffs, negPartialProjGrad, faceBasisGradAtGradEPoints, + Kokkos::parallel_for(policy, functorTypeFace_HGRAD(refBasisCoeffs, negPartialProjGrad, faceBasisGradAtGradEPoints, basisGradAtBasisGradEPoints, basisGradEWeights, wBasisGradAtGradEPoints, targetGradEWeights, basisGradAtTargetGradEPoints,wBasisGradBasisAtTargetGradEPoints, computedDofs, tagToOrdinal, - orts,targetGradTanAtTargetGradEPoints,targetGradAtTargetGradEPoints, - subcellParamFace, faceCardinality, offsetBasisGrad, + targetGradTanAtTargetGradEPoints,targetGradAtTargetGradEPoints, + refFaceNormal, faceCardinality, offsetBasisGrad, offsetTargetGrad, numVertexDofs+numEdgeDofs, numFaces, faceDim, - dim, iface, topoKey)); + dim, iface)); - ScalarViewType faceMassMat_("faceMassMat_", numCells, faceCardinality, faceCardinality), + ScalarViewType faceMassMat_("faceMassMat_", 1, faceCardinality, faceCardinality), faceRhsMat_("rhsMat_", numCells, faceCardinality); - FunctionSpaceTools::integrate(faceMassMat_, faceBasisGradAtGradEPoints, wBasisGradAtGradEPoints); - + FunctionSpaceTools::integrate(faceMassMat_, faceBasisGradAtGradEPoints, Kokkos::subview(wBasisGradAtGradEPoints, std::make_pair(0,1), Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()) ); FunctionSpaceTools::integrate(faceRhsMat_, targetGradTanAtTargetGradEPoints, wBasisGradBasisAtTargetGradEPoints); FunctionSpaceTools::integrate(faceRhsMat_, negPartialProjGrad, wBasisGradAtGradEPoints,true); @@ -519,8 +499,8 @@ ProjectionTools::getHGradBasisCoeffs(BasisCoeffsViewType basisCoeffs auto faceDofs = Kokkos::subview(tagToOrdinal, faceDim, iface, Kokkos::ALL()); - ElemSystem faceSystem("faceSystem", false); - faceSystem.solve(basisCoeffs, faceMassMat_, faceRhsMat_, t_, w_, faceDofs, faceCardinality); + ElemSystem faceSystem("faceSystem", true); + faceSystem.solve(refBasisCoeffs, faceMassMat_, faceRhsMat_, t_, w_, faceDofs, faceCardinality); Kokkos::parallel_for("Compute Face Dofs ", Kokkos::RangePolicy (0, faceCardinality), KOKKOS_LAMBDA (const size_t i) { @@ -540,21 +520,21 @@ ProjectionTools::getHGradBasisCoeffs(BasisCoeffsViewType basisCoeffs auto targetGradEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetDerivEvalWeights(dim,0)); auto basisGradEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisDerivEvalWeights(dim,0)); - ScalarViewType cellBasisGradAtGradEPoints("internalBasisGradAtEPoints",numCells,numElemDofs, numGradEPoints, dim); + ScalarViewType cellBasisGradAtGradEPoints("internalBasisGradAtEPoints",1,numElemDofs, numGradEPoints, dim); ScalarViewType negPartialProjGrad("negPartialProjGrad", numCells, numGradEPoints, dim); ScalarViewType wBasisGradAtGradEPoints("wBasisGradAtGradEPoints",numCells,numElemDofs, numGradEPoints,dim); ScalarViewType wBasisGradBasisAtTargetGradEPoints("wBasisGradAtTargetGradEPoints",numCells,numElemDofs, numTargetGradEPoints,dim); auto elemDof = Kokkos::subview(tagToOrdinal, dim, 0, Kokkos::ALL()); - using functorTypeCell_HGRAD = FunctorsProjectionTools::ComputeBasisCoeffsOnCells_HGRAD; - Kokkos::parallel_for(policy, functorTypeCell_HGRAD(basisCoeffs, negPartialProjGrad, cellBasisGradAtGradEPoints, + using functorTypeCell_HGRAD = FunctorsProjectionTools::ComputeBasisCoeffsOnCells_HGRAD; + Kokkos::parallel_for(policy, functorTypeCell_HGRAD(refBasisCoeffs, negPartialProjGrad, cellBasisGradAtGradEPoints, basisGradAtBasisGradEPoints, basisGradEWeights, wBasisGradAtGradEPoints, targetGradEWeights, basisGradAtTargetGradEPoints, wBasisGradBasisAtTargetGradEPoints, computedDofs, elemDof, dim, numElemDofs, offsetBasisGrad, offsetTargetGrad, numVertexDofs+numEdgeDofs+numFaceDofs)); - ScalarViewType cellMassMat_("cellMassMat_", numCells, numElemDofs, numElemDofs), + ScalarViewType cellMassMat_("cellMassMat_", 1, numElemDofs, numElemDofs), cellRhsMat_("rhsMat_", numCells, numElemDofs); - FunctionSpaceTools::integrate(cellMassMat_, cellBasisGradAtGradEPoints, wBasisGradAtGradEPoints); + FunctionSpaceTools::integrate(cellMassMat_, cellBasisGradAtGradEPoints, Kokkos::subview(wBasisGradAtGradEPoints, std::make_pair(0,1), Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()) ); FunctionSpaceTools::integrate(cellRhsMat_, Kokkos::subview(targetGradAtTargetGradEPoints,Kokkos::ALL(),cellTargetGradEPointsRange,Kokkos::ALL()), wBasisGradBasisAtTargetGradEPoints); FunctionSpaceTools::integrate(cellRhsMat_, negPartialProjGrad, wBasisGradAtGradEPoints, true); @@ -564,8 +544,11 @@ ProjectionTools::getHGradBasisCoeffs(BasisCoeffsViewType basisCoeffs auto cellDofs = Kokkos::subview(tagToOrdinal, dim, 0, Kokkos::ALL()); ElemSystem cellSystem("cellSystem", true); - cellSystem.solve(basisCoeffs, cellMassMat_, cellRhsMat_, t_, w_, cellDofs, numElemDofs); + cellSystem.solve(refBasisCoeffs, cellMassMat_, cellRhsMat_, t_, w_, cellDofs, numElemDofs); } + + OrientationTools::modifyBasisByOrientationInverse(basisCoeffs, refBasisCoeffs, orts, cellBasis, true); + } } // Intrepid2 namespace diff --git a/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefL2.hpp b/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefL2.hpp index 890b57d586d3..d4929a6cb10b 100644 --- a/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefL2.hpp +++ b/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefL2.hpp @@ -59,18 +59,18 @@ namespace Intrepid2 { namespace FunctorsProjectionTools { template +typename ViewType4> struct ComputeBasisCoeffsOnVertices_L2 { ViewType1 basisCoeffs_; const ViewType2 tagToOrdinal_; const ViewType3 targetEPointsRange_; const ViewType4 targetAtTargetEPoints_; - const ViewType5 basisAtTargetEPoints_; + const ViewType1 basisAtTargetEPoints_; ordinal_type numVertices_; ComputeBasisCoeffsOnVertices_L2(ViewType1 basisCoeffs, ViewType2 tagToOrdinal, ViewType3 targetEPointsRange, - ViewType4 targetAtTargetEPoints, ViewType5 basisAtTargetEPoints, ordinal_type numVertices) : + ViewType4 targetAtTargetEPoints, ViewType1 basisAtTargetEPoints, ordinal_type numVertices) : basisCoeffs_(basisCoeffs), tagToOrdinal_(tagToOrdinal), targetEPointsRange_(targetEPointsRange), targetAtTargetEPoints_(targetAtTargetEPoints), basisAtTargetEPoints_(basisAtTargetEPoints), numVertices_(numVertices) {} @@ -81,29 +81,29 @@ struct ComputeBasisCoeffsOnVertices_L2 { ordinal_type idof = tagToOrdinal_(0, iv, 0); ordinal_type pt = targetEPointsRange_(0,iv).first; //the value of the basis at the vertex might be different than 1; HGrad basis, so the function is scalar - basisCoeffs_(ic,idof) = targetAtTargetEPoints_(ic,pt)/basisAtTargetEPoints_(ic,idof,pt,0); + basisCoeffs_(ic,idof) = targetAtTargetEPoints_(ic,pt)/basisAtTargetEPoints_(idof,pt,0); } } }; template +typename ViewType4, typename ViewType5> struct ComputeBasisCoeffsOnEdges_L2 { const ViewType1 basisCoeffs_; - const ViewType2 negPartialProj_; - const ViewType2 basisDofDofAtBasisEPoints_; - const ViewType2 basisAtBasisEPoints_; - const ViewType3 basisEWeights_; - const ViewType2 wBasisDofAtBasisEPoints_; - const ViewType3 targetEWeights_; - const ViewType2 basisAtTargetEPoints_; - const ViewType2 wBasisDofAtTargetEPoints_; - const ViewType4 computedDofs_; - const ViewType5 tagToOrdinal_; - const ViewType6 targetAtTargetEPoints_; - const ViewType2 targetTanAtTargetEPoints_; - const ViewType2 refEdgesVec_; + const ViewType1 negPartialProj_; + const ViewType1 basisDofDofAtBasisEPoints_; + const ViewType1 basisAtBasisEPoints_; + const ViewType2 basisEWeights_; + const ViewType1 wBasisDofAtBasisEPoints_; + const ViewType2 targetEWeights_; + const ViewType1 basisAtTargetEPoints_; + const ViewType1 wBasisDofAtTargetEPoints_; + const ViewType3 computedDofs_; + const ViewType4 tagToOrdinal_; + const ViewType5 targetAtTargetEPoints_; + const ViewType1 targetTanAtTargetEPoints_; + const ViewType1 refEdgesVec_; ordinal_type fieldDim_; ordinal_type edgeCardinality_; ordinal_type offsetBasis_; @@ -112,10 +112,10 @@ struct ComputeBasisCoeffsOnEdges_L2 { ordinal_type edgeDim_; ordinal_type iedge_; - ComputeBasisCoeffsOnEdges_L2(const ViewType1 basisCoeffs, ViewType2 negPartialProj, const ViewType2 basisDofDofAtBasisEPoints, - const ViewType2 basisAtBasisEPoints, const ViewType3 basisEWeights, const ViewType2 wBasisDofAtBasisEPoints, const ViewType3 targetEWeights, - const ViewType2 basisAtTargetEPoints, const ViewType2 wBasisDofAtTargetEPoints, const ViewType4 computedDofs, const ViewType5 tagToOrdinal, - const ViewType6 targetAtTargetEPoints, const ViewType2 targetTanAtTargetEPoints, const ViewType2 refEdgesVec, + ComputeBasisCoeffsOnEdges_L2(const ViewType1 basisCoeffs, ViewType1 negPartialProj, const ViewType1 basisDofDofAtBasisEPoints, + const ViewType1 basisAtBasisEPoints, const ViewType2 basisEWeights, const ViewType1 wBasisDofAtBasisEPoints, const ViewType2 targetEWeights, + const ViewType1 basisAtTargetEPoints, const ViewType1 wBasisDofAtTargetEPoints, const ViewType3 computedDofs, const ViewType4 tagToOrdinal, + const ViewType5 targetAtTargetEPoints, const ViewType1 targetTanAtTargetEPoints, const ViewType1 refEdgesVec, ordinal_type fieldDim, ordinal_type edgeCardinality, ordinal_type offsetBasis, ordinal_type offsetTarget, ordinal_type numVertexDofs, ordinal_type edgeDim, ordinal_type iedge) : basisCoeffs_(basisCoeffs), negPartialProj_(negPartialProj), basisDofDofAtBasisEPoints_(basisDofDofAtBasisEPoints), @@ -133,13 +133,15 @@ struct ComputeBasisCoeffsOnEdges_L2 { for(ordinal_type j=0; j +typename ViewType4, typename ViewType5> struct ComputeBasisCoeffsOnFaces_L2 { const ViewType1 basisCoeffs_; - const ViewType2 negPartialProj_; - const ViewType2 faceBasisDofAtBasisEPoints_; - const ViewType2 basisAtBasisEPoints_; - const ViewType3 basisEWeights_; - const ViewType2 wBasisDofAtBasisEPoints_; - const ViewType3 targetEWeights_; - const ViewType2 basisAtTargetEPoints_; - const ViewType2 wBasisDofAtTargetEPoints_; - const ViewType4 computedDofs_; - const ViewType5 tagToOrdinal_; - const ViewType6 orts_; - const ViewType7 targetAtTargetEPoints_; - const ViewType2 targetDofAtTargetEPoints_; - const ViewType8 faceParametrization_; + const ViewType1 negPartialProj_; + const ViewType1 faceBasisDofAtBasisEPoints_; + const ViewType1 basisAtBasisEPoints_; + const ViewType2 basisEWeights_; + const ViewType1 wBasisDofAtBasisEPoints_; + const ViewType2 targetEWeights_; + const ViewType1 basisAtTargetEPoints_; + const ViewType1 wBasisDofAtTargetEPoints_; + const ViewType3 computedDofs_; + const ViewType4 tagToOrdinal_; + const ViewType5 targetAtTargetEPoints_; + const ViewType1 targetDofAtTargetEPoints_; + const ViewType1 refSideNormal_; ordinal_type fieldDim_; ordinal_type faceCardinality_; ordinal_type offsetBasis_; @@ -183,24 +184,23 @@ struct ComputeBasisCoeffsOnFaces_L2 { ordinal_type faceDim_; ordinal_type dim_; ordinal_type iface_; - unsigned topoKey_; bool isHCurlBasis_, isHDivBasis_; - ComputeBasisCoeffsOnFaces_L2(const ViewType1 basisCoeffs, ViewType2 negPartialProj, const ViewType2 faceBasisDofAtBasisEPoints, - const ViewType2 basisAtBasisEPoints, const ViewType3 basisEWeights, const ViewType2 wBasisDofAtBasisEPoints, const ViewType3 targetEWeights, - const ViewType2 basisAtTargetEPoints, const ViewType2 wBasisDofAtTargetEPoints, const ViewType4 computedDofs, const ViewType5 tagToOrdinal, - const ViewType6 orts, const ViewType7 targetAtTargetEPoints, const ViewType2 targetDofAtTargetEPoints, - const ViewType8 faceParametrization, ordinal_type fieldDim, ordinal_type faceCardinality, ordinal_type offsetBasis, + ComputeBasisCoeffsOnFaces_L2(const ViewType1 basisCoeffs, ViewType1 negPartialProj, const ViewType1 faceBasisDofAtBasisEPoints, + const ViewType1 basisAtBasisEPoints, const ViewType2 basisEWeights, const ViewType1 wBasisDofAtBasisEPoints, const ViewType2 targetEWeights, + const ViewType1 basisAtTargetEPoints, const ViewType1 wBasisDofAtTargetEPoints, const ViewType3 computedDofs, const ViewType4 tagToOrdinal, + const ViewType5 targetAtTargetEPoints, const ViewType1 targetDofAtTargetEPoints, + const ViewType1 refSideNormal, ordinal_type fieldDim, ordinal_type faceCardinality, ordinal_type offsetBasis, ordinal_type offsetTarget, ordinal_type numVertexEdgeDofs, ordinal_type numFaces, ordinal_type faceDim, - ordinal_type dim, ordinal_type iface, unsigned topoKey, bool isHCurlBasis, bool isHDivBasis) : + ordinal_type dim, ordinal_type iface, bool isHCurlBasis, bool isHDivBasis) : basisCoeffs_(basisCoeffs), negPartialProj_(negPartialProj), faceBasisDofAtBasisEPoints_(faceBasisDofAtBasisEPoints), basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisDofAtBasisEPoints_(wBasisDofAtBasisEPoints), targetEWeights_(targetEWeights), basisAtTargetEPoints_(basisAtTargetEPoints), wBasisDofAtTargetEPoints_(wBasisDofAtTargetEPoints), - computedDofs_(computedDofs), tagToOrdinal_(tagToOrdinal), orts_(orts), targetAtTargetEPoints_(targetAtTargetEPoints), - targetDofAtTargetEPoints_(targetDofAtTargetEPoints), faceParametrization_(faceParametrization), + computedDofs_(computedDofs), tagToOrdinal_(tagToOrdinal), targetAtTargetEPoints_(targetAtTargetEPoints), + targetDofAtTargetEPoints_(targetDofAtTargetEPoints), refSideNormal_(refSideNormal), fieldDim_(fieldDim), faceCardinality_(faceCardinality), offsetBasis_(offsetBasis), offsetTarget_(offsetTarget), numVertexEdgeDofs_(numVertexEdgeDofs), numFaces_(numFaces), - faceDim_(faceDim), dim_(dim), iface_(iface), topoKey_(topoKey), + faceDim_(faceDim), dim_(dim), iface_(iface), isHCurlBasis_(isHCurlBasis), isHDivBasis_(isHDivBasis) {} @@ -208,16 +208,9 @@ struct ComputeBasisCoeffsOnFaces_L2 { KOKKOS_INLINE_FUNCTION operator()(const ordinal_type ic) const { - ordinal_type fOrt[6]; - orts_(ic).getFaceOrientation(fOrt, numFaces_); - ordinal_type ort = fOrt[iface_]; - if(isHCurlBasis_) { - typename ViewType3::value_type data[3*3]; - auto tangentsAndNormal = ViewType3(data, dim_, dim_); - Impl::OrientationTools::getRefSideTangentsAndNormal(tangentsAndNormal, faceParametrization_,topoKey_, iface_, ort); //normal - typename ViewType2::value_type n[3] = {tangentsAndNormal(2,0), tangentsAndNormal(2,1), tangentsAndNormal(2,2)}; + typename ViewType1::value_type n[3] = {refSideNormal_(0), refSideNormal_(1), refSideNormal_(2)}; for(ordinal_type j=0; j +template struct ComputeBasisCoeffsOnCells_L2 { const ViewType1 basisCoeffs_; - const ViewType2 negPartialProj_; - const ViewType2 internalBasisAtBasisEPoints_; - const ViewType2 basisAtBasisEPoints_; - const ViewType3 basisEWeights_; - const ViewType2 wBasisAtBasisEPoints_; - const ViewType3 targetEWeights_; - const ViewType2 basisAtTargetEPoints_; - const ViewType2 wBasisDofAtTargetEPoints_; - const ViewType4 computedDofs_; - const ViewType5 elemDof_; + const ViewType1 negPartialProj_; + const ViewType1 internalBasisAtBasisEPoints_; + const ViewType1 basisAtBasisEPoints_; + const ViewType2 basisEWeights_; + const ViewType1 wBasisAtBasisEPoints_; + const ViewType2 targetEWeights_; + const ViewType1 basisAtTargetEPoints_; + const ViewType1 wBasisDofAtTargetEPoints_; + const ViewType3 computedDofs_; + const ViewType4 elemDof_; ordinal_type fieldDim_; ordinal_type numElemDofs_; ordinal_type offsetBasis_; ordinal_type offsetTarget_; ordinal_type numVertexEdgeFaceDofs_; - ComputeBasisCoeffsOnCells_L2(const ViewType1 basisCoeffs, ViewType2 negPartialProj, const ViewType2 internalBasisAtBasisEPoints, - const ViewType2 basisAtBasisEPoints, const ViewType3 basisEWeights, const ViewType2 wBasisAtBasisEPoints, const ViewType3 targetEWeights, - const ViewType2 basisAtTargetEPoints, const ViewType2 wBasisDofAtTargetEPoints, const ViewType4 computedDofs, const ViewType5 elemDof, + ComputeBasisCoeffsOnCells_L2(const ViewType1 basisCoeffs, ViewType1 negPartialProj, const ViewType1 internalBasisAtBasisEPoints, + const ViewType1 basisAtBasisEPoints, const ViewType2 basisEWeights, const ViewType1 wBasisAtBasisEPoints, const ViewType2 targetEWeights, + const ViewType1 basisAtTargetEPoints, const ViewType1 wBasisDofAtTargetEPoints, const ViewType3 computedDofs, const ViewType4 elemDof, ordinal_type fieldDim, ordinal_type numElemDofs, ordinal_type offsetBasis, ordinal_type offsetTarget, ordinal_type numVertexEdgeFaceDofs) : basisCoeffs_(basisCoeffs), negPartialProj_(negPartialProj), internalBasisAtBasisEPoints_(internalBasisAtBasisEPoints), basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisAtBasisEPoints_(wBasisAtBasisEPoints), targetEWeights_(targetEWeights), @@ -338,11 +326,11 @@ struct ComputeBasisCoeffsOnCells_L2 { ordinal_type idof = elemDof_(j); for(ordinal_type d=0; d ::getL2BasisCoeffs(Kokkos::DynRankViewgetBasisPointsRange(); - auto refTopologyKey = projStruct->getTopologyKey(); - ordinal_type numTotalBasisEPoints = projStruct->getNumBasisEvalPoints(); auto basisEPoints = projStruct->getAllEvalPoints(EvalPointsType::BASIS); @@ -449,26 +435,16 @@ ProjectionTools::getL2BasisCoeffs(Kokkos::DynRankViewgetAllDofOrdinal()); - ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, fieldDim); - ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints, fieldDim); + ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",basisCardinality, numTotalBasisEPoints, fieldDim); + ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",basisCardinality, numTotalTargetEPoints, fieldDim); { if(fieldDim == 1) { - ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",basisCardinality, numTotalBasisEPoints); - ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",basisCardinality, numTotalTargetEPoints); - cellBasis->getValues(nonOrientedBasisAtTargetEPoints, targetEPoints); - cellBasis->getValues(nonOrientedBasisAtBasisEPoints,basisEPoints); - OrientationTools::modifyBasisByOrientation(Kokkos::subview(basisAtBasisEPoints, Kokkos::ALL(), Kokkos::ALL(), - Kokkos::ALL(),0), nonOrientedBasisAtBasisEPoints, orts, cellBasis); - OrientationTools::modifyBasisByOrientation(Kokkos::subview(basisAtTargetEPoints, Kokkos::ALL(), - Kokkos::ALL(), Kokkos::ALL(),0), nonOrientedBasisAtTargetEPoints, orts, cellBasis); + cellBasis->getValues(Kokkos::subview(basisAtTargetEPoints, Kokkos::ALL(), Kokkos::ALL(), 0), targetEPoints); + cellBasis->getValues(Kokkos::subview(basisAtBasisEPoints, Kokkos::ALL(), Kokkos::ALL(), 0), basisEPoints); } else { - ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",basisCardinality, numTotalBasisEPoints,fieldDim); - ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",basisCardinality, numTotalTargetEPoints,fieldDim); - cellBasis->getValues(nonOrientedBasisAtTargetEPoints, targetEPoints); - cellBasis->getValues(nonOrientedBasisAtBasisEPoints, basisEPoints); - OrientationTools::modifyBasisByOrientation(basisAtBasisEPoints, nonOrientedBasisAtBasisEPoints, orts, cellBasis); - OrientationTools::modifyBasisByOrientation(basisAtTargetEPoints, nonOrientedBasisAtTargetEPoints, orts, cellBasis); + cellBasis->getValues(basisAtTargetEPoints, targetEPoints); + cellBasis->getValues(basisAtBasisEPoints, basisEPoints); } } @@ -500,13 +476,14 @@ ProjectionTools::getL2BasisCoeffs(Kokkos::DynRankView policy(0, numCells); + ScalarViewType refBasisCoeffs("refBasisCoeffs", basisCoeffs.extent(0), basisCoeffs.extent(1)); if(isHGradBasis) { auto targetEPointsRange = Kokkos::create_mirror_view_and_copy(MemSpaceType(), projStruct->getTargetPointsRange()); - using functorType = FunctorsProjectionTools::ComputeBasisCoeffsOnVertices_L2; - Kokkos::parallel_for(policy, functorType(basisCoeffs, tagToOrdinal, targetEPointsRange, + using functorType = FunctorsProjectionTools::ComputeBasisCoeffsOnVertices_L2; + Kokkos::parallel_for(policy, functorType(refBasisCoeffs, tagToOrdinal, targetEPointsRange, targetAtTargetEPoints, basisAtTargetEPoints, numVertices)); } @@ -527,7 +504,7 @@ ProjectionTools::getL2BasisCoeffs(Kokkos::DynRankView::getL2BasisCoeffs(Kokkos::DynRankView; - Kokkos::parallel_for(policy, functorTypeEdge(basisCoeffs, negPartialProj, basisDofAtBasisEPoints, + Kokkos::parallel_for(policy, functorTypeEdge(refBasisCoeffs, negPartialProj, basisDofAtBasisEPoints, basisAtBasisEPoints, basisEWeights, weightedBasisAtBasisEPoints, targetEWeights, basisAtTargetEPoints, weightedBasisAtTargetEPoints, computedDofs, tagToOrdinal, targetAtTargetEPoints,tragetDofAtTargetEPoints, refEdgesVec, fieldDim, edgeCardinality, offsetBasis, offsetTarget, numVertexDofs, edgeDim, ie)); - ScalarViewType edgeMassMat_("edgeMassMat_", numCells, edgeCardinality, edgeCardinality), + ScalarViewType edgeMassMat_("edgeMassMat_", 1, edgeCardinality, edgeCardinality), edgeRhsMat_("rhsMat_", numCells, edgeCardinality); - FunctionSpaceTools::integrate(edgeMassMat_, basisDofAtBasisEPoints, weightedBasisAtBasisEPoints); + FunctionSpaceTools::integrate(edgeMassMat_, basisDofAtBasisEPoints, Kokkos::subview(weightedBasisAtBasisEPoints, std::make_pair(0,1), Kokkos::ALL(),Kokkos::ALL()) ); FunctionSpaceTools::integrate(edgeRhsMat_, tragetDofAtTargetEPoints, weightedBasisAtTargetEPoints); FunctionSpaceTools::integrate(edgeRhsMat_, negPartialProj, weightedBasisAtBasisEPoints,true); @@ -565,22 +542,17 @@ ProjectionTools::getL2BasisCoeffs(Kokkos::DynRankView::ConstViewType subcellParamFace; - if(numFaces>0) - subcellParamFace = RefSubcellParametrization::get(faceDim, cellBasis->getBaseCellTopology().getKey()); - for(ordinal_type iface=0; ifacegetDofCount(faceDim,iface); ordinal_type numTargetEPoints = range_size(targetEPointsRange(faceDim, iface)); ordinal_type numBasisEPoints = range_size(basisEPointsRange(faceDim, iface)); - ScalarViewType faceBasisDofAtBasisEPoints("faceBasisDofAtBasisEPoints",numCells,faceCardinality, numBasisEPoints,faceDofDim); + ScalarViewType faceBasisDofAtBasisEPoints("faceBasisDofAtBasisEPoints",1,faceCardinality, numBasisEPoints,faceDofDim); ScalarViewType wBasisDofAtBasisEPoints("weightedBasisDofAtBasisEPoints",numCells,faceCardinality, numBasisEPoints,faceDofDim); ScalarViewType faceBasisAtTargetEPoints("faceBasisDofAtTargetEPoints",numCells,faceCardinality, numTargetEPoints,faceDofDim); @@ -594,23 +566,25 @@ ProjectionTools::getL2BasisCoeffs(Kokkos::DynRankViewgetTargetEvalWeights(faceDim,iface)); auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(faceDim,iface)); + ScalarViewType refSideNormal("refSideNormal", dim); + CellTools::getReferenceSideNormal(refSideNormal, iface, cellTopo); - using functorTypeFace = FunctorsProjectionTools::ComputeBasisCoeffsOnFaces_L2; + using functorTypeFace = FunctorsProjectionTools::ComputeBasisCoeffsOnFaces_L2; - Kokkos::parallel_for(policy, functorTypeFace(basisCoeffs, negPartialProj,faceBasisDofAtBasisEPoints, + Kokkos::parallel_for(policy, functorTypeFace(refBasisCoeffs, negPartialProj,faceBasisDofAtBasisEPoints, basisAtBasisEPoints, basisEWeights, wBasisDofAtBasisEPoints, targetEWeights, basisAtTargetEPoints, wBasisDofAtTargetEPoints, computedDofs, tagToOrdinal, - orts, targetAtTargetEPoints,targetDofAtTargetEPoints, - subcellParamFace, fieldDim, faceCardinality, offsetBasis, + targetAtTargetEPoints,targetDofAtTargetEPoints, + refSideNormal, fieldDim, faceCardinality, offsetBasis, offsetTarget, numVertexDofs+numEdgeDofs, numFaces, faceDim, - dim, iface, topoKey, isHCurlBasis, isHDivBasis)); + dim, iface, isHCurlBasis, isHDivBasis)); typedef Kokkos::DynRankView WorkArrayViewType; - ScalarViewType faceMassMat_("faceMassMat_", numCells, faceCardinality, faceCardinality), + ScalarViewType faceMassMat_("faceMassMat_", 1, faceCardinality, faceCardinality), faceRhsMat_("rhsMat_", numCells, faceCardinality); - FunctionSpaceTools::integrate(faceMassMat_, faceBasisDofAtBasisEPoints, wBasisDofAtBasisEPoints); + FunctionSpaceTools::integrate(faceMassMat_, faceBasisDofAtBasisEPoints, Kokkos::subview(wBasisDofAtBasisEPoints, std::make_pair(0,1), Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()) ); FunctionSpaceTools::integrate(faceRhsMat_, targetDofAtTargetEPoints, wBasisDofAtTargetEPoints); FunctionSpaceTools::integrate(faceRhsMat_, negPartialProj, wBasisDofAtBasisEPoints,true); @@ -619,8 +593,8 @@ ProjectionTools::getL2BasisCoeffs(Kokkos::DynRankViewgetDofCount(dim,0); @@ -635,7 +609,7 @@ ProjectionTools::getL2BasisCoeffs(Kokkos::DynRankViewgetTargetEvalWeights(dim,0)); @@ -643,20 +617,19 @@ ProjectionTools::getL2BasisCoeffs(Kokkos::DynRankView; - Kokkos::parallel_for(policy, functorType( basisCoeffs, negPartialProj, internalBasisAtBasisEPoints, + using functorType = FunctorsProjectionTools::ComputeBasisCoeffsOnCells_L2; + Kokkos::parallel_for(policy, functorType( refBasisCoeffs, negPartialProj, internalBasisAtBasisEPoints, basisAtBasisEPoints, basisEWeights, wBasisAtBasisEPoints, targetEWeights, basisAtTargetEPoints, wBasisDofAtTargetEPoints, computedDofs, cellDofs, fieldDim, numElemDofs, offsetBasis, offsetTarget, numVertexDofs+numEdgeDofs+numFaceDofs)); typedef Kokkos::DynRankView WorkArrayViewType; - ScalarViewType cellMassMat_("cellMassMat_", numCells, numElemDofs, numElemDofs), + ScalarViewType cellMassMat_("cellMassMat_", 1, numElemDofs, numElemDofs), cellRhsMat_("rhsMat_", numCells, numElemDofs); - FunctionSpaceTools::integrate(cellMassMat_, internalBasisAtBasisEPoints, wBasisAtBasisEPoints); + FunctionSpaceTools::integrate(cellMassMat_, internalBasisAtBasisEPoints, Kokkos::subview(wBasisAtBasisEPoints, std::make_pair(0,1), Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()) ); if(fieldDim==1) FunctionSpaceTools::integrate(cellRhsMat_, Kokkos::subview(targetAtTargetEPoints,Kokkos::ALL(),cellPointsRange,Kokkos::ALL()), Kokkos::subview(wBasisDofAtTargetEPoints,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),0)); @@ -666,9 +639,11 @@ ProjectionTools::getL2BasisCoeffs(Kokkos::DynRankView::modifyBasisByOrientationInverse(basisCoeffs, refBasisCoeffs, orts, cellBasis, true); } @@ -683,84 +658,14 @@ ProjectionTools::getL2DGBasisCoeffs(Kokkos::DynRankView orts, const BasisType* cellBasis, ProjectionStruct * projStruct){ + + Kokkos::DynRankView refBasisCoeffs("refBasisCoeffs", basisCoeffs.extent(0), basisCoeffs.extent(1)); + getL2DGBasisCoeffs(refBasisCoeffs, targetAtTargetEPoints, cellBasis, projStruct); - typedef typename BasisType::scalarType scalarType; - typedef Kokkos::DynRankView ScalarViewType; - const auto cellTopo = cellBasis->getBaseCellTopology(); - - ordinal_type dim = cellTopo.getDimension(); - - auto basisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(), - projStruct->getEvalPoints(dim,0,EvalPointsType::BASIS)); - auto targetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(), - projStruct->getEvalPoints(dim,0,EvalPointsType::TARGET)); - - - ordinal_type numTotalTargetEPoints(targetAtTargetEPoints.extent(1)); - ordinal_type basisCardinality = cellBasis->getCardinality(); - ordinal_type numCells = targetAtTargetEPoints.extent(0); - const ordinal_type fieldDim = (targetAtTargetEPoints.rank()==2) ? 1 : targetAtTargetEPoints.extent(2); - - ordinal_type numTotalBasisEPoints = projStruct->getNumBasisEvalPoints(); - ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, fieldDim); - ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints, fieldDim); - { - if(fieldDim == 1) { - ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",basisCardinality, numTotalBasisEPoints); - ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",basisCardinality, numTotalTargetEPoints); - cellBasis->getValues(nonOrientedBasisAtTargetEPoints, targetEPoints); - cellBasis->getValues(nonOrientedBasisAtBasisEPoints, basisEPoints); - - OrientationTools::modifyBasisByOrientation(Kokkos::subview(basisAtBasisEPoints, Kokkos::ALL(), Kokkos::ALL(), - Kokkos::ALL(),0), nonOrientedBasisAtBasisEPoints, orts, cellBasis); - OrientationTools::modifyBasisByOrientation(Kokkos::subview(basisAtTargetEPoints, Kokkos::ALL(), - Kokkos::ALL(), Kokkos::ALL(),0), nonOrientedBasisAtTargetEPoints, orts, cellBasis); - } - else { - ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",basisCardinality, numTotalBasisEPoints,fieldDim); - ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",basisCardinality, numTotalTargetEPoints,fieldDim); - cellBasis->getValues(nonOrientedBasisAtTargetEPoints, targetEPoints); - cellBasis->getValues(nonOrientedBasisAtBasisEPoints, basisEPoints); - - OrientationTools::modifyBasisByOrientation(basisAtBasisEPoints, nonOrientedBasisAtBasisEPoints, orts, cellBasis); - OrientationTools::modifyBasisByOrientation(basisAtTargetEPoints, nonOrientedBasisAtTargetEPoints, orts, cellBasis); - } - } - - const Kokkos::RangePolicy policy(0, numCells); - ordinal_type numElemDofs = cellBasis->getCardinality(); - - auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(dim,0)); - auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(dim,0)); - - ScalarViewType wBasisAtBasisEPoints("weightedBasisAtBasisEPoints",numCells,numElemDofs, numTotalBasisEPoints,fieldDim); - ScalarViewType wBasisDofAtTargetEPoints("weightedBasisAtTargetEPoints",numCells,numElemDofs, numTotalTargetEPoints,fieldDim); - - using functorType = FunctorsProjectionTools::MultiplyBasisByWeights; - Kokkos::parallel_for( "Multiply basis by weights", policy, functorType(basisAtBasisEPoints, basisEWeights, - wBasisAtBasisEPoints, targetEWeights, basisAtTargetEPoints, wBasisDofAtTargetEPoints, fieldDim, numElemDofs));// )){ - - typedef Kokkos::DynRankView WorkArrayViewType; - ScalarViewType cellMassMat_("cellMassMat_", numCells, numElemDofs, numElemDofs), - cellRhsMat_("rhsMat_", numCells, numElemDofs); - - FunctionSpaceTools::integrate(cellMassMat_, basisAtBasisEPoints, wBasisAtBasisEPoints); - if(fieldDim==1) - FunctionSpaceTools::integrate(cellRhsMat_, targetAtTargetEPoints, - Kokkos::subview(wBasisDofAtTargetEPoints,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),0)); - else - FunctionSpaceTools::integrate(cellRhsMat_, targetAtTargetEPoints, wBasisDofAtTargetEPoints); - - Kokkos::DynRankView hCellDofs("cellDoFs", numElemDofs); - for(ordinal_type i=0; i::modifyBasisByOrientationInverse(basisCoeffs, refBasisCoeffs, orts, cellBasis, true); } + template template void