56 #ifndef __INTREPID2_ORIENTATIONTOOLS_DEF_COEFF_MATRIX_HGRAD_HPP__
57 #define __INTREPID2_ORIENTATIONTOOLS_DEF_COEFF_MATRIX_HGRAD_HPP__
60 #if defined (__clang__) && !defined (__INTEL_COMPILER)
61 #pragma clang system_header
69 #ifdef HAVE_INTREPID2_DEBUG
70 template<
typename subcellBasisType,
71 typename cellBasisType>
74 check_getCoeffMatrix_HGRAD(
const subcellBasisType& subcellBasis,
75 const cellBasisType& cellBasis,
76 const ordinal_type subcellId,
77 const ordinal_type subcellOrt) {
80 const shards::CellTopology cellTopo = cellBasis.getBaseCellTopology();
81 const shards::CellTopology subcellTopo = subcellBasis.getBaseCellTopology();
83 const ordinal_type cellDim = cellTopo.getDimension();
84 const ordinal_type subcellDim = subcellTopo.getDimension();
86 INTREPID2_TEST_FOR_EXCEPTION( subcellDim > cellDim,
88 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HGRAD): " \
89 "cellDim cannot be smaller than subcellDim.");
91 INTREPID2_TEST_FOR_EXCEPTION( subcellDim > 2,
93 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HGRAD): " \
94 "subCellDim cannot be larger than 2.");
96 const auto subcellBaseKey = subcellTopo.getBaseKey();
98 INTREPID2_TEST_FOR_EXCEPTION( subcellBaseKey != shards::Line<>::key &&
99 subcellBaseKey != shards::Quadrilateral<>::key &&
100 subcellBaseKey != shards::Triangle<>::key,
102 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HGRAD): " \
103 "subcellBasis must have line, quad, or triangle topology.");
111 const bool cellBasisIsHGRAD = cellBasis.getFunctionSpace() == FUNCTION_SPACE_HGRAD;
112 const bool subcellBasisIsHGRAD = subcellBasis.getFunctionSpace() == FUNCTION_SPACE_HGRAD;
113 if (cellBasisIsHGRAD) {
114 INTREPID2_TEST_FOR_EXCEPTION( !subcellBasisIsHGRAD,
116 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HGRAD): " \
117 "subcellBasis function space is not consistent to cellBasis.");
120 INTREPID2_TEST_FOR_EXCEPTION( subcellBasis.getDegree() != cellBasis.getDegree(),
122 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HGRAD): " \
123 "subcellBasis has a different polynomial degree from cellBasis' degree.");
129 template<
typename OutputViewType,
130 typename subcellBasisHostType,
131 typename cellBasisHostType>
136 const subcellBasisHostType& subcellBasis,
137 const cellBasisHostType& cellBasis,
138 const ordinal_type subcellId,
139 const ordinal_type subcellOrt,
140 const bool inverse) {
142 #ifdef HAVE_INTREPID2_DEBUG
143 Debug::check_getCoeffMatrix_HGRAD(subcellBasis,cellBasis,subcellId,subcellOrt);
146 using host_device_type =
typename Kokkos::HostSpace::device_type;
147 using value_type =
typename OutputViewType::non_const_value_type;
153 const shards::CellTopology cellTopo = cellBasis.getBaseCellTopology();
154 const shards::CellTopology subcellTopo = subcellBasis.getBaseCellTopology();
155 const ordinal_type cellDim = cellTopo.getDimension();
156 const ordinal_type subcellDim = subcellTopo.getDimension();
157 const auto subcellBaseKey = subcellTopo.getBaseKey();
158 const ordinal_type numCellBasis = cellBasis.getCardinality();
159 const ordinal_type numSubcellBasis = subcellBasis.getCardinality();
160 const ordinal_type ndofSubcell = cellBasis.getDofCount(subcellDim,subcellId);
167 Kokkos::DynRankView<value_type,host_device_type> refPtsSubcell(
"refPtsSubcell", ndofSubcell, subcellDim);
170 INTREPID2_TEST_FOR_EXCEPTION( latticeSize != ndofSubcell,
172 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HGRAD): " \
173 "Lattice size should be equal to the number of subcell internal DoFs");
177 Kokkos::DynRankView<value_type,host_device_type> refPtsCell(
"refPtsCell", ndofSubcell, cellDim);
179 if(cellDim == subcellDim)
191 Kokkos::DynRankView<value_type,host_device_type> cellBasisValues(
"cellBasisValues", numCellBasis, ndofSubcell);
194 Kokkos::DynRankView<value_type,host_device_type> subcellBasisValues(
"subcellBasisValues", numSubcellBasis, ndofSubcell);
196 cellBasis.getValues(cellBasisValues, refPtsCell, OPERATOR_VALUE);
197 subcellBasis.getValues(subcellBasisValues, refPtsSubcell, OPERATOR_VALUE);
206 Kokkos::DynRankView<value_type,Kokkos::LayoutLeft,host_device_type>
207 PsiMat(
"PsiMat", ndofSubcell, ndofSubcell),
208 PhiMat(
"PhiMat", ndofSubcell, ndofSubcell),
212 auto cellTagToOrdinal = cellBasis.getAllDofOrdinal();
213 auto subcellTagToOrdinal = subcellBasis.getAllDofOrdinal();
215 for (ordinal_type i=0;i<ndofSubcell;++i) {
216 const ordinal_type ic = cellTagToOrdinal(subcellDim, subcellId, i);
217 const ordinal_type isc = subcellTagToOrdinal(subcellDim, 0, i);
218 for (ordinal_type j=0;j<ndofSubcell;++j) {
219 PsiMat(j, i) = cellBasisValues(ic,j);
220 PhiMat(j, i) = subcellBasisValues(isc,j);
224 RefMat = inverse ? PhiMat : PsiMat;
225 OrtMat = inverse ? PsiMat : PhiMat;
229 Teuchos::LAPACK<ordinal_type,value_type> lapack;
230 ordinal_type info = 0;
232 Kokkos::DynRankView<ordinal_type,Kokkos::LayoutLeft,host_device_type> pivVec(
"pivVec", ndofSubcell);
233 lapack.GESV(ndofSubcell, ndofSubcell,
242 std::stringstream ss;
243 ss <<
">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HGRAD): "
244 <<
"LAPACK return with error code: "
246 INTREPID2_TEST_FOR_EXCEPTION(
true, std::runtime_error, ss.str().c_str() );
252 const double eps = tolerence();
253 for (ordinal_type i=0;i<ndofSubcell;++i) {
254 auto intmatii = std::round(OrtMat(i,i));
255 OrtMat(i,i) = (std::abs(OrtMat(i,i) - intmatii) < eps) ? intmatii : OrtMat(i,i);
256 for (ordinal_type j=i+1;j<ndofSubcell;++j) {
257 auto matij = OrtMat(i,j);
259 auto intmatji = std::round(OrtMat(j,i));
260 OrtMat(i,j) = (std::abs(OrtMat(j,i) - intmatji) < eps) ? intmatji : OrtMat(j,i);
262 auto intmatij = std::round(matij);
263 OrtMat(j,i) = (std::abs(matij - intmatij) < eps) ? intmatij : matij;
285 const Kokkos::pair<ordinal_type,ordinal_type> range(0, ndofSubcell);
286 auto suboutput = Kokkos::subview(output, range, range);
287 auto tmp = Kokkos::create_mirror_view_and_copy(
typename OutputViewType::device_type::memory_space(), OrtMat);
288 Kokkos::deep_copy(suboutput, tmp);
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...