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  this->addDependentField(basisValues->basis_ref_scalar);
99  this->addDependentField(basisValues->basis_scalar);
100  }
101  else if(basis->isVectorBasis()) {
102  dof_ip_vector = PHX::MDField<ScalarT,Cell,Point,Dim>(
103  evalName,
104  pointRule->dl_vector);
105  this->addEvaluatedField(dof_ip_vector);
106  this->addDependentField(basisValues->basis_ref_vector);
107  this->addDependentField(basisValues->basis_vector);
108  }
109  else
110  { TEUCHOS_ASSERT(false); }
111 
112  std::string n = "DOF_PointValues: " + dof_basis.fieldTag().name();
113  this->setName(n);
114 }
115 
116 //**********************************************************************
117 template<typename EvalT, typename TRAITS>
119 postRegistrationSetup(typename TRAITS::SetupData /* sd */,
121 {
122  if(!is_vector_basis) {
123  this->utils.setFieldData(basisValues->basis_ref_scalar,fm);
124  this->utils.setFieldData(basisValues->basis_scalar,fm);
125  }
126  else {
127  this->utils.setFieldData(basisValues->basis_ref_vector,fm);
128  this->utils.setFieldData(basisValues->basis_vector,fm);
129  }
130 }
131 
132 //**********************************************************************
133 template<typename EvalT, typename TRAITS>
135 evaluateFields(typename TRAITS::EvalData workset)
136 {
137  const int vector_size = panzer::HP::inst().vectorSize<ScalarT>();
138 
139  if(is_vector_basis) {
140  int spaceDim = basisValues->basis_vector.extent(3);
141  if(spaceDim==3) {
142  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);
143  Kokkos::parallel_for(Kokkos::TeamPolicy<PHX::Device>(workset.num_cells,Kokkos::AUTO(),vector_size),functor);
144  }
145  else {
146  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);
147  Kokkos::parallel_for(Kokkos::TeamPolicy<PHX::Device>(workset.num_cells,Kokkos::AUTO(),vector_size),functor);
148  }
149  }
150  else {
151  dof_functors::EvaluateDOFWithSens_Scalar<ScalarT,typename BasisValues2<double>::Array_CellBasisIP> functor(dof_basis,dof_ip_scalar,basisValues->basis_scalar);
152  Kokkos::parallel_for(workset.num_cells,functor);
153  }
154 }
155 
156 //**********************************************************************
157 
158 //**********************************************************************
159 // JACOBIAN EVALUATION TYPES
160 //**********************************************************************
161 
162 //**********************************************************************
163 template<typename TRAITS>
166 {
167  const std::string fieldName = p.get<std::string>("Name");
168  basis = p.get< Teuchos::RCP<const PureBasis> >("Basis");
169  Teuchos::RCP<const PointRule> pointRule = p.get< Teuchos::RCP<const PointRule> >("Point Rule");
170  is_vector_basis = basis->isVectorBasis();
171 
172  if(p.isType<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector")) {
173  const std::vector<int> & offsets = *p.get<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector");
174 
175  // allocate and copy offsets vector to Kokkos array
176  offsets_array = Kokkos::View<int*,PHX::Device>("offsets",offsets.size());
177  for(std::size_t i=0;i<offsets.size();i++)
178  offsets_array(i) = offsets[i];
179 
180  accelerate_jacobian = true; // short cut for identity matrix
181  }
182  else
183  accelerate_jacobian = false; // don't short cut for identity matrix
184 
185  std::string evalName = fieldName+"_"+pointRule->getName();
186  if(p.isType<bool>("Use DOF Name")) {
187  if(p.get<bool>("Use DOF Name"))
188  evalName = fieldName;
189  }
190 
191  dof_basis = PHX::MDField<const ScalarT,Cell,Point>(fieldName, basis->functional);
192 
193  this->addDependentField(dof_basis);
194 
195  // setup all basis fields that are required
196  Teuchos::RCP<BasisIRLayout> layout = Teuchos::rcp(new BasisIRLayout(basis,*pointRule));
197  basisValues = Teuchos::rcp(new BasisValues2<double>(basis->name()+"_"+pointRule->getName()+"_"));
198  basisValues->setupArrays(layout,false);
199 
200  // the field manager will allocate all of these field
201  // swap between scalar basis value, or vector basis value
202  if(basis->isScalarBasis()) {
203  dof_ip_scalar = PHX::MDField<ScalarT,Cell,Point>(
204  evalName,
205  pointRule->dl_scalar);
206  this->addEvaluatedField(dof_ip_scalar);
207  this->addDependentField(basisValues->basis_ref_scalar);
208  this->addDependentField(basisValues->basis_scalar);
209  }
210  else if(basis->isVectorBasis()) {
211  dof_ip_vector = PHX::MDField<ScalarT,Cell,Point,Dim>(
212  evalName,
213  pointRule->dl_vector);
214  this->addEvaluatedField(dof_ip_vector);
215  this->addDependentField(basisValues->basis_ref_vector);
216  this->addDependentField(basisValues->basis_vector);
217  }
218  else
219  { TEUCHOS_ASSERT(false); }
220 
221  std::string n = "DOF_PointValues: " + dof_basis.fieldTag().name() + " Jacobian";
222  this->setName(n);
223 }
224 
225 //**********************************************************************
226 template<typename TRAITS>
228 postRegistrationSetup(typename TRAITS::SetupData /* sd */,
230 {
231  if(!is_vector_basis) {
232  this->utils.setFieldData(basisValues->basis_ref_scalar,fm);
233  this->utils.setFieldData(basisValues->basis_scalar,fm);
234  }
235  else {
236  this->utils.setFieldData(basisValues->basis_ref_vector,fm);
237  this->utils.setFieldData(basisValues->basis_vector,fm);
238  }
239 }
240 
241 //**********************************************************************
242 template<typename TRAITS>
244 evaluateFields(typename TRAITS::EvalData workset)
245 {
246  const int vector_size = panzer::HP::inst().vectorSize<ScalarT>();
247 
248  if(is_vector_basis) {
249  if(accelerate_jacobian) {
250  int spaceDim = basisValues->basis_vector.extent(3);
251  if(spaceDim==3) {
253  Kokkos::parallel_for(workset.num_cells,functor);
254  }
255  else {
257  Kokkos::parallel_for(workset.num_cells,functor);
258  }
259  }
260  else {
261  int spaceDim = basisValues->basis_vector.extent(3);
262  if(spaceDim==3) {
264  Kokkos::parallel_for(Kokkos::TeamPolicy<PHX::Device>(workset.num_cells,Kokkos::AUTO(),vector_size),functor);
265  }
266  else {
268  Kokkos::parallel_for(Kokkos::TeamPolicy<PHX::Device>(workset.num_cells,Kokkos::AUTO(),vector_size),functor);
269  }
270  }
271  }
272  else {
273  if(accelerate_jacobian) {
275  Kokkos::parallel_for(workset.num_cells,functor);
276  }
277  else {
279  Kokkos::parallel_for(workset.num_cells,functor);
280  }
281  }
282 }
283 
284 }
285 
286 #endif
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 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)
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)
Kokkos::View< const int *, PHX::Device > offsets