62 #ifndef __INTREPID2_ORIENTATIONTOOLS_DEF_COEFF_MATRIX_HCURL_HPP__ 63 #define __INTREPID2_ORIENTATIONTOOLS_DEF_COEFF_MATRIX_HCURL_HPP__ 70 #if defined (__clang__) && !defined (__INTEL_COMPILER) 71 #pragma clang system_header 80 #ifdef HAVE_INTREPID2_DEBUG 81 template<
typename subcellBasisType,
82 typename cellBasisType>
85 check_getCoeffMatrix_HCURL(
const subcellBasisType& subcellBasis,
86 const cellBasisType& cellBasis,
87 const ordinal_type subcellId,
88 const ordinal_type subcellOrt) {
89 const shards::CellTopology cellTopo = cellBasis.getBaseCellTopology();
90 const shards::CellTopology subcellTopo = subcellBasis.getBaseCellTopology();
92 const ordinal_type cellDim = cellTopo.getDimension();
93 const ordinal_type subcellDim = subcellTopo.getDimension();
97 INTREPID2_TEST_FOR_EXCEPTION( subcellDim >= cellDim,
99 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): " \
100 "cellDim must be greater than subcellDim.");
102 const auto subcellBaseKey = subcellTopo.getBaseKey();
103 const auto cellBaseKey = cellTopo.getBaseKey();
105 INTREPID2_TEST_FOR_EXCEPTION( subcellBaseKey != shards::Line<>::key &&
106 subcellBaseKey != shards::Quadrilateral<>::key &&
107 subcellBaseKey != shards::Triangle<>::key,
109 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): " \
110 "subcellBasis must have line, quad, or triangle topology.");
112 INTREPID2_TEST_FOR_EXCEPTION( cellBaseKey != shards::Quadrilateral<>::key &&
113 cellBaseKey != shards::Triangle<>::key &&
114 cellBaseKey != shards::Hexahedron<>::key &&
115 cellBaseKey != shards::Wedge<>::key &&
116 cellBaseKey != shards::Tetrahedron<>::key,
118 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): " \
119 "cellBasis must have quad, triangle, hexhedron or tetrahedron topology.");
125 const bool cellBasisIsHCURL = cellBasis.getFunctionSpace() == FUNCTION_SPACE_HCURL;
128 if (cellBasisIsHCURL) {
129 const bool subcellBasisIsHGRAD = subcellBasis.getFunctionSpace() == FUNCTION_SPACE_HGRAD;
130 const bool subcellBasisIsHVOL = subcellBasis.getFunctionSpace() == FUNCTION_SPACE_HVOL;
131 const bool subcellBasisIsHCURL = subcellBasis.getFunctionSpace() == FUNCTION_SPACE_HCURL;
132 const bool cellIsTri = cellBaseKey == shards::Triangle<>::key;
133 const bool cellIsTet = cellBaseKey == shards::Tetrahedron<>::key;
134 const bool cellIsHex = cellBaseKey == shards::Hexahedron<>::key;
135 const bool cellIsQuad = cellBaseKey == shards::Quadrilateral<>::key;
139 switch (subcellDim) {
142 if (cellIsHex || cellIsQuad) {
143 INTREPID2_TEST_FOR_EXCEPTION( !(subcellBasisIsHGRAD||subcellBasisIsHVOL),
145 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): " 146 "subcellBasis function space (1d) is not consistent to cellBasis.");
147 }
else if (cellIsTet || cellIsTri) {
148 INTREPID2_TEST_FOR_EXCEPTION( !subcellBasisIsHVOL,
150 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): " 151 "subcellBasis function space (1d) is not consistent to cellBasis.");
156 INTREPID2_TEST_FOR_EXCEPTION( !subcellBasisIsHCURL,
158 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): " 159 "subcellBasis function space (2d) is not consistent to cellBasis.");
169 template<
typename OutputViewType,
170 typename subcellBasisHostType,
171 typename cellBasisHostType>
176 const subcellBasisHostType& subcellBasis,
177 const cellBasisHostType& cellBasis,
178 const ordinal_type subcellId,
179 const ordinal_type subcellOrt) {
181 #ifdef HAVE_INTREPID2_DEBUG 182 Debug::check_getCoeffMatrix_HCURL(subcellBasis,cellBasis,subcellId,subcellOrt);
185 using value_type =
typename OutputViewType::non_const_value_type;
186 using host_device_type =
typename Kokkos::HostSpace::device_type;
192 const shards::CellTopology cellTopo = cellBasis.getBaseCellTopology();
193 const shards::CellTopology subcellTopo = subcellBasis.getBaseCellTopology();
194 const ordinal_type cellDim = cellTopo.getDimension();
195 const ordinal_type subcellDim = subcellTopo.getDimension();
196 const auto subcellBaseKey = subcellTopo.getBaseKey();
197 const ordinal_type numCellBasis = cellBasis.getCardinality();
198 const ordinal_type numSubcellBasis = subcellBasis.getCardinality();
199 const ordinal_type ndofSubcell = cellBasis.getDofCount(subcellDim,subcellId);
206 Kokkos::DynRankView<value_type,host_device_type> subcellTangents(
"subcellTangents", numSubcellBasis, subcellDim);
207 auto degree = subcellBasis.getDegree();
208 BasisPtr<host_device_type, value_type, value_type> basisPtr;
209 if(subcellBaseKey == shards::Line<>::key) {
211 basisPtr->getDofCoeffs(Kokkos::subview(subcellTangents, Kokkos::ALL(),0));
212 }
else if (subcellBaseKey == shards::Triangle<>::key) {
214 basisPtr->getDofCoeffs(subcellTangents);
215 }
else if (subcellBaseKey == shards::Quadrilateral<>::key) {
217 basisPtr->getDofCoeffs(subcellTangents);
221 Kokkos::DynRankView<value_type,host_device_type> subcellDofCoords(
"subcellDofCoords", basisPtr->getCardinality(), subcellDim);
222 basisPtr->getDofCoords(subcellDofCoords);
223 INTREPID2_TEST_FOR_EXCEPTION( basisPtr->getDofCount(subcellDim,0) != ndofSubcell,
225 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): " \
226 "the number of basisPtr internal DoFs should equate those of the subcell");
229 Kokkos::DynRankView<value_type,host_device_type> refPtsSubcell(
"refPtsSubcell", ndofSubcell, subcellDim);
230 Kokkos::DynRankView<value_type,host_device_type> refSubcellTangents(
"subcellTangents", ndofSubcell, subcellDim);
231 auto tagToOrdinal = basisPtr->getAllDofOrdinal();
232 for (ordinal_type i=0;i<ndofSubcell;++i) {
233 ordinal_type isc = tagToOrdinal(subcellDim, 0, i);
234 for(ordinal_type d=0; d <subcellDim; ++d){
235 refPtsSubcell(i,d) = subcellDofCoords(isc,d);
236 for(ordinal_type k=0; k <subcellDim; ++k)
237 refSubcellTangents(i,d) = subcellTangents(isc,d);
246 Kokkos::DynRankView<value_type,host_device_type> subCellValues(
"subCellValues", numSubcellBasis, ndofSubcell, subcellDim);
248 auto lineValues = Kokkos::subview(subCellValues, Kokkos::ALL(), Kokkos::ALL(), 0);
249 subcellBasis.getValues(lineValues, refPtsSubcell, OPERATOR_VALUE);
251 subcellBasis.getValues(subCellValues, refPtsSubcell, OPERATOR_VALUE);
261 Kokkos::DynRankView<value_type,host_device_type> refPtsCell(
"refPtsCell", ndofSubcell, cellDim);
266 Kokkos::DynRankView<value_type,host_device_type> trJacobianF(
"trJacobianF", subcellDim, cellDim );
272 Kokkos::DynRankView<value_type,host_device_type> cellBasisValues(
"cellBasisValues", numCellBasis, ndofSubcell, cellDim);
273 cellBasis.getValues(cellBasisValues, refPtsCell, OPERATOR_VALUE);
282 Kokkos::DynRankView<value_type,Kokkos::LayoutLeft,host_device_type>
283 PsiMat(
"PsiMat", ndofSubcell, ndofSubcell),
284 PhiMat(
"PhiMat", ndofSubcell, ndofSubcell);
286 auto cellTagToOrdinal = cellBasis.getAllDofOrdinal();
287 auto subcellTagToOrdinal = subcellBasis.getAllDofOrdinal();
289 for (ordinal_type i=0;i<ndofSubcell;++i) {
290 const ordinal_type ic = cellTagToOrdinal(subcellDim, subcellId, i);
291 for (ordinal_type j=0;j<ndofSubcell;++j) {
292 const ordinal_type isc = subcellTagToOrdinal(subcellDim, 0, i);
293 value_type refEntry = 0, ortEntry =0;
294 for (ordinal_type k=0;k<subcellDim;++k) {
295 ortEntry += subCellValues(isc,j,k)*refSubcellTangents(j,k);
296 for (ordinal_type d=0; d<cellDim; ++d)
297 refEntry += cellBasisValues(ic,j,d)*trJacobianF(k,d)*refSubcellTangents(j,k);
299 PsiMat(j,i) = refEntry;
300 PhiMat(j,i) = ortEntry;
306 Teuchos::LAPACK<ordinal_type,value_type> lapack;
307 ordinal_type info = 0;
321 Kokkos::DynRankView<ordinal_type,host_device_type> pivVec(
"pivVec", ndofSubcell);
322 lapack.GESV(ndofSubcell, ndofSubcell,
331 std::stringstream ss;
332 ss <<
">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): " 333 <<
"LAPACK return with error code: " 335 INTREPID2_TEST_FOR_EXCEPTION(
true, std::runtime_error, ss.str().c_str() );
341 const double eps = tolerence();
342 for (ordinal_type i=0;i<ndofSubcell;++i) {
343 auto intmatii = std::round(PhiMat(i,i));
344 PhiMat(i,i) = (std::abs(PhiMat(i,i) - intmatii) < eps) ? intmatii : PhiMat(i,i);
345 for (ordinal_type j=i+1;j<ndofSubcell;++j) {
346 auto matij = PhiMat(i,j);
348 auto intmatji = std::round(PhiMat(j,i));
349 PhiMat(i,j) = (std::abs(PhiMat(j,i) - intmatji) < eps) ? intmatji : PhiMat(j,i);
351 auto intmatij = std::round(matij);
352 PhiMat(j,i) = (std::abs(matij - intmatij) < eps) ? intmatij : matij;
375 const Kokkos::pair<ordinal_type,ordinal_type> range(0, ndofSubcell);
376 auto suboutput = Kokkos::subview(output, range, range);
377 auto tmp = Kokkos::create_mirror_view_and_copy(
typename OutputViewType::device_type::memory_space(), PhiMat);
378 Kokkos::deep_copy(suboutput, tmp);
Implementation of the default H(curl)-compatible FEM basis on Quadrilateral cell. ...
static ConstViewType get(const ordinal_type subcellDim, const unsigned parentCellKey)
Returns a Kokkos view with the coefficients of the parametrization maps for the edges or faces of a r...
Header file for the Intrepid2::Basis_HVOL_LINE_Cn_FEM class.
Header file for the Intrepid2::Basis_HCURL_QUAD_In_FEM class.
Implementation of the locally HVOL-compatible FEM basis of variable order on the [-1,1] reference line cell, using Lagrange polynomials.
Implementation of the default H(curl)-compatible Nedelec (first kind) basis of arbitrary degree on Tr...
Header file for the Intrepid2::Basis_HCURL_TRI_In_FEM class.