Panzer  Version of the Day
Panzer_DOF_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_IMPL_HPP
44 #define PANZER_DOF_IMPL_HPP
45 
46 #include <algorithm>
48 #include "Panzer_BasisIRLayout.hpp"
51 #include "Panzer_DOF_Functors.hpp"
52 
53 #include "Intrepid2_FunctionSpaceTools.hpp"
54 
55 namespace panzer {
56 
57 //**********************************************************************
58 //* DOF evaluator
59 //**********************************************************************
60 
61 //**********************************************************************
62 // MOST EVALUATION TYPES
63 //**********************************************************************
64 
65 //**********************************************************************
66 template<typename EvalT, typename TRAITS>
69  use_descriptors_(false),
70  dof_basis( p.get<std::string>("Name"),
71  p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->functional),
72  basis_name(p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->name())
73 {
75  = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
76  is_vector_basis = basis->isVectorBasis();
77 
78  // swap between scalar basis value, or vector basis value
79  if(basis->isScalarBasis()) {
80  dof_ip_scalar = PHX::MDField<ScalarT,Cell,Point>(
81  p.get<std::string>("Name"),
83  this->addEvaluatedField(dof_ip_scalar);
84  }
85  else if(basis->isVectorBasis()) {
86  dof_ip_vector = PHX::MDField<ScalarT,Cell,Point,Dim>(
87  p.get<std::string>("Name"),
89  this->addEvaluatedField(dof_ip_vector);
90  }
91  else
92  { TEUCHOS_ASSERT(false); }
93 
94  this->addDependentField(dof_basis);
95 
96  std::string n = "DOF: " + dof_basis.fieldTag().name() + " ("+PHX::typeAsString<EvalT>()+")";
97  this->setName(n);
98 }
99 
100 //**********************************************************************
101 template<typename EvalT, typename TRAITS>
103 DOF(const PHX::FieldTag & input,
104  const PHX::FieldTag & output,
105  const panzer::BasisDescriptor & bd,
107  : use_descriptors_(true)
108  , bd_(bd)
109  , id_(id)
110  , dof_basis(input)
111 {
112  TEUCHOS_ASSERT(bd.getType()=="HGrad" || bd.getType()=="HCurl" ||
113  bd.getType()=="HDiv" || bd.getType()=="Const")
114 
115  is_vector_basis = (bd.getType()=="HCurl" || bd.getType()=="HDiv");
116 
117  // swap between scalar basis value, or vector basis value
118  if(not is_vector_basis) {
119  dof_ip_scalar = output;
120  this->addEvaluatedField(dof_ip_scalar);
121  }
122  else {
123  dof_ip_vector = output;
124  this->addEvaluatedField(dof_ip_vector);
125  }
126 
127  this->addDependentField(dof_basis);
128 
129  std::string n = "DOF: " + dof_basis.fieldTag().name() + " ("+PHX::typeAsString<EvalT>()+")";
130  this->setName(n);
131 }
132 
133 //**********************************************************************
134 template<typename EvalT, typename TRAITS>
136 postRegistrationSetup(typename TRAITS::SetupData sd,
138 {
139  this->utils.setFieldData(dof_basis,fm);
140  if(is_vector_basis)
141  this->utils.setFieldData(dof_ip_vector,fm);
142  else
143  this->utils.setFieldData(dof_ip_scalar,fm);
144 
145  // descriptors don't access the basis values in the same way
146  if(not use_descriptors_)
147  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
148 }
149 
150 //**********************************************************************
151 template<typename EvalT, typename TRAITS>
153 evaluateFields(typename TRAITS::EvalData workset)
154 {
155  const panzer::BasisValues2<double> & basisValues = use_descriptors_ ? this->wda(workset).getBasisValues(bd_,id_)
156  : *this->wda(workset).bases[basis_index];
157 
158  if(is_vector_basis) {
159  int spaceDim = basisValues.basis_vector.extent(3);
160  if(spaceDim==3) {
162  Kokkos::parallel_for(workset.num_cells,functor);
163  }
164  else {
166  Kokkos::parallel_for(workset.num_cells,functor);
167  }
168  }
169  else {
171  Kokkos::parallel_for(workset.num_cells,functor);
172  }
173 }
174 
175 //**********************************************************************
176 
177 //**********************************************************************
178 // JACOBIAN EVALUATION TYPES
179 //**********************************************************************
180 
181 //**********************************************************************
182 template<typename TRAITS>
185  use_descriptors_(false),
186  dof_basis( p.get<std::string>("Name"),
187  p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->functional),
188  basis_name(p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->name())
189 {
191  = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
192  is_vector_basis = basis->isVectorBasis();
193 
194  if(p.isType<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector")) {
195  const std::vector<int> & offsets = *p.get<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector");
196 
197  // allocate and copy offsets vector to Kokkos array
198  offsets_array = Kokkos::View<int*,PHX::Device>("offsets",offsets.size());
199  for(std::size_t i=0;i<offsets.size();i++)
200  offsets_array(i) = offsets[i];
201 
202  accelerate_jacobian_enabled = true; // short cut for identity matrix
203 
204  // get the sensitivities name that is valid for accelerated jacobians
205  sensitivities_name = true;
206  if (p.isType<std::string>("Sensitivities Name"))
207  sensitivities_name = p.get<std::string>("Sensitivities Name");
208  }
209  else
210  accelerate_jacobian_enabled = false; // don't short cut for identity matrix
211 
212  // swap between scalar basis value, or vector basis value
213  if(basis->isScalarBasis()) {
214  dof_ip_scalar = PHX::MDField<ScalarT,Cell,Point>(
215  p.get<std::string>("Name"),
217  this->addEvaluatedField(dof_ip_scalar);
218  }
219  else if(basis->isVectorBasis()) {
220  dof_ip_vector = PHX::MDField<ScalarT,Cell,Point,Dim>(
221  p.get<std::string>("Name"),
223  this->addEvaluatedField(dof_ip_vector);
224  }
225  else
226  { TEUCHOS_ASSERT(false); }
227 
228  this->addDependentField(dof_basis);
229 
230  std::string n = "DOF: " + dof_basis.fieldTag().name()
231  + ( accelerate_jacobian_enabled ? " accel_jac " : "slow_jac" )
232  + " ("+PHX::typeAsString<panzer::Traits::Jacobian>()+")";
233  this->setName(n);
234 }
235 
236 //**********************************************************************
237 template<typename TRAITS>
239 DOF(const PHX::FieldTag & input,
240  const PHX::FieldTag & output,
241  const panzer::BasisDescriptor & bd,
243  : use_descriptors_(true)
244  , bd_(bd)
245  , id_(id)
246  , dof_basis(input)
247 {
248  TEUCHOS_ASSERT(bd.getType()=="HGrad" || bd.getType()=="HCurl" ||
249  bd.getType()=="HDiv" || bd.getType()=="Const")
250 
251  accelerate_jacobian_enabled = false; // don't short cut for identity matrix
252 
253  is_vector_basis = (bd.getType()=="HCurl" || bd.getType()=="HDiv");
254 
255  // swap between scalar basis value, or vector basis value
256  if(not is_vector_basis) {
257  dof_ip_scalar = output;
258  this->addEvaluatedField(dof_ip_scalar);
259  }
260  else {
261  dof_ip_vector = output;
262  this->addEvaluatedField(dof_ip_vector);
263  }
264 
265  this->addDependentField(dof_basis);
266 
267  std::string n = "DOF: " + dof_basis.fieldTag().name() + " slow_jac(descriptor) ("+PHX::typeAsString<typename TRAITS::Jacobian>()+")";
268  this->setName(n);
269 }
270 
271 //**********************************************************************
272 template<typename TRAITS>
274 postRegistrationSetup(typename TRAITS::SetupData sd,
276 {
277  this->utils.setFieldData(dof_basis,fm);
278  if(is_vector_basis)
279  this->utils.setFieldData(dof_ip_vector,fm);
280  else
281  this->utils.setFieldData(dof_ip_scalar,fm);
282 
283  // descriptors don't access the basis values in the same way
284  if(not use_descriptors_)
285  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
286 }
287 
288 // **********************************************************************
289 template<typename TRAITS>
291 preEvaluate(typename TRAITS::PreEvalData d)
292 {
293  // if sensitivities were requrested for this field enable accelerated
294  // jacobian calculations
295  accelerate_jacobian = false;
296  if(accelerate_jacobian_enabled && d.first_sensitivities_name==sensitivities_name) {
297  accelerate_jacobian = true;
298  }
299 }
300 
301 //**********************************************************************
302 template<typename TRAITS>
304 evaluateFields(typename TRAITS::EvalData workset)
305 {
306  const panzer::BasisValues2<double> & basisValues = use_descriptors_ ? this->wda(workset).getBasisValues(bd_,id_)
307  : *this->wda(workset).bases[basis_index];
308 
309  if(is_vector_basis) {
310  if(accelerate_jacobian) {
311  int spaceDim = basisValues.basis_vector.extent(3);
312  if(spaceDim==3) {
314  Kokkos::parallel_for(workset.num_cells,functor);
315  }
316  else {
318  Kokkos::parallel_for(workset.num_cells,functor);
319  }
320  }
321  else {
322  int spaceDim = basisValues.basis_vector.extent(3);
323  if(spaceDim==3) {
325  Kokkos::parallel_for(workset.num_cells,functor);
326  }
327  else {
329  Kokkos::parallel_for(workset.num_cells,functor);
330  }
331  }
332  }
333  else {
334  if(accelerate_jacobian) {
336  Kokkos::parallel_for(workset.num_cells,functor);
337  }
338  else {
340  Kokkos::parallel_for(workset.num_cells,functor);
341  }
342  }
343 }
344 
345 }
346 
347 #endif
PHX::MDField< ScalarT, Cell, Point > dof_ip_scalar
Definition: Panzer_DOF.hpp:89
DOF(const Teuchos::ParameterList &p)
T & get(const std::string &name, T def_value)
bool use_descriptors_
Definition: Panzer_DOF.hpp:83
Array_CellBasisIP basis_scalar
panzer::BasisDescriptor bd_
Definition: Panzer_DOF.hpp:84
Array_CellBasisIPDim basis_vector
Teuchos::RCP< PHX::DataLayout > dl_scalar
Data layout for scalar fields.
PHX::MDField< const ScalarT, Cell, Point > dof_basis
Definition: Panzer_DOF.hpp:87
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.
panzer::IntegrationDescriptor id_
Definition: Panzer_DOF.hpp:85
bool is_vector_basis
Definition: Panzer_DOF.hpp:95
Teuchos::RCP< PHX::DataLayout > dl_vector
Data layout for vector fields.
void evaluateFields(typename TRAITS::EvalData d)
Interpolates basis DOF values to IP DOF values.
Definition: Panzer_DOF.hpp:56
PHX::MDField< ScalarT, Cell, Point, Dim > dof_ip_vector
Definition: Panzer_DOF.hpp:90
#define TEUCHOS_ASSERT(assertion_test)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
std::string basis_name
Definition: Panzer_DOF.hpp:92
std::size_t basis_index
Definition: Panzer_DOF.hpp:93
const std::string & getType() const
Get type of basis.
Kokkos::View< const int *, PHX::Device > offsets