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>
59 const Teuchos::ParameterList& p)
61 dof_name = (p.get< std::string >(
"DOF Name"));
63 if(p.isType< Teuchos::RCP<PureBasis> >(
"Basis"))
64 basis = p.get< Teuchos::RCP<PureBasis> >(
"Basis");
66 basis = p.get< Teuchos::RCP<const PureBasis> >(
"Basis");
69 if(p.isType<
int>(
"Quadrature Order"))
70 quad_degree = p.get<
int>(
"Quadrature Order");
72 Teuchos::RCP<PHX::DataLayout> basis_layout =
basis->functional;
73 Teuchos::RCP<PHX::DataLayout> vector_layout =
basis->functional_grad;
76 TEUCHOS_ASSERT(
basis->isVectorBasis());
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);
85 const shards::CellTopology & parentCell = *
basis->getCellTopology();
86 Intrepid2::DefaultCubatureFactory<double,Kokkos::DynRankView<double,PHX::Device>,Kokkos::DynRankView<double,PHX::Device>> quadFactory;
87 Teuchos::RCP< Intrepid2::Cubature<double,Kokkos::DynRankView<double,PHX::Device>,Kokkos::DynRankView<double,PHX::Device>> > quadRule
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);
134 TEUCHOS_ASSERT(result.dimension(1) ==
normals.dimension(1));
139 template<
typename EvalT,
typename Traits>
150 const shards::CellTopology & parentCell = *
basis->getCellTopology();
151 Intrepid2::DefaultCubatureFactory<double,Kokkos::DynRankView<double,PHX::Device>,Kokkos::DynRankView<double,PHX::Device>> quadFactory;
152 Teuchos::RCP< Intrepid2::Cubature<double,Kokkos::DynRankView<double,PHX::Device>,Kokkos::DynRankView<double,PHX::Device>> > faceQuad;
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);
217 TEUCHOS_ASSERT(faceQuad->getNumPoints() == vector_values.size());
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;
void evaluateFields(typename Traits::EvalData d)
void postRegistrationSetup(typename Traits::SetupData d, PHX::FieldManager< Traits > &vm)
PHX::MDField< ScalarT > normal
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
PHX::MDField< ScalarT, Cell, Point, Dim > normals
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.