Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
53 
54 #include "Intrepid2_FunctionSpaceTools.hpp"
55 
56 namespace panzer {
57 
58 //**********************************************************************
59 //* DOF evaluator
60 //**********************************************************************
61 
62 //**********************************************************************
63 // MOST EVALUATION TYPES
64 //**********************************************************************
65 
66 //**********************************************************************
67 template<typename EvalT, typename TRAITS>
70  use_descriptors_(false),
71  dof_basis( p.get<std::string>("Name"),
72  p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->functional),
73  basis_name(p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->name())
74 {
76  = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
77  is_vector_basis = basis->isVectorBasis();
78 
79  // swap between scalar basis value, or vector basis value
80  if(basis->isScalarBasis()) {
81  dof_ip_scalar = PHX::MDField<ScalarT,Cell,Point>(
82  p.get<std::string>("Name"),
83  p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_scalar);
84  this->addEvaluatedField(dof_ip_scalar);
85  }
86  else if(basis->isVectorBasis()) {
87  dof_ip_vector = PHX::MDField<ScalarT,Cell,Point,Dim>(
88  p.get<std::string>("Name"),
89  p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_vector);
90  this->addEvaluatedField(dof_ip_vector);
91  }
92  else
93  { TEUCHOS_ASSERT(false); }
94 
95  this->addDependentField(dof_basis);
96 
97  std::string n = "DOF: " + dof_basis.fieldTag().name() + " ("+PHX::typeAsString<EvalT>()+")";
98  this->setName(n);
99 }
100 
101 //**********************************************************************
102 template<typename EvalT, typename TRAITS>
104 DOF(const PHX::FieldTag & input,
105  const PHX::FieldTag & output,
106  const panzer::BasisDescriptor & bd,
108  : use_descriptors_(true)
109  , bd_(bd)
110  , id_(id)
111  , dof_basis(input)
112 {
113  TEUCHOS_ASSERT(bd.getType()=="HGrad" || bd.getType()=="HCurl" ||
114  bd.getType()=="HDiv" || bd.getType()=="Const")
115 
116  is_vector_basis = (bd.getType()=="HCurl" || bd.getType()=="HDiv");
117 
118  // swap between scalar basis value, or vector basis value
119  if(not is_vector_basis) {
120  dof_ip_scalar = output;
121  this->addEvaluatedField(dof_ip_scalar);
122  }
123  else {
124  dof_ip_vector = output;
125  this->addEvaluatedField(dof_ip_vector);
126  }
127 
128  this->addDependentField(dof_basis);
129 
130  std::string n = "DOF: " + dof_basis.fieldTag().name() + " ("+PHX::typeAsString<EvalT>()+")";
131  this->setName(n);
132 }
133 
134 //**********************************************************************
135 template<typename EvalT, typename TRAITS>
137 postRegistrationSetup(typename TRAITS::SetupData sd,
139 {
140  this->utils.setFieldData(dof_basis,fm);
141  if(is_vector_basis)
142  this->utils.setFieldData(dof_ip_vector,fm);
143  else
144  this->utils.setFieldData(dof_ip_scalar,fm);
145 
146  // descriptors don't access the basis values in the same way
147  if(not use_descriptors_)
148  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
149 }
150 
151 //**********************************************************************
152 template<typename EvalT, typename TRAITS>
154 evaluateFields(typename TRAITS::EvalData workset)
155 {
156  const panzer::BasisValues2<double> & basisValues = use_descriptors_ ? this->wda(workset).getBasisValues(bd_,id_)
157  : *this->wda(workset).bases[basis_index];
158 
159  const auto policy = panzer::HP::inst().teamPolicy<ScalarT,PHX::exec_space>(workset.num_cells);
160  const bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
161 
162  if(is_vector_basis) {
163  int spaceDim = basisValues.basis_vector.extent(3);
164  if(spaceDim==3) {
165  dof_functors::EvaluateDOFWithSens_Vector<ScalarT,typename BasisValues2<double>::Array_CellBasisIPDim,3> functor(dof_basis.get_static_view(),dof_ip_vector.get_static_view(),basisValues.basis_vector,use_shared_memory);
166  Kokkos::parallel_for(policy,functor,this->getName());
167  }
168  else {
169  dof_functors::EvaluateDOFWithSens_Vector<ScalarT,typename BasisValues2<double>::Array_CellBasisIPDim,2> functor(dof_basis.get_static_view(),dof_ip_vector.get_static_view(),basisValues.basis_vector,use_shared_memory);
170  Kokkos::parallel_for(policy,functor,this->getName());
171  }
172  }
173  else {
175  Kokkos::parallel_for(workset.num_cells,functor);
176  }
177 }
178 
179 //**********************************************************************
180 
181 //**********************************************************************
182 // JACOBIAN EVALUATION TYPES
183 //**********************************************************************
184 
185 //**********************************************************************
186 template<typename TRAITS>
189  use_descriptors_(false),
190  dof_basis( p.get<std::string>("Name"),
191  p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->functional),
192  basis_name(p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->name())
193 {
195  = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
196  is_vector_basis = basis->isVectorBasis();
197 
198  if(p.isType<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector")) {
199  const std::vector<int> & offsets = *p.get<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector");
200 
201  // allocate and copy offsets vector to Kokkos array
202  offsets_array = Kokkos::View<int*,PHX::Device>("offsets",offsets.size());
203  for(std::size_t i=0;i<offsets.size();i++)
204  offsets_array(i) = offsets[i];
205 
206  accelerate_jacobian_enabled = true; // short cut for identity matrix
207 
208  // get the sensitivities name that is valid for accelerated jacobians
209  sensitivities_name = true;
210  if (p.isType<std::string>("Sensitivities Name"))
211  sensitivities_name = p.get<std::string>("Sensitivities Name");
212  }
213  else
214  accelerate_jacobian_enabled = false; // don't short cut for identity matrix
215 
216  // swap between scalar basis value, or vector basis value
217  if(basis->isScalarBasis()) {
218  dof_ip_scalar = PHX::MDField<ScalarT,Cell,Point>(
219  p.get<std::string>("Name"),
220  p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_scalar);
221  this->addEvaluatedField(dof_ip_scalar);
222  }
223  else if(basis->isVectorBasis()) {
224  dof_ip_vector = PHX::MDField<ScalarT,Cell,Point,Dim>(
225  p.get<std::string>("Name"),
226  p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_vector);
227  this->addEvaluatedField(dof_ip_vector);
228  }
229  else
230  { TEUCHOS_ASSERT(false); }
231 
232  this->addDependentField(dof_basis);
233 
234  std::string n = "DOF: " + dof_basis.fieldTag().name()
235  + ( accelerate_jacobian_enabled ? " accel_jac " : "slow_jac" )
236  + " ("+PHX::typeAsString<panzer::Traits::Jacobian>()+")";
237  this->setName(n);
238 }
239 
240 //**********************************************************************
241 template<typename TRAITS>
243 DOF(const PHX::FieldTag & input,
244  const PHX::FieldTag & output,
245  const panzer::BasisDescriptor & bd,
247  : use_descriptors_(true)
248  , bd_(bd)
249  , id_(id)
250  , dof_basis(input)
251 {
252  TEUCHOS_ASSERT(bd.getType()=="HGrad" || bd.getType()=="HCurl" ||
253  bd.getType()=="HDiv" || bd.getType()=="Const")
254 
255  accelerate_jacobian_enabled = false; // don't short cut for identity matrix
256 
257  is_vector_basis = (bd.getType()=="HCurl" || bd.getType()=="HDiv");
258 
259  // swap between scalar basis value, or vector basis value
260  if(not is_vector_basis) {
261  dof_ip_scalar = output;
262  this->addEvaluatedField(dof_ip_scalar);
263  }
264  else {
265  dof_ip_vector = output;
266  this->addEvaluatedField(dof_ip_vector);
267  }
268 
269  this->addDependentField(dof_basis);
270 
271  std::string n = "DOF: " + dof_basis.fieldTag().name() + " slow_jac(descriptor) ("+PHX::typeAsString<typename TRAITS::Jacobian>()+")";
272  this->setName(n);
273 }
274 
275 //**********************************************************************
276 template<typename TRAITS>
278 postRegistrationSetup(typename TRAITS::SetupData sd,
280 {
281  this->utils.setFieldData(dof_basis,fm);
282  if(is_vector_basis)
283  this->utils.setFieldData(dof_ip_vector,fm);
284  else
285  this->utils.setFieldData(dof_ip_scalar,fm);
286 
287  // descriptors don't access the basis values in the same way
288  if(not use_descriptors_)
289  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
290 }
291 
292 // **********************************************************************
293 template<typename TRAITS>
295 preEvaluate(typename TRAITS::PreEvalData d)
296 {
297  // if sensitivities were requrested for this field enable accelerated
298  // jacobian calculations
299  accelerate_jacobian = false;
300  if(accelerate_jacobian_enabled && d.first_sensitivities_name==sensitivities_name) {
301  accelerate_jacobian = true;
302  }
303 }
304 
305 //**********************************************************************
306 template<typename TRAITS>
308 evaluateFields(typename TRAITS::EvalData workset)
309 {
310  const panzer::BasisValues2<double> & basisValues = use_descriptors_ ? this->wda(workset).getBasisValues(bd_,id_)
311  : *this->wda(workset).bases[basis_index];
312 
313  if(is_vector_basis) {
314  if(accelerate_jacobian) {
315  int spaceDim = basisValues.basis_vector.extent(3);
316  if(spaceDim==3) {
318  Kokkos::parallel_for(workset.num_cells,functor);
319  }
320  else {
322  Kokkos::parallel_for(workset.num_cells,functor);
323  }
324  }
325  else {
326  const bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
327  const auto policy = panzer::HP::inst().teamPolicy<ScalarT,PHX::exec_space>(workset.num_cells);
328  const int spaceDim = basisValues.basis_vector.extent(3);
329  if(spaceDim==3) {
331  Kokkos::parallel_for(policy,functor,this->getName());
332  }
333  else {
335  Kokkos::parallel_for(policy,functor,this->getName());
336  }
337  }
338  }
339  else {
340  if(accelerate_jacobian) {
342  Kokkos::parallel_for(workset.num_cells,functor);
343  }
344  else {
346  Kokkos::parallel_for(workset.num_cells,functor);
347  }
348  }
349 }
350 
351 }
352 
353 #endif
EvalT::ScalarT ScalarT
Definition: Panzer_DOF.hpp:81
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)
Kokkos::TeamPolicy< TeamPolicyProperties...> teamPolicy(const int &league_size)
Returns a TeamPolicy for hierarchic parallelism.
bool use_descriptors_
Definition: Panzer_DOF.hpp:83
Array_CellBasisIP basis_scalar
panzer::BasisDescriptor bd_
Definition: Panzer_DOF.hpp:84
Array_CellBasisIPDim basis_vector
const std::string & getType() const
Get type of basis.
panzer::Traits::Jacobian::ScalarT ScalarT
Definition: Panzer_DOF.hpp:123
PHX::MDField< const ScalarT, Cell, Point > dof_basis
Definition: Panzer_DOF.hpp:87
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
void evaluateFields(typename TRAITS::EvalData d)
bool isType(const std::string &name) const
static HP & inst()
Private ctor.
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
Kokkos::View< const int *, PHX::Device > offsets