Panzer  Version of the Day
Panzer_DOFDiv_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_DOF_DIV_IMPL_HPP
44 #define PANZER_DOF_DIV_IMPL_HPP
45 
47 #include "Panzer_BasisIRLayout.hpp"
49 #include "Intrepid2_FunctionSpaceTools.hpp"
50 
51 namespace panzer {
52 
53 namespace {
54 
55 //**********************************************************************
56 template<typename ScalarT,typename ArrayT>
57 void evaluateDiv_withSens(int numCells,
58  PHX::MDField<ScalarT,Cell,IP> & dof_div,
59  PHX::MDField<const ScalarT,Cell,Point> & dof_value,
60  const ArrayT & div_basis)
61 {
62  if(numCells>0) {
63  // evaluate at quadrature points
64 
65  int numFields = div_basis.extent(1);
66  int numPoints = div_basis.extent(2);
67 
68  for (int cell=0; cell<numCells; cell++) {
69  for (int pt=0; pt<numPoints; pt++) {
70  // first initialize to the right thing (prevents over writing with 0)
71  // then loop over one less basis function
72  // ScalarT & div = dof_div(cell,pt);
73  dof_div(cell,pt) = dof_value(cell, 0) * div_basis(cell, 0, pt);
74  for (int bf=1; bf<numFields; bf++)
75  dof_div(cell,pt) += dof_value(cell, bf) * div_basis(cell, bf, pt);
76  }
77  }
78  }
79 }
80 
81 }
82 
83 //**********************************************************************
84 // MOST EVALUATION TYPES
85 //**********************************************************************
86 
87 //**********************************************************************
88 template<typename EvalT, typename TRAITS>
91  use_descriptors_(false),
92  dof_value( p.get<std::string>("Name"),
93  p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->functional),
94  basis_name(p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->name())
95 {
97  = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
98 
99  // Verify that this basis supports the div operation
100  TEUCHOS_TEST_FOR_EXCEPTION(!basis->supportsDiv(),std::logic_error,
101  "DOFDiv: Basis of type \"" << basis->name() << "\" does not support DIV");
102  TEUCHOS_TEST_FOR_EXCEPTION(!basis->requiresOrientations(),std::logic_error,
103  "DOFDiv: Basis of type \"" << basis->name() << "\" in DOF Div should require orientations. So we are throwing.");
104 
105  // build dof_div
106  dof_div = PHX::MDField<ScalarT,Cell,IP>(p.get<std::string>("Div Name"),
108 
109  // add to evaluation graph
110  this->addEvaluatedField(dof_div);
111  this->addDependentField(dof_value);
112 
113  std::string n = "DOFDiv: " + dof_div.fieldTag().name() + " ("+PHX::typeAsString<EvalT>()+")";
114  this->setName(n);
115 }
116 
117 //**********************************************************************
118 template<typename EvalT, typename TRAITS>
120 DOFDiv(const PHX::FieldTag & input,
121  const PHX::FieldTag & output,
122  const panzer::BasisDescriptor & bd,
124  : use_descriptors_(true)
125  , bd_(bd)
126  , id_(id)
127  , dof_value(input)
128 {
129  TEUCHOS_ASSERT(bd.getType()=="HDiv");
130 
131  // build dof_div
132  dof_div = output;
133 
134  // add to evaluation graph
135  this->addEvaluatedField(dof_div);
136  this->addDependentField(dof_value);
137 
138  std::string n = "DOFDiv: " + dof_div.fieldTag().name() + " ("+PHX::typeAsString<EvalT>()+")";
139  this->setName(n);
140 }
141 
142 //**********************************************************************
143 template<typename EvalT, typename TRAITS>
145 postRegistrationSetup(typename TRAITS::SetupData sd,
147 {
148  this->utils.setFieldData(dof_value,fm);
149  this->utils.setFieldData(dof_div,fm);
150 
151  if(not use_descriptors_)
152  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
153 }
154 
155 //**********************************************************************
156 template<typename EvalT, typename TRAITS>
158 evaluateFields(typename TRAITS::EvalData workset)
159 {
160  const panzer::BasisValues2<double> & basisValues = use_descriptors_ ? this->wda(workset).getBasisValues(bd_,id_)
161  : *this->wda(workset).bases[basis_index];
162 
163  evaluateDiv_withSens<ScalarT>(workset.num_cells,dof_div,dof_value,basisValues.div_basis);
164 }
165 
166 //**********************************************************************
167 
168 //**********************************************************************
169 // JACOBIAN EVALUATION TYPES
170 //**********************************************************************
171 
172 //**********************************************************************
173 template<typename TRAITS>
176  use_descriptors_(false),
177  dof_value( p.get<std::string>("Name"),
178  p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->functional),
179  basis_name(p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->name())
180 {
182  = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
183 
184  // do you specialize because you know where the basis functions are and can
185  // skip a large number of AD calculations?
186  if(p.isType<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector")) {
187  offsets = *p.get<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector");
188  accelerate_jacobian = true; // short cut for identity matrix
189  }
190  else
191  accelerate_jacobian = false; // don't short cut for identity matrix
192  accelerate_jacobian = false; // don't short cut for identity matrix
193 
194  // Verify that this basis supports the div operation
195  TEUCHOS_TEST_FOR_EXCEPTION(!basis->supportsDiv(),std::logic_error,
196  "DOFDiv: Basis of type \"" << basis->name() << "\" does not support DIV");
197  TEUCHOS_TEST_FOR_EXCEPTION(!basis->requiresOrientations(),std::logic_error,
198  "DOFDiv: Basis of type \"" << basis->name() << "\" in DOF Div should require orientations. So we are throwing.");
199 
200  // build dof_div
201  dof_div = PHX::MDField<ScalarT,Cell,IP>(p.get<std::string>("Div Name"),
203 
204  // add to evaluation graph
205  this->addEvaluatedField(dof_div);
206  this->addDependentField(dof_value);
207 
208  std::string n = "DOFDiv: " + dof_div.fieldTag().name() + " ("+PHX::typeAsString<panzer::Traits::Jacobian>()+")";
209  this->setName(n);
210 }
211 
212 //**********************************************************************
213 template<typename TRAITS>
215 DOFDiv(const PHX::FieldTag & input,
216  const PHX::FieldTag & output,
217  const panzer::BasisDescriptor & bd,
219  : use_descriptors_(true)
220  , bd_(bd)
221  , id_(id)
222  , dof_value(input)
223 {
224  TEUCHOS_ASSERT(bd.getType()=="HDiv");
225 
226  // build dof_div
227  dof_div = output;
228 
229  accelerate_jacobian = false; // don't short cut for identity matrix
230 
231  // add to evaluation graph
232  this->addEvaluatedField(dof_div);
233  this->addDependentField(dof_value);
234 
235  std::string n = "DOFDiv: " + dof_div.fieldTag().name() + " ("+PHX::typeAsString<panzer::Traits::Jacobian>()+")";
236  this->setName(n);
237 }
238 
239 //**********************************************************************
240 template<typename TRAITS>
242 postRegistrationSetup(typename TRAITS::SetupData sd,
244 {
245  this->utils.setFieldData(dof_value,fm);
246  this->utils.setFieldData(dof_div,fm);
247 
248  if(not use_descriptors_)
249  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
250 }
251 
252 template<typename TRAITS>
254 evaluateFields(typename TRAITS::EvalData workset)
255 {
256  const panzer::BasisValues2<double> & basisValues = use_descriptors_ ? this->wda(workset).getBasisValues(bd_,id_)
257  : *this->wda(workset).bases[basis_index];
258 
259  if(!accelerate_jacobian) {
260  // do the case where we use the AD types to determine the derivatives
261  evaluateDiv_withSens(workset.num_cells,dof_div,dof_value,basisValues.div_basis);
262  return;
263  }
264 
265  TEUCHOS_ASSERT(false);
266 }
267 
268 }
269 
270 #endif
panzer::BasisDescriptor bd_
std::size_t basis_index
panzer::IntegrationDescriptor id_
DOFDiv(const Teuchos::ParameterList &p)
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
int numPoints
int numFields
Array_CellBasisIP div_basis
Teuchos::RCP< PHX::DataLayout > dl_scalar
Data layout for scalar fields.
PHX::MDField< const ScalarT, Cell, Point > dof_value
bool isType(const std::string &name) const
std::vector< std::string >::size_type getBasisIndex(std::string basis_name, const panzer::Workset &workset, WorksetDetailsAccessor &wda)
Returns the index in the workset bases for a particular BasisIRLayout name.
void evaluateFields(typename TRAITS::EvalData d)
#define TEUCHOS_ASSERT(assertion_test)
PHX::MDField< const ScalarT, Cell, Point > dof_value
std::string basis_name
const std::string & getType() const
Get type of basis.
Kokkos::View< const int *, PHX::Device > offsets
PHX::MDField< ScalarT, Cell, IP > dof_div