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" 63 template<
typename EvalT,
typename Traits>
68 dof_name = (p.
get< std::string >(
"DOF Name"));
76 if(p.
isType<
int>(
"Quadrature Order"))
77 quad_degree = p.
get<
int>(
"Quadrature Order");
85 result = PHX::MDField<ScalarT,Cell,BASIS>(dof_name,basis_layout);
86 this->addEvaluatedField(
result);
88 normals = PHX::MDField<const ScalarT,Cell,BASIS,Dim>(dof_name+
"_Normals",vector_layout);
89 this->addDependentField(normals);
92 const shards::CellTopology & parentCell = *basis->getCellTopology();
93 Intrepid2::DefaultCubatureFactory quadFactory;
95 = quadFactory.create<PHX::exec_space,double,
double>(parentCell.getCellTopologyData(2,0), quad_degree);
96 int numQPoints = quadRule->getNumPoints();
98 vector_values.resize(numQPoints);
99 for (
int qp(0); qp < numQPoints; ++qp)
101 vector_values[qp] = PHX::MDField<const ScalarT,Cell,BASIS,Dim>(dof_name+
"_Vector"+
"_qp_"+std::to_string(qp),vector_layout);
102 this->addDependentField(vector_values[qp]);
105 gatherFieldNormals = PHX::MDField<ScalarT,Cell,NODE,Dim>(dof_name+
"_Normals",basis->functional_grad);
106 this->addEvaluatedField(gatherFieldNormals);
109 vector_values.resize(1);
110 vector_values[0] = PHX::MDField<const ScalarT,Cell,BASIS,Dim>(dof_name+
"_Vector",vector_layout);
111 this->addDependentField(vector_values[0]);
114 this->setName(
"Project To Faces");
118 template<
typename EvalT,
typename Traits>
125 num_pts =
result.extent(1);
126 num_dim = vector_values[0].extent(2);
133 template<
typename EvalT,
typename Traits>
144 const shards::CellTopology & parentCell = *basis->getCellTopology();
145 Intrepid2::DefaultCubatureFactory quadFactory;
149 if (quad_degree == 0){
152 const unsigned num_faces = parentCell.getFaceCount();
153 std::vector<double> refFaceWt(num_faces, 0.0);
154 for (
unsigned f=0; f<num_faces; f++) {
155 faceQuad = quadFactory.create<PHX::exec_space,double,
double>(parentCell.getCellTopologyData(2,f), 1);
156 const int numQPoints = faceQuad->getNumPoints();
157 Kokkos::DynRankView<double,PHX::Device> quadWts(
"quadWts",numQPoints);
158 Kokkos::DynRankView<double,PHX::Device> quadPts(
"quadPts",numQPoints,num_dim);
159 faceQuad->getCubature(quadPts,quadWts);
160 for (
int q=0; q<numQPoints; q++)
161 refFaceWt[f] += quadWts(q);
166 for (index_t cell = 0; cell < workset.
num_cells; ++cell) {
167 for (
int p = 0; p < num_pts; ++p) {
169 for (
int dim = 0; dim < num_dim; ++dim)
170 result(cell,p) += vector_values[0](cell,p,dim) * normals(cell,p,dim);
171 result(cell,p) *= refFaceWt[p];
179 int cellDim = parentCell.getDimension();
180 int numFaces = Teuchos::as<int>(parentCell.getFaceCount());
189 for (index_t cell = 0; cell < workset.
num_cells; ++cell) {
192 Kokkos::DynRankView<double,PHX::Device> physicalNodes(
"physicalNodes",1,vertex_coords.extent(1),num_dim);
193 for (
int point(0); point < vertex_coords.extent_int(1); ++point)
195 for (
int ict(0); ict < num_dim; ict++)
196 physicalNodes(0,point,ict) = vertex_coords(cell,point,ict);
199 int faceOrts[6] = {};
200 orientations->at(
details.cell_local_ids[cell]).getFaceOrientation(faceOrts, numFaces);
203 for (
int p = 0; p < num_pts; ++p){
206 auto ortEdgeTan_U = Kokkos::subview(refEdges, 0, Kokkos::ALL());
207 auto ortEdgeTan_V = Kokkos::subview(refEdges, 1, Kokkos::ALL());
210 Intrepid2::Orientation::getReferenceFaceTangents(ortEdgeTan_U,
217 const shards::CellTopology & subcell = parentCell.getCellTopologyData(subcell_dim,p);
218 faceQuad = quadFactory.create<PHX::exec_space,double,
double>(subcell, quad_degree);
220 faceQuad->getNumPoints() ==
static_cast<int>(vector_values.size()));
221 Kokkos::DynRankView<double,PHX::Device> quadWts(
"quadWts",faceQuad->getNumPoints());
222 Kokkos::DynRankView<double,PHX::Device> quadPts(
"quadPts",faceQuad->getNumPoints(),subcell_dim);
223 faceQuad->getCubature(quadPts,quadWts);
226 Kokkos::DynRankView<double,PHX::Device> refQuadPts(
"refQuadPts",faceQuad->getNumPoints(),num_dim);
227 Intrepid2::CellTools<PHX::exec_space>::mapToReferenceSubcell(refQuadPts, quadPts, subcell_dim, p, parentCell);
230 Kokkos::DynRankView<double,PHX::Device> jacobianSide(
"jacobianSide", 1, faceQuad->getNumPoints(), num_dim, num_dim);
231 Intrepid2::CellTools<PHX::exec_space>::setJacobian(jacobianSide, refQuadPts, physicalNodes, parentCell);
234 Kokkos::DynRankView<double,PHX::Device> weighted_measure(
"weighted_measure",1,faceQuad->getNumPoints());
235 Kokkos::DynRankView<double,PHX::Device> scratch_space(
"scratch_space",jacobianSide.span());
236 Intrepid2::FunctionSpaceTools<PHX::exec_space>::computeFaceMeasure(weighted_measure, jacobianSide, quadWts, p, parentCell, scratch_space);
239 for (
int qp = 0; qp < faceQuad->getNumPoints(); ++qp) {
241 auto phyEdgeTan_U = Kokkos::subview(phyEdges, 0, Kokkos::ALL());
242 auto phyEdgeTan_V = Kokkos::subview(phyEdges, 1, Kokkos::ALL());
243 auto J = Kokkos::subview(jacobianSide, 0, qp, Kokkos::ALL(), Kokkos::ALL());
245 Intrepid2::Kernels::Serial::matvec_product(phyEdgeTan_U, J, ortEdgeTan_U);
246 Intrepid2::Kernels::Serial::matvec_product(phyEdgeTan_V, J, ortEdgeTan_V);
249 std::vector<ScalarT> normal(3,0.0);
250 normal[0] = (phyEdgeTan_U(1)*phyEdgeTan_V(2) - phyEdgeTan_U(2)*phyEdgeTan_V(1));
251 normal[1] = (phyEdgeTan_U(2)*phyEdgeTan_V(0) - phyEdgeTan_U(0)*phyEdgeTan_V(2));
252 normal[2] = (phyEdgeTan_U(0)*phyEdgeTan_V(1) - phyEdgeTan_U(1)*phyEdgeTan_V(0));
256 for(
int dim = 0; dim < num_dim; ++dim){
257 nnorm += normal[dim]*normal[dim];
259 nnorm = std::sqrt(nnorm);
263 for (
int dim = 0; dim < num_dim; ++dim)
264 result(cell,p) += weighted_measure(0,qp) * vector_values[qp](cell,p,dim) * normal[dim] / nnorm;
Teuchos::RCP< const std::vector< Intrepid2::Orientation > > orientations_
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, panzer::Cell, panzer::IP > result
A field that will be used to build up the result of the integral we're performing.
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.
bool isType(const std::string &name) const
#define TEUCHOS_ASSERT(assertion_test)