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
70template<
typename subcellBasisType,
71typename cellBasisType>
74check_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.");
129template<
typename OutputViewType,
130typename subcellBasisHostType,
131typename cellBasisHostType>
136 const subcellBasisHostType& subcellBasis,
137 const cellBasisHostType& cellBasis,
138 const ordinal_type subcellId,
139 const ordinal_type subcellOrt) {
141#ifdef HAVE_INTREPID2_DEBUG
142 Debug::check_getCoeffMatrix_HGRAD(subcellBasis,cellBasis,subcellId,subcellOrt);
145 using host_device_type =
typename Kokkos::HostSpace::device_type;
146 using value_type =
typename OutputViewType::non_const_value_type;
152 const shards::CellTopology cellTopo = cellBasis.getBaseCellTopology();
153 const shards::CellTopology subcellTopo = subcellBasis.getBaseCellTopology();
154 const ordinal_type cellDim = cellTopo.getDimension();
155 const ordinal_type subcellDim = subcellTopo.getDimension();
156 const auto subcellBaseKey = subcellTopo.getBaseKey();
157 const ordinal_type numCellBasis = cellBasis.getCardinality();
158 const ordinal_type numSubcellBasis = subcellBasis.getCardinality();
159 const ordinal_type ndofSubcell = cellBasis.getDofCount(subcellDim,subcellId);
166 Kokkos::DynRankView<value_type,host_device_type> refPtsSubcell(
"refPtsSubcell", ndofSubcell, subcellDim);
169 INTREPID2_TEST_FOR_EXCEPTION( latticeSize != ndofSubcell,
171 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HGRAD): " \
172 "Lattice size should be equal to the onber of subcell internal DoFs");
176 Kokkos::DynRankView<value_type,host_device_type> refPtsCell(
"refPtsCell", ndofSubcell, cellDim);
178 if(cellDim == subcellDim)
190 Kokkos::DynRankView<value_type,host_device_type> cellBasisValues(
"cellBasisValues", numCellBasis, ndofSubcell);
193 Kokkos::DynRankView<value_type,host_device_type> subcellBasisValues(
"subcellBasisValues", numSubcellBasis, ndofSubcell);
195 cellBasis.getValues(cellBasisValues, refPtsCell, OPERATOR_VALUE);
196 subcellBasis.getValues(subcellBasisValues, refPtsSubcell, OPERATOR_VALUE);
205 Kokkos::DynRankView<value_type,Kokkos::LayoutLeft,host_device_type>
206 PsiMat(
"PsiMat", ndofSubcell, ndofSubcell),
207 PhiMat(
"PhiMat", ndofSubcell, ndofSubcell);
209 auto cellTagToOrdinal = cellBasis.getAllDofOrdinal();
210 auto subcellTagToOrdinal = subcellBasis.getAllDofOrdinal();
212 for (ordinal_type i=0;i<ndofSubcell;++i) {
213 const ordinal_type ic = cellTagToOrdinal(subcellDim, subcellId, i);
214 const ordinal_type isc = subcellTagToOrdinal(subcellDim, 0, i);
215 for (ordinal_type j=0;j<ndofSubcell;++j) {
216 PsiMat(j, i) = cellBasisValues(ic,j);
217 PhiMat(j, i) = subcellBasisValues(isc,j);
223 Teuchos::LAPACK<ordinal_type,value_type> lapack;
224 ordinal_type info = 0;
226 Kokkos::DynRankView<ordinal_type,Kokkos::LayoutLeft,host_device_type> pivVec(
"pivVec", ndofSubcell);
227 lapack.GESV(ndofSubcell, ndofSubcell,
236 std::stringstream ss;
237 ss <<
">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HGRAD): "
238 <<
"LAPACK return with error code: "
240 INTREPID2_TEST_FOR_EXCEPTION(
true, std::runtime_error, ss.str().c_str() );
246 const double eps = tolerence();
247 for (ordinal_type i=0;i<ndofSubcell;++i) {
248 auto intmatii = std::round(PhiMat(i,i));
249 PhiMat(i,i) = (std::abs(PhiMat(i,i) - intmatii) < eps) ? intmatii : PhiMat(i,i);
250 for (ordinal_type j=i+1;j<ndofSubcell;++j) {
251 auto matij = PhiMat(i,j);
253 auto intmatji = std::round(PhiMat(j,i));
254 PhiMat(i,j) = (std::abs(PhiMat(j,i) - intmatji) < eps) ? intmatji : PhiMat(j,i);
256 auto intmatij = std::round(matij);
257 PhiMat(j,i) = (std::abs(matij - intmatij) < eps) ? intmatij : matij;
279 const Kokkos::pair<ordinal_type,ordinal_type> range(0, ndofSubcell);
280 auto suboutput = Kokkos::subview(output, range, range);
281 auto tmp = Kokkos::create_mirror_view_and_copy(
typename OutputViewType::device_type::memory_space(), PhiMat);
282 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...