43 #ifndef PANZER_PROJECT_TO_FACES_IMPL_HPP
44 #define PANZER_PROJECT_TO_FACES_IMPL_HPP
46 #include "Teuchos_Assert.hpp"
47 #include "Phalanx_DataLayout.hpp"
49 #include "Intrepid2_Cubature.hpp"
50 #include "Intrepid2_DefaultCubatureFactory.hpp"
51 #include "Intrepid2_FunctionSpaceTools.hpp"
52 #include "Intrepid2_Kernels.hpp"
53 #include "Intrepid2_CellTools.hpp"
54 #include "Intrepid2_OrientationTools.hpp"
59 #include "Kokkos_ViewFactory.hpp"
61 #include "Teuchos_FancyOStream.hpp"
65 template<
typename EvalT,
typename Traits>
69 use_fast_method_on_rectangular_hex_mesh(false)
78 if(p.
isType<
int>(
"Quadrature Order"))
81 if(p.
isType<
bool>(
"Use Fast Method for Rectangular Hex Mesh"))
91 this->addEvaluatedField(result);
93 normals = PHX::MDField<const ScalarT,Cell,BASIS,Dim>(dof_name+
"_Normals",vector_layout);
94 this->addDependentField(
normals);
97 const shards::CellTopology & parentCell = *
basis->getCellTopology();
98 Intrepid2::DefaultCubatureFactory quadFactory;
100 = quadFactory.create<PHX::exec_space,double,
double>(parentCell.getCellTopologyData(2,0),
quad_degree);
101 int numQPoints = quadRule->getNumPoints();
104 for (
int qp(0); qp < numQPoints; ++qp)
106 vector_values[qp] = PHX::MDField<const ScalarT,Cell,BASIS,Dim>(dof_name+
"_Vector"+
"_qp_"+std::to_string(qp),vector_layout);
115 TEUCHOS_ASSERT(std::strstr(
basis->getCellTopology()->getBaseName(),
"Hexahedron") !=
nullptr);
118 vector_values[0] = PHX::MDField<const ScalarT,Cell,BASIS,Dim>(dof_name+
"_Vector",vector_layout);
122 if (use_fast_method_on_rectangular_hex_mesh)
123 this->setName(
"Project To Faces (Fast Rectangular Hex)");
125 this->setName(
"Project To Faces");
129 template<
typename EvalT,
typename Traits>
136 num_pts =
result.extent(1);
137 num_dim = vector_values[0].extent(2);
144 template<
typename EvalT,
typename Traits>
155 const shards::CellTopology & parentCell = *basis->getCellTopology();
156 Intrepid2::DefaultCubatureFactory quadFactory;
163 if (use_fast_method_on_rectangular_hex_mesh){
166 const unsigned num_faces = parentCell.getFaceCount();
167 std::vector<double> refFaceWt(num_faces, 0.0);
168 for (
unsigned f=0; f<num_faces; f++) {
169 faceQuad = quadFactory.create<PHX::exec_space,double,
double>(parentCell.getCellTopologyData(2,f), 1);
170 const int numQPoints = faceQuad->getNumPoints();
171 Kokkos::DynRankView<double,PHX::Device> quadWts(
"quadWts",numQPoints);
172 Kokkos::DynRankView<double,PHX::Device> quadPts(
"quadPts",numQPoints,num_dim);
173 faceQuad->getCubature(quadPts,quadWts);
174 for (
int q=0; q<numQPoints; q++)
175 refFaceWt[f] += quadWts(q);
179 for (index_t cell = 0; cell < workset.
num_cells; ++cell) {
180 for (
int p = 0; p < num_pts; ++p) {
182 for (
int dim = 0; dim < num_dim; ++dim)
183 result(cell,p) += vector_values[0](cell,p,dim) * normals(cell,p,dim);
184 result(cell,p) *= refFaceWt[p];
192 int cellDim = parentCell.getDimension();
193 int numFaces = Teuchos::as<int>(parentCell.getFaceCount());
202 for (index_t cell = 0; cell < workset.
num_cells; ++cell) {
205 Kokkos::DynRankView<double,PHX::Device> physicalNodes(
"physicalNodes",1,vertex_coords.extent(1),num_dim);
206 for (
int point(0); point < vertex_coords.extent_int(1); ++point)
208 for (
int ict(0); ict < num_dim; ict++)
209 physicalNodes(0,point,ict) = vertex_coords(cell,point,ict);
212 int faceOrts[6] = {};
213 orientations->at(details.cell_local_ids[cell]).getFaceOrientation(faceOrts, numFaces);
216 for (
int p = 0; p < num_pts; ++p){
219 auto ortEdgeTan_U = Kokkos::subview(refEdges, 0, Kokkos::ALL());
220 auto ortEdgeTan_V = Kokkos::subview(refEdges, 1, Kokkos::ALL());
223 Intrepid2::Orientation::getReferenceFaceTangents(ortEdgeTan_U,
230 const shards::CellTopology & subcell = parentCell.getCellTopologyData(subcell_dim,p);
231 faceQuad = quadFactory.create<PHX::exec_space,double,
double>(subcell, quad_degree);
233 faceQuad->getNumPoints() ==
static_cast<int>(vector_values.size()));
234 Kokkos::DynRankView<double,PHX::Device> quadWts(
"quadWts",faceQuad->getNumPoints());
235 Kokkos::DynRankView<double,PHX::Device> quadPts(
"quadPts",faceQuad->getNumPoints(),subcell_dim);
236 faceQuad->getCubature(quadPts,quadWts);
239 Kokkos::DynRankView<double,PHX::Device> refQuadPts(
"refQuadPts",faceQuad->getNumPoints(),num_dim);
240 Intrepid2::CellTools<PHX::exec_space>::mapToReferenceSubcell(refQuadPts, quadPts, subcell_dim, p, parentCell);
243 Kokkos::DynRankView<double,PHX::Device> jacobianSide(
"jacobianSide", 1, faceQuad->getNumPoints(), num_dim, num_dim);
244 Intrepid2::CellTools<PHX::exec_space>::setJacobian(jacobianSide, refQuadPts, physicalNodes, parentCell);
247 Kokkos::DynRankView<double,PHX::Device> weighted_measure(
"weighted_measure",1,faceQuad->getNumPoints());
248 Kokkos::DynRankView<double,PHX::Device> scratch_space(
"scratch_space",jacobianSide.span());
249 Intrepid2::FunctionSpaceTools<PHX::exec_space>::computeFaceMeasure(weighted_measure, jacobianSide, quadWts, p, parentCell, scratch_space);
252 for (
int qp = 0; qp < faceQuad->getNumPoints(); ++qp) {
254 auto phyEdgeTan_U = Kokkos::subview(phyEdges, 0, Kokkos::ALL());
255 auto phyEdgeTan_V = Kokkos::subview(phyEdges, 1, Kokkos::ALL());
256 auto J = Kokkos::subview(jacobianSide, 0, qp, Kokkos::ALL(), Kokkos::ALL());
258 Intrepid2::Kernels::Serial::matvec_product(phyEdgeTan_U, J, ortEdgeTan_U);
259 Intrepid2::Kernels::Serial::matvec_product(phyEdgeTan_V, J, ortEdgeTan_V);
262 std::vector<ScalarT> normal(3,0.0);
263 normal[0] = (phyEdgeTan_U(1)*phyEdgeTan_V(2) - phyEdgeTan_U(2)*phyEdgeTan_V(1));
264 normal[1] = (phyEdgeTan_U(2)*phyEdgeTan_V(0) - phyEdgeTan_U(0)*phyEdgeTan_V(2));
265 normal[2] = (phyEdgeTan_U(0)*phyEdgeTan_V(1) - phyEdgeTan_U(1)*phyEdgeTan_V(0));
269 for(
int dim = 0; dim < num_dim; ++dim){
270 nnorm += normal[dim]*normal[dim];
272 nnorm = std::sqrt(nnorm);
276 for (
int dim = 0; dim < num_dim; ++dim)
277 result(cell,p) += weighted_measure(0,qp) * vector_values[qp](cell,p,dim) * normal[dim] / nnorm;
Kokkos::DynRankView< typename InputArray::value_type, PHX::Device > createDynRankView(const InputArray &a, const std::string &name, const DimensionPack...dims)
Wrapper to simplify Panzer use of Sacado ViewFactory.
Teuchos::RCP< const std::vector< Intrepid2::Orientation > > orientations_
PHX::MDField< ScalarT, Cell, NODE, Dim > gatherFieldNormals
bool use_fast_method_on_rectangular_hex_mesh
If true, a fast algorithm can be used, but this requires a rectangular/square structured hex mesh...
PHX::MDField< const ScalarT, Cell, BASIS, Dim > normals
void evaluateFields(typename Traits::EvalData d)
T & get(const std::string &name, T def_value)
void postRegistrationSetup(typename Traits::SetupData d, PHX::FieldManager< Traits > &vm)
CellCoordArray cell_vertex_coordinates
PHX::MDField< ScalarT, Cell, BASIS > result
PHX::MDField< ScalarT, panzer::Cell, panzer::IP > result
A field that will be used to build up the result of the integral we're performing.
Teuchos::RCP< const PureBasis > basis
bool isType(const std::string &name) const
#define TEUCHOS_ASSERT(assertion_test)
std::vector< PHX::MDField< const ScalarT, Cell, BASIS, Dim > > vector_values