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
70template<
typename cellBasisType>
73check_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.");
109template<
typename OutputViewType,
110typename cellBasisHostType>
115 const cellBasisHostType& cellBasis,
116 const ordinal_type cellOrt) {
118#ifdef HAVE_INTREPID2_DEBUG
119 Debug::check_getCoeffMatrix_HVOL(cellBasis,cellOrt);
122 using host_device_type =
typename Kokkos::HostSpace::device_type;
123 using value_type =
typename OutputViewType::non_const_value_type;
129 const shards::CellTopology cellTopo = cellBasis.getBaseCellTopology();
130 const ordinal_type cellDim = cellTopo.getDimension();
131 const auto cellBaseKey = cellTopo.getBaseKey();
132 const ordinal_type cardinality = cellBasis.getCardinality();
139 Kokkos::DynRankView<value_type,host_device_type> refPtsCell(
"refPtsCell", cardinality, cellDim),refPtsCellNotOriented(
"refPtsCellNotOriented", cardinality, cellDim);
141 ordinal_type latticeOffset(1);
144 ordinal_type latticeOrder = (cellTopo.getBaseKey() == shards::Triangle<>::key) ?
145 cellBasis.getDegree() + 3 * latticeOffset :
146 cellBasis.getDegree() + 2 * latticeOffset;
158 Kokkos::DynRankView<value_type,host_device_type> cellBasisValues(
"cellBasisValues", cardinality, cardinality);
161 Kokkos::DynRankView<value_type,host_device_type> nonOrientedBasisValues(
"subcellBasisValues", cardinality, cardinality);
163 cellBasis.getValues(cellBasisValues, refPtsCell, OPERATOR_VALUE);
164 cellBasis.getValues(nonOrientedBasisValues, refPtsCellNotOriented, OPERATOR_VALUE);
173 Kokkos::DynRankView<value_type,Kokkos::LayoutLeft,host_device_type>
174 PsiMat(
"PsiMat", cardinality, cardinality),
175 PhiMat(
"PhiMat", cardinality, cardinality);
177 auto cellTagToOrdinal = cellBasis.getAllDofOrdinal();
179 Kokkos::DynRankView<value_type,host_device_type> jac(
"jacobian",cellDim,cellDim);
181 value_type jacDet(0);
183 jacDet = jac(0,0)*jac(1,1)-jac(0,1)*jac(1,0);
188 for (ordinal_type i=0;i<cardinality;++i) {
189 const ordinal_type ic = cellTagToOrdinal(cellDim, 0, i);
190 for (ordinal_type j=0;j<cardinality;++j) {
191 PsiMat(j, i) = cellBasisValues(ic,j)*jacDet;
192 PhiMat(j, i) = nonOrientedBasisValues(ic,j);
198 Teuchos::LAPACK<ordinal_type,value_type> lapack;
199 ordinal_type info = 0;
201 Kokkos::DynRankView<ordinal_type,Kokkos::LayoutLeft,host_device_type> pivVec(
"pivVec", cardinality);
202 lapack.GESV(cardinality, cardinality,
211 std::stringstream ss;
212 ss <<
">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HVOL): "
213 <<
"LAPACK return with error code: "
215 INTREPID2_TEST_FOR_EXCEPTION(
true, std::runtime_error, ss.str().c_str() );
221 const double eps = tolerence();
222 for (ordinal_type i=0;i<cardinality;++i) {
223 auto intmatii = std::round(PhiMat(i,i));
224 PhiMat(i,i) = (std::abs(PhiMat(i,i) - intmatii) < eps) ? intmatii : PhiMat(i,i);
225 for (ordinal_type j=i+1;j<cardinality;++j) {
226 auto matij = PhiMat(i,j);
228 auto intmatji = std::round(PhiMat(j,i));
229 PhiMat(i,j) = (std::abs(PhiMat(j,i) - intmatji) < eps) ? intmatji : PhiMat(j,i);
231 auto intmatij = std::round(matij);
232 PhiMat(j,i) = (std::abs(matij - intmatij) < eps) ? intmatij : matij;
254 const Kokkos::pair<ordinal_type,ordinal_type> range(0, cardinality);
255 auto suboutput = Kokkos::subview(output, range, range);
256 auto tmp = Kokkos::create_mirror_view_and_copy(
typename OutputViewType::device_type::memory_space(), PhiMat);
257 Kokkos::deep_copy(suboutput, tmp);