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" 52 #include "Kokkos_ViewFactory.hpp" 54 #include "Teuchos_FancyOStream.hpp" 56 template<
typename EvalT,
typename Traits>
61 dof_name = (p.
get< std::string >(
"DOF Name"));
69 if(p.
isType<
int>(
"Quadrature Order"))
70 quad_degree = p.
get<
int>(
"Quadrature Order");
78 result = PHX::MDField<ScalarT,Cell,BASIS>(dof_name,basis_layout);
79 this->addEvaluatedField(result);
81 normals = PHX::MDField<ScalarT,Cell,BASIS,Dim>(dof_name+
"_Normals",vector_layout);
82 this->addDependentField(
normals);
86 Intrepid2::DefaultCubatureFactory<double,Kokkos::DynRankView<double,PHX::Device>,Kokkos::DynRankView<double,PHX::Device>> quadFactory;
88 = quadFactory.create(parentCell.getCellTopologyData(2,0), quad_degree);
89 int numQPoints = quadRule->getNumPoints();
91 vector_values.resize(numQPoints);
92 for(
unsigned qp = 0; qp < numQPoints; ++qp){
93 vector_values[qp] = PHX::MDField<ScalarT,Cell,BASIS,Dim>(dof_name+
"_Vector"+
"_qp_"+std::to_string(qp),vector_layout);
94 this->addDependentField(vector_values[qp]);
98 std::string orientationFieldName =
basis->
name() +
" Orientation";
99 dof_orientation = PHX::MDField<ScalarT,Cell,NODE>(orientationFieldName,basis_layout);
102 gatherFieldNormals = PHX::MDField<ScalarT,Cell,NODE,Dim>(dof_name+
"_Normals",
basis->
functional_grad);
103 this->addEvaluatedField(gatherFieldNormals);
106 vector_values.resize(1);
107 vector_values[0] = PHX::MDField<ScalarT,Cell,BASIS,Dim>(dof_name+
"_Vector",vector_layout);
108 this->addDependentField(vector_values[0]);
111 this->setName(
"Project To Faces");
115 template<
typename EvalT,
typename Traits>
121 this->utils.setFieldData(result,fm);
122 for(
unsigned qp = 0; qp < vector_values.size(); ++qp)
123 this->utils.setFieldData(vector_values[qp],fm);
124 this->utils.setFieldData(
normals,fm);
128 this->utils.setFieldData(gatherFieldNormals,fm);
132 num_dim = vector_values[0].dimension(2);
139 template<
typename EvalT,
typename Traits>
151 Intrepid2::DefaultCubatureFactory<double,Kokkos::DynRankView<double,PHX::Device>,Kokkos::DynRankView<double,PHX::Device>> quadFactory;
157 if (quad_degree == 0){
160 const unsigned num_faces = parentCell.getFaceCount();
161 std::vector<double> refFaceWt(num_faces, 0.0);
162 for (
unsigned f=0; f<num_faces; f++) {
163 faceQuad = quadFactory.create(parentCell.getCellTopologyData(2,f), 1);
164 const int numQPoints = faceQuad->getNumPoints();
165 Kokkos::DynRankView<double,PHX::Device> quadWts(
"quadWts",numQPoints);
166 Kokkos::DynRankView<double,PHX::Device> quadPts(
"quadPts",numQPoints,
num_dim);
167 faceQuad->getCubature(quadPts,quadWts);
168 for (
int q=0; q<numQPoints; q++)
169 refFaceWt[f] += quadWts(q);
174 for (index_t cell = 0; cell < workset.
num_cells; ++cell) {
175 for (
int p = 0; p <
num_pts; ++p) {
177 for (
int dim = 0; dim <
num_dim; ++dim)
178 result(cell,p) += vector_values[0](cell,p,dim) *
normals(cell,p,dim);
179 result(cell,p) *= refFaceWt[p];
186 int numFaces = gatherFieldNormals.dimension(1);
189 for(
int i=0;i<numFaces;i++) {
190 Kokkos::DynRankView<double,PHX::Device> refTanU = Kokkos::DynRankView<double,PHX::Device>(
"refTanU",
num_dim);
191 Kokkos::DynRankView<double,PHX::Device> refTanV = Kokkos::DynRankView<double,PHX::Device>(
"refTanV",
num_dim);
192 Intrepid2::CellTools<double>::getReferenceFaceTangents(refTanU, refTanV, i, parentCell);
194 refFaceTanU(i,d) = refTanU(d);
195 refFaceTanV(i,d) = refTanV(d);
200 for (index_t cell = 0; cell < workset.
num_cells; ++cell) {
203 Kokkos::DynRankView<double,PHX::Device> physicalNodes(
"physicalNodes",1,vertex_coords.dimension(1),
num_dim);
204 for (
int point = 0; point < vertex_coords.dimension(1); ++point){
205 for(
int ict = 0; ict <
num_dim; ict++){
206 physicalNodes(0,point,ict) = vertex_coords(cell,point,ict);
211 for (
int p = 0; p <
num_pts; ++p){
215 const shards::CellTopology & subcell = parentCell.getCellTopologyData(subcell_dim,p);
216 faceQuad = quadFactory.create(subcell, quad_degree);
218 Kokkos::DynRankView<double,PHX::Device> quadWts(
"quadWts",faceQuad->getNumPoints());
219 Kokkos::DynRankView<double,PHX::Device> quadPts(
"quadPts",faceQuad->getNumPoints(),subcell_dim);
220 faceQuad->getCubature(quadPts,quadWts);
223 Kokkos::DynRankView<double,PHX::Device> refQuadPts(
"refQuadPts",faceQuad->getNumPoints(),
num_dim);
224 Intrepid2::CellTools<double>::mapToReferenceSubcell(refQuadPts, quadPts, subcell_dim, p, parentCell);
228 Kokkos::DynRankView<double,PHX::Device> jacobianSide(
"jacobianSide", 1, faceQuad->getNumPoints(),
num_dim,
num_dim);
229 Intrepid2::CellTools<double>::setJacobian(jacobianSide, refQuadPts, physicalNodes, parentCell);
232 Kokkos::DynRankView<double,PHX::Device> weighted_measure(
"weighted_measure",1,faceQuad->getNumPoints());
233 Intrepid2::FunctionSpaceTools::computeFaceMeasure<double> (weighted_measure, jacobianSide, quadWts, p, parentCell);
236 for (
int qp = 0; qp < faceQuad->getNumPoints(); ++qp) {
239 std::vector<double> faceTanU(3);
240 std::vector<double> faceTanV(3);
241 for(
int i = 0; i < 3; i++) {
242 faceTanU[i] = Sacado::ScalarValue<ScalarT>::eval(jacobianSide(0,qp,i,0)*refFaceTanU(p,0))
243 + Sacado::ScalarValue<ScalarT>::eval(jacobianSide(0,qp,i,1)*refFaceTanU(p,1))
244 + Sacado::ScalarValue<ScalarT>::eval(jacobianSide(0,qp,i,2)*refFaceTanU(p,2));
245 faceTanV[i] = Sacado::ScalarValue<ScalarT>::eval(jacobianSide(0,qp,i,0)*refFaceTanV(p,0))
246 + Sacado::ScalarValue<ScalarT>::eval(jacobianSide(0,qp,i,1)*refFaceTanV(p,1))
247 + Sacado::ScalarValue<ScalarT>::eval(jacobianSide(0,qp,i,2)*refFaceTanV(p,2));
251 std::vector<ScalarT>
normal(3,0.0);
258 for(
int dim = 0; dim <
num_dim; ++dim){
261 nnorm = std::sqrt(nnorm);
265 for (
int dim = 0; dim <
num_dim; ++dim)
266 result(cell,p) += weighted_measure(0,qp) * vector_values[qp](cell,p,dim) *
normal[dim] / nnorm;
std::string name() const
A unique key that is the combination of the basis type and basis order.
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)
PHX::MDField< ScalarT > normal
bool isVectorBasis() const
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.
PHX::MDField< const ScalarT, Cell, BASIS > dof_orientation
CellCoordArray cell_vertex_coordinates
Teuchos::RCP< PHX::DataLayout > functional_grad
<Cell,Basis,Dim>
bool isType(const std::string &name) const
PHX::MDField< ScalarT, Cell, Point, Dim > normals
Teuchos::RCP< const shards::CellTopology > getCellTopology() const
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< PHX::DataLayout > functional
<Cell,Basis> or <Cell,Basis>