56 #ifndef __INTREPID2_ORIENTATIONTOOLS_DEF_COEFF_MATRIX_HVOL_HPP__
57 #define __INTREPID2_ORIENTATIONTOOLS_DEF_COEFF_MATRIX_HVOL_HPP__
60 #if defined (__clang__) && !defined (__INTEL_COMPILER)
61 #pragma clang system_header
69 #ifdef HAVE_INTREPID2_DEBUG
70 template<
typename cellBasisType>
73 check_getCoeffMatrix_HVOL(
const cellBasisType& cellBasis,
74 const ordinal_type cellOrt) {
77 const shards::CellTopology cellTopo = cellBasis.getBaseCellTopology();
78 const ordinal_type cellDim = cellTopo.getDimension();
80 INTREPID2_TEST_FOR_EXCEPTION( cellDim > 2,
82 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HVOL): " \
83 "HVOL orientation supported only for (side) cells with dimension less than 3.");
85 const auto cellBaseKey = cellTopo.getBaseKey();
87 INTREPID2_TEST_FOR_EXCEPTION( cellBaseKey != shards::Line<>::key &&
88 cellBaseKey != shards::Quadrilateral<>::key &&
89 cellBaseKey != shards::Triangle<>::key,
91 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HVOL): " \
92 "cellBasis must have line, quad, or triangle topology.");
100 INTREPID2_TEST_FOR_EXCEPTION( cellBasis.getFunctionSpace() != FUNCTION_SPACE_HVOL,
102 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HVOL): " \
103 "cellBasis function space is not HVOL.");
109 template<
typename OutputViewType,
110 typename cellBasisHostType>
115 const cellBasisHostType& cellBasis,
116 const ordinal_type cellOrt,
117 const bool inverse) {
119 #ifdef HAVE_INTREPID2_DEBUG
120 Debug::check_getCoeffMatrix_HVOL(cellBasis,cellOrt);
123 using host_device_type =
typename Kokkos::HostSpace::device_type;
124 using value_type =
typename OutputViewType::non_const_value_type;
130 const shards::CellTopology cellTopo = cellBasis.getBaseCellTopology();
131 const ordinal_type cellDim = cellTopo.getDimension();
132 const auto cellBaseKey = cellTopo.getBaseKey();
133 const ordinal_type cardinality = cellBasis.getCardinality();
140 Kokkos::DynRankView<value_type,host_device_type> refPtsCell(
"refPtsCell", cardinality, cellDim),refPtsCellNotOriented(
"refPtsCellNotOriented", cardinality, cellDim);
142 ordinal_type latticeOffset(1);
145 ordinal_type latticeOrder = (cellTopo.getBaseKey() == shards::Triangle<>::key) ?
146 cellBasis.getDegree() + 3 * latticeOffset :
147 cellBasis.getDegree() + 2 * latticeOffset;
159 Kokkos::DynRankView<value_type,host_device_type> cellBasisValues(
"cellBasisValues", cardinality, cardinality);
162 Kokkos::DynRankView<value_type,host_device_type> nonOrientedBasisValues(
"subcellBasisValues", cardinality, cardinality);
164 cellBasis.getValues(cellBasisValues, refPtsCell, OPERATOR_VALUE);
165 cellBasis.getValues(nonOrientedBasisValues, refPtsCellNotOriented, OPERATOR_VALUE);
174 Kokkos::DynRankView<value_type,Kokkos::LayoutLeft,host_device_type>
175 PsiMat(
"PsiMat", cardinality, cardinality),
176 PhiMat(
"PhiMat", cardinality, cardinality),
180 auto cellTagToOrdinal = cellBasis.getAllDofOrdinal();
182 Kokkos::DynRankView<value_type,host_device_type> jac(
"jacobian",cellDim,cellDim);
184 value_type jacDet(0);
186 jacDet = jac(0,0)*jac(1,1)-jac(0,1)*jac(1,0);
191 for (ordinal_type i=0;i<cardinality;++i) {
192 const ordinal_type ic = cellTagToOrdinal(cellDim, 0, i);
193 for (ordinal_type j=0;j<cardinality;++j) {
194 PsiMat(j, i) = cellBasisValues(ic,j)*jacDet;
195 PhiMat(j, i) = nonOrientedBasisValues(ic,j);
199 RefMat = inverse ? PhiMat : PsiMat;
200 OrtMat = inverse ? PsiMat : PhiMat;
204 Teuchos::LAPACK<ordinal_type,value_type> lapack;
205 ordinal_type info = 0;
207 Kokkos::DynRankView<ordinal_type,Kokkos::LayoutLeft,host_device_type> pivVec(
"pivVec", cardinality);
208 lapack.GESV(cardinality, cardinality,
217 std::stringstream ss;
218 ss <<
">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HVOL): "
219 <<
"LAPACK return with error code: "
221 INTREPID2_TEST_FOR_EXCEPTION(
true, std::runtime_error, ss.str().c_str() );
227 const double eps = tolerence();
228 for (ordinal_type i=0;i<cardinality;++i) {
229 auto intmatii = std::round(OrtMat(i,i));
230 OrtMat(i,i) = (std::abs(OrtMat(i,i) - intmatii) < eps) ? intmatii : OrtMat(i,i);
231 for (ordinal_type j=i+1;j<cardinality;++j) {
232 auto matij = OrtMat(i,j);
234 auto intmatji = std::round(OrtMat(j,i));
235 OrtMat(i,j) = (std::abs(OrtMat(j,i) - intmatji) < eps) ? intmatji : OrtMat(j,i);
237 auto intmatij = std::round(matij);
238 OrtMat(j,i) = (std::abs(matij - intmatij) < eps) ? intmatij : matij;
260 const Kokkos::pair<ordinal_type,ordinal_type> range(0, cardinality);
261 auto suboutput = Kokkos::subview(output, range, range);
262 auto tmp = Kokkos::create_mirror_view_and_copy(
typename OutputViewType::device_type::memory_space(), OrtMat);
263 Kokkos::deep_copy(suboutput, tmp);