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 must be greater than subcellDim.");
91 const auto subcellBaseKey = subcellTopo.getBaseKey();
93 INTREPID2_TEST_FOR_EXCEPTION( subcellBaseKey != shards::Line<>::key &&
94 subcellBaseKey != shards::Quadrilateral<>::key &&
95 subcellBaseKey != shards::Triangle<>::key,
97 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HGRAD): " \
98 "subcellBasis must have line, quad, or triangle topology.");
106 const bool cellBasisIsHGRAD = cellBasis.getFunctionSpace() == FUNCTION_SPACE_HGRAD;
107 const bool subcellBasisIsHGRAD = subcellBasis.getFunctionSpace() == FUNCTION_SPACE_HGRAD;
108 if (cellBasisIsHGRAD) {
109 INTREPID2_TEST_FOR_EXCEPTION( !subcellBasisIsHGRAD,
111 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HGRAD): " \
112 "subcellBasis function space is not consistent to cellBasis.");
115 INTREPID2_TEST_FOR_EXCEPTION( subcellBasis.getDegree() != cellBasis.getDegree(),
117 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HGRAD): " \
118 "subcellBasis has a different polynomial degree from cellBasis' degree.");
124 template<
typename OutputViewType,
125 typename subcellBasisHostType,
126 typename cellBasisHostType>
131 const subcellBasisHostType& subcellBasis,
132 const cellBasisHostType& cellBasis,
133 const ordinal_type subcellId,
134 const ordinal_type subcellOrt) {
136 #ifdef HAVE_INTREPID2_DEBUG 137 Debug::check_getCoeffMatrix_HGRAD(subcellBasis,cellBasis,subcellId,subcellOrt);
140 using host_device_type =
typename Kokkos::HostSpace::device_type;
141 using value_type =
typename OutputViewType::non_const_value_type;
147 const shards::CellTopology cellTopo = cellBasis.getBaseCellTopology();
148 const shards::CellTopology subcellTopo = subcellBasis.getBaseCellTopology();
149 const ordinal_type cellDim = cellTopo.getDimension();
150 const ordinal_type subcellDim = subcellTopo.getDimension();
151 const auto subcellBaseKey = subcellTopo.getBaseKey();
152 const ordinal_type numCellBasis = cellBasis.getCardinality();
153 const ordinal_type numSubcellBasis = subcellBasis.getCardinality();
154 const ordinal_type ndofSubcell = cellBasis.getDofCount(subcellDim,subcellId);
161 Kokkos::DynRankView<value_type,host_device_type> refPtsSubcell(
"refPtsSubcell", ndofSubcell, subcellDim);
164 INTREPID2_TEST_FOR_EXCEPTION( latticeSize != ndofSubcell,
166 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HGRAD): " \
167 "Lattice size should be equal to the onber of subcell internal DoFs");
172 Kokkos::DynRankView<value_type,host_device_type> refPtsCell(
"refPtsCell", ndofSubcell, cellDim);
181 Kokkos::DynRankView<value_type,host_device_type> cellBasisValues(
"cellBasisValues", numCellBasis, ndofSubcell);
184 Kokkos::DynRankView<value_type,host_device_type> subcellBasisValues(
"subcellBasisValues", numSubcellBasis, ndofSubcell);
186 cellBasis.getValues(cellBasisValues, refPtsCell, OPERATOR_VALUE);
187 subcellBasis.getValues(subcellBasisValues, refPtsSubcell, OPERATOR_VALUE);
196 Kokkos::DynRankView<value_type,Kokkos::LayoutLeft,host_device_type>
197 PsiMat(
"PsiMat", ndofSubcell, ndofSubcell),
198 PhiMat(
"PhiMat", ndofSubcell, ndofSubcell);
200 auto cellTagToOrdinal = cellBasis.getAllDofOrdinal();
201 auto subcellTagToOrdinal = subcellBasis.getAllDofOrdinal();
203 for (ordinal_type i=0;i<ndofSubcell;++i) {
204 const ordinal_type ic = cellTagToOrdinal(subcellDim, subcellId, i);
205 const ordinal_type isc = subcellTagToOrdinal(subcellDim, 0, i);
206 for (ordinal_type j=0;j<ndofSubcell;++j) {
207 PsiMat(j, i) = cellBasisValues(ic,j);
208 PhiMat(j, i) = subcellBasisValues(isc,j);
214 Teuchos::LAPACK<ordinal_type,value_type> lapack;
215 ordinal_type info = 0;
217 Kokkos::DynRankView<ordinal_type,Kokkos::LayoutLeft,host_device_type> pivVec(
"pivVec", ndofSubcell);
218 lapack.GESV(ndofSubcell, ndofSubcell,
227 std::stringstream ss;
228 ss <<
">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HGRAD): " 229 <<
"LAPACK return with error code: " 231 INTREPID2_TEST_FOR_EXCEPTION(
true, std::runtime_error, ss.str().c_str() );
237 const double eps = tolerence();
238 for (ordinal_type i=0;i<ndofSubcell;++i) {
239 auto intmatii = std::round(PhiMat(i,i));
240 PhiMat(i,i) = (std::abs(PhiMat(i,i) - intmatii) < eps) ? intmatii : PhiMat(i,i);
241 for (ordinal_type j=i+1;j<ndofSubcell;++j) {
242 auto matij = PhiMat(i,j);
244 auto intmatji = std::round(PhiMat(j,i));
245 PhiMat(i,j) = (std::abs(PhiMat(j,i) - intmatji) < eps) ? intmatji : PhiMat(j,i);
247 auto intmatij = std::round(matij);
248 PhiMat(j,i) = (std::abs(matij - intmatij) < eps) ? intmatij : matij;
270 const Kokkos::pair<ordinal_type,ordinal_type> range(0, ndofSubcell);
271 auto suboutput = Kokkos::subview(output, range, range);
272 auto tmp = Kokkos::create_mirror_view_and_copy(
typename OutputViewType::device_type::memory_space(), PhiMat);
273 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...