Panzer  Version of the Day
Panzer_Dirichlet_Residual_EdgeBasis_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_EDGEBASIS_IMPL_HPP
44 #define PANZER_DIRICHLET_RESIDUAL_EDGEBASIS_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_EdgeBasis,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(Teuchos::as<unsigned>(basis->dimension())==vector_layout_dof->dimension(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<const 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 Edge Basis Evaluator";
107  this->setName(n);
108 }
109 
110 //**********************************************************************
111 PHX_POST_REGISTRATION_SETUP(DirichletResidual_EdgeBasis,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  edgeTan = Kokkos::createDynRankView(residual.get_static_view(),"edgeTan",dof.dimension(0),dof.dimension(1),dof.dimension(2));
120 }
121 
122 //**********************************************************************
123 PHX_EVALUATE_FIELDS(DirichletResidual_EdgeBasis,workset)
124 {
125  if(workset.num_cells<=0)
126  return;
127 
128  residual.deep_copy(ScalarT(0.0));
129 
130  if(workset.subcell_dim==1) {
131  Intrepid2::CellTools<ScalarT>::getPhysicalEdgeTangents(edgeTan,
133  this->wda(workset).subcell_index,
134  *basis->getCellTopology());
135 
136  for(index_t c=0;c<workset.num_cells;c++) {
137  for(int b=0;b<dof.extent_int(1);b++) {
138  for(int d=0;d<dof.extent_int(2);d++)
139  residual(c,b) += (dof(c,b,d)-value(c,b,d))*edgeTan(c,b,d);
140  }
141  }
142  }
143  else if(workset.subcell_dim==2) {
144  // we need to compute the tangents on each edge for each cell.
145  // how do we do this????
146  const shards::CellTopology & parentCell = *basis->getCellTopology();
147  int cellDim = parentCell.getDimension();
148  int numEdges = dof.extent_int(1);
149 
150  refEdgeTan = Kokkos::createDynRankView(residual.get_static_view(),"refEdgeTan",numEdges,cellDim);
151 
152  for(int i=0;i<numEdges;i++) {
153  Kokkos::DynRankView<double,PHX::Device> refEdgeTan_local("refEdgeTan_local",cellDim);
154  Intrepid2::CellTools<double>::getReferenceEdgeTangent(refEdgeTan_local, i, parentCell);
155 
156  for(int d=0;d<cellDim;d++)
157  refEdgeTan(i,d) = refEdgeTan_local(d);
158  }
159 
160  // Loop over workset faces and edge points
161  for(index_t c=0;c<workset.num_cells;c++) {
162  for(int pt = 0; pt < numEdges; pt++) {
163 
164  // Apply parent cell Jacobian to ref. edge tangent
165  for(int i = 0; i < cellDim; i++) {
166  edgeTan(c, pt, i) = 0.0;
167  for(int j = 0; j < cellDim; j++){
168  edgeTan(c, pt, i) += pointValues.jac(c, pt, i, j)*refEdgeTan(pt,j);
169  }// for j
170  }// for i
171  }// for pt
172  }// for pCell
173 
174  for(index_t c=0;c<workset.num_cells;c++) {
175  for(int b=0;b<dof.extent_int(1);b++) {
176  for(int d=0;d<dof.extent_int(2);d++)
177  residual(c,b) += (dof(c,b,d)-value(c,b,d))*edgeTan(c,b,d);
178  }
179  }
180 
181  }
182  else {
183  // don't know what to do
184  TEUCHOS_ASSERT(false);
185  }
186 
187  // loop over residuals scaling by orientation. This gurantees
188  // everything is oriented in the "positive" direction, this allows
189  // sums acrossed processor to be oriented in the same way (right?)
190  for(index_t c=0;c<workset.num_cells;c++) {
191  for(int b=0;b<dof.extent_int(1);b++) {
192  residual(c,b) *= dof_orientation(c,b);
193  }
194  }
195 }
196 
197 //**********************************************************************
198 
199 }
200 
201 #endif
PHX::MDField< ScalarT > residual
Evaluates a Dirichlet BC residual corresponding to a field value.
std::string name() const
A unique key that is the combination of the basis type and basis order.
PointValues2< ScalarT, PHX::MDField > pointValues
Interpolates basis DOF values to IP DOF values.
T * get() const
bool isVectorBasis() const
Kokkos::DynRankView< ScalarT, PHX::Device > edgeTan
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
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)
Kokkos::DynRankView< ScalarT, PHX::Device > refEdgeTan
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)