Panzer  Version of the Day
Panzer_Dirichlet_Residual_FaceBasis_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_DIRICHLET_RESIDUAL_FACEBASIS_IMPL_HPP
44 #define PANZER_DIRICHLET_RESIDUAL_FACEBASIS_IMPL_HPP
45 
46 #include <cstddef>
47 #include <string>
48 #include <vector>
49 
50 #include "Intrepid2_CellTools.hpp"
51 
53 #include "Kokkos_ViewFactory.hpp"
54 
55 namespace panzer {
56 
57 //**********************************************************************
58 PHX_EVALUATOR_CTOR(DirichletResidual_FaceBasis,p)
59 {
60  std::string residual_name = p.get<std::string>("Residual Name");
61 
64 
65  std::string field_name = p.get<std::string>("DOF Name");
66  std::string dof_name = field_name+"_"+pointRule->getName();
67  std::string value_name = p.get<std::string>("Value Name");
68 
71  Teuchos::RCP<PHX::DataLayout> vector_layout_vector = basis->functional_grad;
72 
73  // some sanity checks
75  TEUCHOS_ASSERT(basis_layout->dimension(0)==vector_layout_dof->dimension(0));
76  TEUCHOS_ASSERT(basis_layout->dimension(1)==vector_layout_dof->dimension(1));
77  TEUCHOS_ASSERT(basis->dimension()==vector_layout_dof->extent_int(2));
78  TEUCHOS_ASSERT(vector_layout_vector->dimension(0)==vector_layout_dof->dimension(0));
79  TEUCHOS_ASSERT(vector_layout_vector->dimension(1)==vector_layout_dof->dimension(1));
80  TEUCHOS_ASSERT(vector_layout_vector->dimension(2)==vector_layout_dof->dimension(2));
81 
82  residual = PHX::MDField<ScalarT,Cell,BASIS>(residual_name, basis_layout);
83  dof = PHX::MDField<const ScalarT,Cell,Point,Dim>(dof_name, vector_layout_dof);
84  value = PHX::MDField<const ScalarT,Cell,Point,Dim>(value_name, vector_layout_vector);
85 
86  // setup the orientation field
87  // std::string orientationFieldName = basis->name() + " Orientation";
88  // dof_orientation = PHX::MDField<ScalarT,Cell,BASIS>(orientationFieldName,
89  // basis_layout);
90 
91  // setup all basis fields that are required
93 
94  // setup all fields to be evaluated and constructed
96 
97  // the field manager will allocate all of these field
98  this->addDependentField(pointValues.jac);
99 
100 
101  this->addEvaluatedField(residual);
102  this->addDependentField(dof);
103  // this->addDependentField(dof_orientation);
104  this->addDependentField(value);
105 
106  std::string n = "Dirichlet Residual Face Basis Evaluator";
107  this->setName(n);
108 }
109 
110 //**********************************************************************
111 PHX_POST_REGISTRATION_SETUP(DirichletResidual_FaceBasis,worksets,fm)
112 {
113  this->utils.setFieldData(residual,fm);
114  this->utils.setFieldData(dof,fm);
115  // this->utils.setFieldData(dof_orientation,fm);
116  this->utils.setFieldData(value,fm);
117  this->utils.setFieldData(pointValues.jac,fm);
118 
119  faceNormal = Kokkos::createDynRankView(residual.get_static_view(),"faceNormal",dof.dimension(0),dof.dimension(1),dof.dimension(2));
120 }
121 
122 //**********************************************************************
123 PHX_EVALUATE_FIELDS(DirichletResidual_FaceBasis,workset)
124 {
125  if(workset.num_cells<=0)
126  return;
127  else {
128  Intrepid2::CellTools<ScalarT>::getPhysicalFaceNormals(faceNormal,
130  this->wda(workset).subcell_index,
131  *basis->getCellTopology());
132 
133  for(index_t c=0;c<workset.num_cells;c++) {
134  for(int b=0;b<dof.extent_int(1);b++) {
135  residual(c,b) = ScalarT(0.0);
136  for(int d=0;d<dof.extent_int(2);d++)
137  residual(c,b) += (dof(c,b,d)-value(c,b,d))*faceNormal(c,b,d);
138  }
139  }
140  }
141 /*
142  if(workset.subcell_dim==1) {
143  Intrepid2::CellTools<ScalarT>::getPhysicalFaceNormals(faceNormal,
144  pointValues.jac,
145  this->wda(workset).subcell_index,
146  *basis->getCellTopology());
147 
148  for(index_t c=0;c<workset.num_cells;c++) {
149  for(int b=0;b<dof.extent_int(1);b++) {
150  residual(c,b) = ScalarT(0.0);
151  for(int d=0;d<dof.extent_int(2);d++)
152  residual(c,b) += (dof(c,b,d)-value(c,b,d))*faceNormal(c,b,d);
153  }
154  }
155  }
156  else if(workset.subcell_dim==2) {
157  // we need to compute the tangents on each edge for each cell.
158  // how do we do this????
159  const shards::CellTopology & parentCell = *basis->getCellTopology();
160  int cellDim = parentCell.getDimension();
161  int numFaces = dof.extent_int(1);
162 
163  refFaceNormal = Kokkos::createDynRankView(residual.get_static_view(),"refFaceNormal",numFaces,cellDim);
164 
165  for(int i=0;i<numFaces;i++) {
166  Kokkos::DynRankView<double,PHX::Device> refFaceNormal_local(cellDim);
167  Intrepid2::CellTools<double>::getReferenceFaceNormal(refFaceNormal_local, i, parentCell);
168 
169  for(int d=0;d<cellDim;d++)
170  refFaceNormal(i,d) = refFaceNormal_local(d);
171  }
172 
173  // Loop over workset faces and edge points
174  for(index_t c=0;c<workset.num_cells;c++) {
175  for(int pt = 0; pt < numFaces; pt++) {
176 
177  // Apply parent cell Jacobian to ref. edge tangent
178  for(int i = 0; i < cellDim; i++) {
179  faceNormal(c, pt, i) = 0.0;
180  for(int j = 0; j < cellDim; j++){
181  faceNormal(c, pt, i) += pointValues.jac(c, pt, i, j)*refFaceNormal(pt,j);
182  }// for j
183  }// for i
184  }// for pt
185  }// for pCell
186 
187  for(index_t c=0;c<workset.num_cells;c++) {
188  for(int b=0;b<dof.extent_int(1);b++) {
189  residual(c,b) = ScalarT(0.0);
190  for(int d=0;d<dof.extent_int(2);d++)
191  residual(c,b) += (dof(c,b,d)-value(c,b,d))*faceNormal(c,b,d);
192  }
193  }
194 
195  }
196  else {
197  // don't know what to do
198  TEUCHOS_ASSERT(false);
199  }
200 */
201 /*
202  // loop over residuals scaling by orientation. This gurantees
203  // everything is oriented in the "positive" direction, this allows
204  // sums acrossed processor to be oriented in the same way (right?)
205  for(index_t c=0;c<workset.num_cells;c++) {
206  for(int b=0;b<dof.extent_int(1);b++) {
207  residual(c,b) *= dof_orientation(c,b);
208  }
209  }
210 */
211 }
212 
213 //**********************************************************************
214 
215 }
216 
217 #endif
PHX::MDField< ScalarT > residual
Evaluates a Dirichlet BC residual corresponding to a field value.
PointValues2< ScalarT, PHX::MDField > pointValues
Interpolates basis DOF values to IP DOF values.
T * get() const
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 > dof
Teuchos::RCP< const panzer::PointRule > pointRule
Teuchos::RCP< PHX::DataLayout > functional_grad
<Cell,Basis,Dim>
int dimension() const
Returns the dimension of the basis from the topology.
PHX_EVALUATOR_CTOR(BasisValues_Evaluator, p)
PHX_EVALUATE_FIELDS(BasisValues_Evaluator, workset)
Teuchos::RCP< PHX::DataLayout > dl_vector
Data layout for vector fields.
Teuchos::RCP< const shards::CellTopology > getCellTopology() const
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
Array< Scalar, Cell, IP, Dim, Dim, void, void, void, void > jac
void setupArrays(const Teuchos::RCP< const panzer::PointRule > &pr, const ArrayFactory &af)
Sizes/allocates memory for arrays.
const std::string & getName() const
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< PHX::DataLayout > functional
<Cell,Basis> or <Cell,Basis>
PHX_POST_REGISTRATION_SETUP(BasisValues_Evaluator, sd, fm)
Kokkos::DynRankView< ScalarT, PHX::Device > faceNormal