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
81template<
typename subcellBasisType,
82typename cellBasisType>
85check_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 cannot be smaller 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.");
169template<
typename OutputViewType,
170typename subcellBasisHostType,
171typename 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);
258 Kokkos::DynRankView<value_type,host_device_type> refPtsCell(
"refPtsCell", ndofSubcell, cellDim);
259 Kokkos::DynRankView<value_type,host_device_type> trJacobianF(
"trJacobianF", subcellDim, cellDim );
261 if(cellDim == subcellDim) {
266 Kokkos::DynRankView<value_type,host_device_type> jac(
"data", subcellDim, subcellDim);
268 for(ordinal_type d=0; d<subcellDim; ++d)
269 for(ordinal_type j=0; j<cellDim; ++j) {
270 trJacobianF(j,d) = jac(d,j);
283 Kokkos::DynRankView<value_type,host_device_type> cellBasisValues(
"cellBasisValues", numCellBasis, ndofSubcell, cellDim);
284 cellBasis.getValues(cellBasisValues, refPtsCell, OPERATOR_VALUE);
293 Kokkos::DynRankView<value_type,Kokkos::LayoutLeft,host_device_type>
294 PsiMat(
"PsiMat", ndofSubcell, ndofSubcell),
295 PhiMat(
"PhiMat", ndofSubcell, ndofSubcell);
297 auto cellTagToOrdinal = cellBasis.getAllDofOrdinal();
298 auto subcellTagToOrdinal = subcellBasis.getAllDofOrdinal();
300 for (ordinal_type i=0;i<ndofSubcell;++i) {
301 const ordinal_type ic = cellTagToOrdinal(subcellDim, subcellId, i);
302 for (ordinal_type j=0;j<ndofSubcell;++j) {
303 const ordinal_type isc = subcellTagToOrdinal(subcellDim, 0, i);
304 value_type refEntry = 0, ortEntry =0;
305 for (ordinal_type k=0;k<subcellDim;++k) {
306 ortEntry += subCellValues(isc,j,k)*refSubcellTangents(j,k);
307 for (ordinal_type d=0; d<cellDim; ++d)
308 refEntry += cellBasisValues(ic,j,d)*trJacobianF(k,d)*refSubcellTangents(j,k);
310 PsiMat(j,i) = refEntry;
311 PhiMat(j,i) = ortEntry;
317 Teuchos::LAPACK<ordinal_type,value_type> lapack;
318 ordinal_type info = 0;
332 Kokkos::DynRankView<ordinal_type,host_device_type> pivVec(
"pivVec", ndofSubcell);
333 lapack.GESV(ndofSubcell, ndofSubcell,
342 std::stringstream ss;
343 ss <<
">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): "
344 <<
"LAPACK return with error code: "
346 INTREPID2_TEST_FOR_EXCEPTION(
true, std::runtime_error, ss.str().c_str() );
352 const double eps = tolerence();
353 for (ordinal_type i=0;i<ndofSubcell;++i) {
354 auto intmatii = std::round(PhiMat(i,i));
355 PhiMat(i,i) = (std::abs(PhiMat(i,i) - intmatii) < eps) ? intmatii : PhiMat(i,i);
356 for (ordinal_type j=i+1;j<ndofSubcell;++j) {
357 auto matij = PhiMat(i,j);
359 auto intmatji = std::round(PhiMat(j,i));
360 PhiMat(i,j) = (std::abs(PhiMat(j,i) - intmatji) < eps) ? intmatji : PhiMat(j,i);
362 auto intmatij = std::round(matij);
363 PhiMat(j,i) = (std::abs(matij - intmatij) < eps) ? intmatij : matij;
386 const Kokkos::pair<ordinal_type,ordinal_type> range(0, ndofSubcell);
387 auto suboutput = Kokkos::subview(output, range, range);
388 auto tmp = Kokkos::create_mirror_view_and_copy(
typename OutputViewType::device_type::memory_space(), PhiMat);
389 Kokkos::deep_copy(suboutput, tmp);
Header file for the Intrepid2::Basis_HCURL_QUAD_In_FEM class.
Header file for the Intrepid2::Basis_HCURL_TRI_In_FEM class.
Header file for the Intrepid2::Basis_HVOL_LINE_Cn_FEM class.
Implementation of the default H(curl)-compatible FEM basis on Quadrilateral cell.
Implementation of the default H(curl)-compatible Nedelec (first kind) basis of arbitrary degree on Tr...
Implementation of the locally HVOL-compatible FEM basis of variable order on the [-1,...
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...