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>
68 DOF(const Teuchos::ParameterList & p) :
69  dof_basis( p.get<std::string>("Name"),
70  p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->functional),
71  basis_name(p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->name())
72 {
73  Teuchos::RCP<const PureBasis> basis
74  = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
75  is_vector_basis = basis->isVectorBasis();
76 
77  // swap between scalar basis value, or vector basis value
78  if(basis->isScalarBasis()) {
79  dof_ip_scalar = PHX::MDField<ScalarT,Cell,Point>(
80  p.get<std::string>("Name"),
81  p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_scalar);
82  this->addEvaluatedField(dof_ip_scalar);
83  }
84  else if(basis->isVectorBasis()) {
85  dof_ip_vector = PHX::MDField<ScalarT,Cell,Point,Dim>(
86  p.get<std::string>("Name"),
87  p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_vector);
88  this->addEvaluatedField(dof_ip_vector);
89  }
90  else
91  { TEUCHOS_ASSERT(false); }
92 
93  this->addDependentField(dof_basis);
94 
95  std::string n = "DOF: " + dof_basis.fieldTag().name() + " ("+PHX::typeAsString<EvalT>()+")";
96  this->setName(n);
97 }
98 
99 //**********************************************************************
100 template<typename EvalT, typename TRAITS>
102 postRegistrationSetup(typename TRAITS::SetupData sd,
104 {
105  this->utils.setFieldData(dof_basis,fm);
106  if(is_vector_basis)
107  this->utils.setFieldData(dof_ip_vector,fm);
108  else
109  this->utils.setFieldData(dof_ip_scalar,fm);
110 
111  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
112 }
113 
114 //**********************************************************************
115 template<typename EvalT, typename TRAITS>
117 evaluateFields(typename TRAITS::EvalData workset)
118 {
119  panzer::BasisValues2<double> & basisValues = *this->wda(workset).bases[basis_index];
120 
121  if(is_vector_basis) {
122  int spaceDim = basisValues.basis_vector.dimension(3);
123  if(spaceDim==3) {
125  Kokkos::parallel_for(workset.num_cells,functor);
126  }
127  else {
129  Kokkos::parallel_for(workset.num_cells,functor);
130  }
131  }
132  else {
134  Kokkos::parallel_for(workset.num_cells,functor);
135  }
136 }
137 
138 //**********************************************************************
139 
140 //**********************************************************************
141 // JACOBIAN EVALUATION TYPES
142 //**********************************************************************
143 
144 //**********************************************************************
145 template<typename TRAITS>
147 DOF(const Teuchos::ParameterList & p) :
148  dof_basis( p.get<std::string>("Name"),
149  p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->functional),
150  basis_name(p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->name())
151 {
152  Teuchos::RCP<const PureBasis> basis
153  = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
154  is_vector_basis = basis->isVectorBasis();
155 
156  if(p.isType<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector")) {
157  const std::vector<int> & offsets = *p.get<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector");
158 
159  // allocate and copy offsets vector to Kokkos array
160  offsets_array = Kokkos::View<int*,PHX::Device>("offsets",offsets.size());
161  for(std::size_t i=0;i<offsets.size();i++)
162  offsets_array(i) = offsets[i];
163 
164  accelerate_jacobian_enabled = true; // short cut for identity matrix
165 
166  // get the sensitivities name that is valid for accelerated jacobians
167  sensitivities_name = true;
168  if (p.isType<std::string>("Sensitivities Name"))
169  sensitivities_name = p.get<std::string>("Sensitivities Name");
170  }
171  else
172  accelerate_jacobian_enabled = false; // don't short cut for identity matrix
173 
174  // swap between scalar basis value, or vector basis value
175  if(basis->isScalarBasis()) {
176  dof_ip_scalar = PHX::MDField<ScalarT,Cell,Point>(
177  p.get<std::string>("Name"),
178  p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_scalar);
179  this->addEvaluatedField(dof_ip_scalar);
180  }
181  else if(basis->isVectorBasis()) {
182  dof_ip_vector = PHX::MDField<ScalarT,Cell,Point,Dim>(
183  p.get<std::string>("Name"),
184  p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_vector);
185  this->addEvaluatedField(dof_ip_vector);
186  }
187  else
188  { TEUCHOS_ASSERT(false); }
189 
190  this->addDependentField(dof_basis);
191 
192  std::string n = "DOF: " + dof_basis.fieldTag().name()
193  + ( accelerate_jacobian_enabled ? " accel_jac " : "slow_jac" )
194  + " ("+PHX::typeAsString<panzer::Traits::Jacobian>()+")";
195  this->setName(n);
196 }
197 
198 //**********************************************************************
199 template<typename TRAITS>
201 postRegistrationSetup(typename TRAITS::SetupData sd,
203 {
204  this->utils.setFieldData(dof_basis,fm);
205  if(is_vector_basis)
206  this->utils.setFieldData(dof_ip_vector,fm);
207  else
208  this->utils.setFieldData(dof_ip_scalar,fm);
209 
210  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
211 }
212 
213 // **********************************************************************
214 template<typename TRAITS>
216 preEvaluate(typename TRAITS::PreEvalData d)
217 {
218  // if sensitivities were requrested for this field enable accelerated
219  // jacobian calculations
220  accelerate_jacobian = false;
221  if(accelerate_jacobian_enabled && d.first_sensitivities_name==sensitivities_name) {
222  accelerate_jacobian = true;
223  }
224 }
225 
226 //**********************************************************************
227 template<typename TRAITS>
229 evaluateFields(typename TRAITS::EvalData workset)
230 {
231  panzer::BasisValues2<double> & basisValues = *this->wda(workset).bases[basis_index];
232 
233  if(is_vector_basis) {
234  if(accelerate_jacobian) {
235  int spaceDim = basisValues.basis_vector.dimension(3);
236  if(spaceDim==3) {
238  Kokkos::parallel_for(workset.num_cells,functor);
239  }
240  else {
242  Kokkos::parallel_for(workset.num_cells,functor);
243  }
244  }
245  else {
246  int spaceDim = basisValues.basis_vector.dimension(3);
247  if(spaceDim==3) {
249  Kokkos::parallel_for(workset.num_cells,functor);
250  }
251  else {
253  Kokkos::parallel_for(workset.num_cells,functor);
254  }
255  }
256  }
257  else {
258  if(accelerate_jacobian) {
260  Kokkos::parallel_for(workset.num_cells,functor);
261  }
262  else {
264  Kokkos::parallel_for(workset.num_cells,functor);
265  }
266  }
267 }
268 
269 }
270 
271 #endif
std::size_t basis_index
std::vector< std::string >::size_type getBasisIndex(std::string basis_name, panzer::Workset &workset, WorksetDetailsAccessor &wda)
Returns the index in the workset bases for a particular BasisIRLayout name.
std::string basis_name
PHX::MDField< ScalarT, Cell, Point > dof_ip_scalar
DOF(const Teuchos::ParameterList &p)
Teuchos::RCP< BasisValues2< ScalarT > > basisValues
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
void evaluateFields(typename TRAITS::EvalData d)
Interpolates basis DOF values to IP DOF values.
PHX::MDField< ScalarT, Cell, Point, Dim > dof_ip_vector
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
std::string basis_name
std::size_t basis_index
PHX::MDField< ScalarT, Cell, Point > dof_basis
Kokkos::View< const int *, PHX::Device > offsets