Panzer  Version of the Day
Panzer_ProjectToFaces_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Panzer: A partial differential equation assembly
5 // engine for strongly coupled complex multiphysics systems
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39 // Eric C. Cyr (eccyr@sandia.gov)
40 // ***********************************************************************
41 // @HEADER
42 
43 #ifndef PANZER_PROJECT_TO_FACES_IMPL_HPP
44 #define PANZER_PROJECT_TO_FACES_IMPL_HPP
45 
46 #include "Teuchos_Assert.hpp"
47 #include "Phalanx_DataLayout.hpp"
48 
50 #include "Panzer_PureBasis.hpp"
52 #include "Kokkos_ViewFactory.hpp"
53 
54 #include "Teuchos_FancyOStream.hpp"
55 
56 template<typename EvalT,typename Traits>
59  const Teuchos::ParameterList& p)
60 {
61  dof_name = (p.get< std::string >("DOF Name"));
62 
63  if(p.isType< Teuchos::RCP<PureBasis> >("Basis"))
64  basis = p.get< Teuchos::RCP<PureBasis> >("Basis");
65  else
66  basis = p.get< Teuchos::RCP<const PureBasis> >("Basis");
67 
68  quad_degree = 0;
69  if(p.isType<int>("Quadrature Order"))
70  quad_degree = p.get<int>("Quadrature Order");
71 
74 
75  // some sanity checks
77 
78  result = PHX::MDField<ScalarT,Cell,BASIS>(dof_name,basis_layout);
79  this->addEvaluatedField(result);
80 
81  normals = PHX::MDField<ScalarT,Cell,BASIS,Dim>(dof_name+"_Normals",vector_layout);
82  this->addDependentField(normals);
83 
84  if(quad_degree > 0){
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();
90 
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]);
95  }
96 
97  // setup the orientation field
98  std::string orientationFieldName = basis->name() + " Orientation";
99  dof_orientation = PHX::MDField<ScalarT,Cell,NODE>(orientationFieldName,basis_layout);
100  this->addDependentField(dof_orientation);
101 
102  gatherFieldNormals = PHX::MDField<ScalarT,Cell,NODE,Dim>(dof_name+"_Normals",basis->functional_grad);
103  this->addEvaluatedField(gatherFieldNormals);
104 
105  } else {
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]);
109  }
110 
111  this->setName("Project To Faces");
112 }
113 
114 // **********************************************************************
115 template<typename EvalT,typename Traits>
119 {
120  // setup the field data object
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);
125 
126  if(quad_degree > 0){
127  this->utils.setFieldData(dof_orientation,fm);
128  this->utils.setFieldData(gatherFieldNormals,fm);
129  }
130 
131  num_pts = result.dimension(1);
132  num_dim = vector_values[0].dimension(2);
133 
134  TEUCHOS_ASSERT(result.dimension(1) == normals.dimension(1));
135  TEUCHOS_ASSERT(vector_values[0].dimension(2) == normals.dimension(2));
136 }
137 
138 // **********************************************************************
139 template<typename EvalT,typename Traits>
142 {
143 
144  // The coefficients being calculated here in the projection to the face basis
145  // are defined as the integral over the face of the field dotted with the face
146  // normal vector. For a first-order face basis, single point integration is
147  // adequate, so the cubature here just provides the proper weighting.
148  // For higher order, a distinction between "cell" and Gauss points will need
149  // to be made so the field is appropriately projected.
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;
153  PHX::MDField<double,Cell,panzer::NODE,Dim> vertex_coords = workset.cell_vertex_coordinates;
154  int subcell_dim = 2;
155 
156  // One point quadrature if higher order quadrature not requested
157  if (quad_degree == 0){
158 
159  // Collect the reference face information. For now, do nothing with the quadPts.
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);
170  }
171 
172 
173  // Loop over the faces of the workset cells.
174  for (index_t cell = 0; cell < workset.num_cells; ++cell) {
175  for (int p = 0; p < num_pts; ++p) {
176  result(cell,p) = ScalarT(0.0);
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];
180  }
181  }
182 
183  } else {
184 
185  // to compute normals at qps (copied from GatherNormals)
186  int numFaces = gatherFieldNormals.dimension(1);
187  Kokkos::DynRankView<ScalarT,PHX::Device> refFaceTanU = Kokkos::createDynRankView(gatherFieldNormals.get_static_view(),"refFaceTanU",numFaces,num_dim);
188  Kokkos::DynRankView<ScalarT,PHX::Device> refFaceTanV = Kokkos::createDynRankView(gatherFieldNormals.get_static_view(),"refFaceTanV",numFaces,num_dim);
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);
193  for(int d=0;d<num_dim;d++) {
194  refFaceTanU(i,d) = refTanU(d);
195  refFaceTanV(i,d) = refTanV(d);
196  }
197  }
198 
199  // Loop over the faces of the workset cells
200  for (index_t cell = 0; cell < workset.num_cells; ++cell) {
201 
202  // get nodal coordinates for this 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);
207  }
208  }
209 
210  // loop over faces
211  for (int p = 0; p < num_pts; ++p){
212  result(cell,p) = ScalarT(0.0);
213 
214  // get quad weights/pts on reference 2d cell
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);
221 
222  // map 2d quad pts to reference cell (3d)
223  Kokkos::DynRankView<double,PHX::Device> refQuadPts("refQuadPts",faceQuad->getNumPoints(),num_dim);
224  Intrepid2::CellTools<double>::mapToReferenceSubcell(refQuadPts, quadPts, subcell_dim, p, parentCell);
225 
226 
227  // Calculate side jacobian
228  Kokkos::DynRankView<double,PHX::Device> jacobianSide("jacobianSide", 1, faceQuad->getNumPoints(), num_dim, num_dim);
229  Intrepid2::CellTools<double>::setJacobian(jacobianSide, refQuadPts, physicalNodes, parentCell);
230 
231  // Calculate weighted measure at quadrature points
232  Kokkos::DynRankView<double,PHX::Device> weighted_measure("weighted_measure",1,faceQuad->getNumPoints());
233  Intrepid2::FunctionSpaceTools::computeFaceMeasure<double> (weighted_measure, jacobianSide, quadWts, p, parentCell);
234 
235  // loop over quadrature points
236  for (int qp = 0; qp < faceQuad->getNumPoints(); ++qp) {
237 
238  // get normal vector at quad points
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));
248  }
249 
250  // normal = TanU x TanV
251  std::vector<ScalarT> normal(3,0.0);
252  normal[0] = (faceTanU[1]*faceTanV[2] - faceTanU[2]*faceTanV[1])*dof_orientation(cell,p);
253  normal[1] = (faceTanU[2]*faceTanV[0] - faceTanU[0]*faceTanV[2])*dof_orientation(cell,p);
254  normal[2] = (faceTanU[0]*faceTanV[1] - faceTanU[1]*faceTanV[0])*dof_orientation(cell,p);
255 
256  // compute the magnitude of the normal vector
257  ScalarT nnorm(0.0);
258  for(int dim = 0; dim < num_dim; ++dim){
259  nnorm += normal[dim]*normal[dim];
260  }
261  nnorm = std::sqrt(nnorm);
262 
263  // integrate vector dot n
264  // normalize n since jacobian information is factored into both weighted measure and normal
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;
267 
268  }
269  }
270 
271  }
272 
273  }
274 
275 }
276 
277 #endif
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>