Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_DOF_PointValues_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_PointValues evaluator
60 //**********************************************************************
61 
62 //**********************************************************************
63 // MOST EVALUATION TYPES
64 //**********************************************************************
65 
66 //**********************************************************************
67 template<typename EvalT, typename TRAITS>
70 {
71  const std::string fieldName = p.get<std::string>("Name");
72  basis = p.get< Teuchos::RCP<const PureBasis> >("Basis");
74  is_vector_basis = basis->isVectorBasis();
75 
76  std::string evalName = fieldName+"_"+pointRule->getName();
77  if(p.isType<bool>("Use DOF Name")) {
78  if(p.get<bool>("Use DOF Name"))
79  evalName = fieldName;
80  }
81 
82  dof_basis = PHX::MDField<const ScalarT,Cell,Point>(fieldName, basis->functional);
83 
84  this->addDependentField(dof_basis);
85 
86  // setup all basis fields that are required
87  Teuchos::RCP<BasisIRLayout> layout = Teuchos::rcp(new BasisIRLayout(basis,*pointRule));
88  basisValues = Teuchos::rcp(new BasisValues2<double>(basis->name()+"_"+pointRule->getName()+"_"));
89  basisValues->setupArrays(layout,false);
90 
91  // the field manager will allocate all of these field
92  // swap between scalar basis value, or vector basis value
93  if(basis->isScalarBasis()) {
94  dof_ip_scalar = PHX::MDField<ScalarT,Cell,Point>(
95  evalName,
96  pointRule->dl_scalar);
97  this->addEvaluatedField(dof_ip_scalar);
98  constBasisRefScalar_ = basisValues->basis_ref_scalar;
99  constBasisScalar_ = basisValues->basis_scalar;
100  this->addDependentField(constBasisRefScalar_);
101  this->addDependentField(constBasisScalar_);
102  }
103  else if(basis->isVectorBasis()) {
104  dof_ip_vector = PHX::MDField<ScalarT,Cell,Point,Dim>(
105  evalName,
106  pointRule->dl_vector);
107  this->addEvaluatedField(dof_ip_vector);
108  constBasisRefVector_ = basisValues->basis_ref_vector;
109  constBasisVector_ = basisValues->basis_vector;
110  this->addDependentField(constBasisRefVector_);
111  this->addDependentField(constBasisVector_);
112  }
113  else
114  { TEUCHOS_ASSERT(false); }
115 
116  std::string n = "DOF_PointValues: " + dof_basis.fieldTag().name();
117  this->setName(n);
118 }
119 
120 //**********************************************************************
121 template<typename EvalT, typename TRAITS>
123 postRegistrationSetup(typename TRAITS::SetupData /* sd */,
125 {
126  if(!is_vector_basis) {
127  this->utils.setFieldData(basisValues->basis_ref_scalar,fm);
128  this->utils.setFieldData(basisValues->basis_scalar,fm);
129  }
130  else {
131  this->utils.setFieldData(basisValues->basis_ref_vector,fm);
132  this->utils.setFieldData(basisValues->basis_vector,fm);
133  }
134 }
135 
136 //**********************************************************************
137 template<typename EvalT, typename TRAITS>
139 evaluateFields(typename TRAITS::EvalData workset)
140 {
141  const int vector_size = panzer::HP::inst().vectorSize<ScalarT>();
142 
143  if(is_vector_basis) {
144  int spaceDim = basisValues->basis_vector.extent(3);
145  if(spaceDim==3) {
146  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);
147  Kokkos::parallel_for(Kokkos::TeamPolicy<PHX::Device>(workset.num_cells,Kokkos::AUTO(),vector_size),functor);
148  }
149  else {
150  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);
151  Kokkos::parallel_for(Kokkos::TeamPolicy<PHX::Device>(workset.num_cells,Kokkos::AUTO(),vector_size),functor);
152  }
153  }
154  else {
155  dof_functors::EvaluateDOFWithSens_Scalar<ScalarT,typename BasisValues2<double>::Array_CellBasisIP> functor(dof_basis,dof_ip_scalar,basisValues->basis_scalar);
156  Kokkos::parallel_for(workset.num_cells,functor);
157  }
158 }
159 
160 //**********************************************************************
161 
162 //**********************************************************************
163 // JACOBIAN EVALUATION TYPES
164 //**********************************************************************
165 
166 //**********************************************************************
167 template<typename TRAITS>
170 {
171  const std::string fieldName = p.get<std::string>("Name");
172  basis = p.get< Teuchos::RCP<const PureBasis> >("Basis");
173  Teuchos::RCP<const PointRule> pointRule = p.get< Teuchos::RCP<const PointRule> >("Point Rule");
174  is_vector_basis = basis->isVectorBasis();
175 
176  if(p.isType<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector")) {
177  const std::vector<int> & offsets = *p.get<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector");
178 
179  // allocate and copy offsets vector to Kokkos array
180  offsets_array = Kokkos::View<int*,PHX::Device>("offsets",offsets.size());
181  for(std::size_t i=0;i<offsets.size();i++)
182  offsets_array(i) = offsets[i];
183 
184  accelerate_jacobian = true; // short cut for identity matrix
185  }
186  else
187  accelerate_jacobian = false; // don't short cut for identity matrix
188 
189  std::string evalName = fieldName+"_"+pointRule->getName();
190  if(p.isType<bool>("Use DOF Name")) {
191  if(p.get<bool>("Use DOF Name"))
192  evalName = fieldName;
193  }
194 
195  dof_basis = PHX::MDField<const ScalarT,Cell,Point>(fieldName, basis->functional);
196 
197  this->addDependentField(dof_basis);
198 
199  // setup all basis fields that are required
200  Teuchos::RCP<BasisIRLayout> layout = Teuchos::rcp(new BasisIRLayout(basis,*pointRule));
201  basisValues = Teuchos::rcp(new BasisValues2<double>(basis->name()+"_"+pointRule->getName()+"_"));
202  basisValues->setupArrays(layout,false);
203 
204  // the field manager will allocate all of these field
205  // swap between scalar basis value, or vector basis value
206  if(basis->isScalarBasis()) {
207  dof_ip_scalar = PHX::MDField<ScalarT,Cell,Point>(
208  evalName,
209  pointRule->dl_scalar);
210  this->addEvaluatedField(dof_ip_scalar);
211  constBasisRefScalar_ = basisValues->basis_ref_scalar;
212  constBasisScalar_ = basisValues->basis_scalar;
213  this->addDependentField(constBasisRefScalar_);
214  this->addDependentField(constBasisScalar_);
215  }
216  else if(basis->isVectorBasis()) {
217  dof_ip_vector = PHX::MDField<ScalarT,Cell,Point,Dim>(
218  evalName,
219  pointRule->dl_vector);
220  this->addEvaluatedField(dof_ip_vector);
221  constBasisRefVector_ = basisValues->basis_ref_vector;
222  constBasisVector_ = basisValues->basis_vector;
223  this->addDependentField(constBasisRefVector_);
224  this->addDependentField(constBasisVector_);
225  }
226  else
227  { TEUCHOS_ASSERT(false); }
228 
229  std::string n = "DOF_PointValues: " + dof_basis.fieldTag().name() + " Jacobian";
230  this->setName(n);
231 }
232 
233 //**********************************************************************
234 template<typename TRAITS>
236 postRegistrationSetup(typename TRAITS::SetupData /* sd */,
238 {
239  if(!is_vector_basis) {
240  this->utils.setFieldData(basisValues->basis_ref_scalar,fm);
241  this->utils.setFieldData(basisValues->basis_scalar,fm);
242  }
243  else {
244  this->utils.setFieldData(basisValues->basis_ref_vector,fm);
245  this->utils.setFieldData(basisValues->basis_vector,fm);
246  }
247 }
248 
249 //**********************************************************************
250 template<typename TRAITS>
252 evaluateFields(typename TRAITS::EvalData workset)
253 {
254  const int vector_size = panzer::HP::inst().vectorSize<ScalarT>();
255 
256  if(is_vector_basis) {
257  if(accelerate_jacobian) {
258  int spaceDim = basisValues->basis_vector.extent(3);
259  if(spaceDim==3) {
261  Kokkos::parallel_for(workset.num_cells,functor);
262  }
263  else {
265  Kokkos::parallel_for(workset.num_cells,functor);
266  }
267  }
268  else {
269  int spaceDim = basisValues->basis_vector.extent(3);
270  if(spaceDim==3) {
272  Kokkos::parallel_for(Kokkos::TeamPolicy<PHX::Device>(workset.num_cells,Kokkos::AUTO(),vector_size),functor);
273  }
274  else {
276  Kokkos::parallel_for(Kokkos::TeamPolicy<PHX::Device>(workset.num_cells,Kokkos::AUTO(),vector_size),functor);
277  }
278  }
279  }
280  else {
281  if(accelerate_jacobian) {
283  Kokkos::parallel_for(workset.num_cells,functor);
284  }
285  else {
287  Kokkos::parallel_for(workset.num_cells,functor);
288  }
289  }
290 }
291 
292 }
293 
294 #endif
PHX::MDField< const double, BASIS, IP, void, void, void, void, void, void > constBasisRefScalar_
T & get(const std::string &name, T def_value)
void evaluateFields(typename TRAITS::EvalData d)
PHX::MDField< ScalarT, Cell, Point, Dim > dof_ip_vector
PHX::MDField< const double, Cell, BASIS, IP, void, void, void, void, void > constBasisScalar_
PHX::MDField< const ScalarT, Cell, Point > dof_basis
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< const PureBasis > basis
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
PHX::MDField< const double, Cell, BASIS, IP, Dim, void, void, void, void > constBasisVector_
int vectorSize() const
Returns the vector size. Specialized for AD scalar types.
bool isType(const std::string &name) const
static HP & inst()
Private ctor.
PHX::MDField< ScalarT, Cell, Point > dof_ip_scalar
Teuchos::RCP< BasisValues2< double > > basisValues
#define TEUCHOS_ASSERT(assertion_test)
DOF_PointValues(const Teuchos::ParameterList &p)
PHX::MDField< const double, BASIS, IP, Dim, void, void, void, void, void > constBasisRefVector_
Kokkos::View< const int *, PHX::Device > offsets