Panzer  Version of the Day
Panzer_DOFDiv_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_DIV_IMPL_HPP
44 #define PANZER_DOF_DIV_IMPL_HPP
45 
47 #include "Panzer_BasisIRLayout.hpp"
49 #include "Intrepid2_FunctionSpaceTools.hpp"
50 
51 namespace panzer {
52 
53 namespace {
54 
55 //**********************************************************************
56 template<typename ScalarT,typename ArrayT>
57 void evaluateDiv_withSens(int numCells,
58  int basis_dimension,
59  PHX::MDField<ScalarT,Cell,IP> & dof_div,
60  PHX::MDField<ScalarT,Cell,Point> & dof_value,
61  const ArrayT & div_basis)
62 {
63  if(numCells>0) {
64  // evaluate at quadrature points
65 
66  int numFields = div_basis.dimension(1);
67  int numPoints = div_basis.dimension(2);
68 
69  for (int cell=0; cell<numCells; cell++) {
70  for (int pt=0; pt<numPoints; pt++) {
71  // first initialize to the right thing (prevents over writing with 0)
72  // then loop over one less basis function
73  // ScalarT & div = dof_div(cell,pt);
74  dof_div(cell,pt) = dof_value(cell, 0) * div_basis(cell, 0, pt);
75  for (int bf=1; bf<numFields; bf++)
76  dof_div(cell,pt) += dof_value(cell, bf) * div_basis(cell, bf, pt);
77  }
78  }
79  }
80 }
81 
82 }
83 
84 //**********************************************************************
85 // MOST EVALUATION TYPES
86 //**********************************************************************
87 
88 //**********************************************************************
89 template<typename EvalT, typename TRAITS>
92  dof_value( p.get<std::string>("Name"),
93  p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->functional),
94  basis_name(p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->name())
95 {
97  = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
98 
99  // Verify that this basis supports the div operation
100  TEUCHOS_TEST_FOR_EXCEPTION(!basis->supportsDiv(),std::logic_error,
101  "DOFDiv: Basis of type \"" << basis->name() << "\" does not support DIV");
103  "DOFDiv: Basis of type \"" << basis->name() << "\" in DOF Div should require orientations. So we are throwing.");
104 
105  // build dof_div
107  dof_div = PHX::MDField<ScalarT,Cell,IP>(p.get<std::string>("Div Name"),
109 
110  // add to evaluation graph
111  this->addEvaluatedField(dof_div);
112  this->addDependentField(dof_value);
113 
114  std::string n = "DOFDiv: " + dof_div.fieldTag().name() + " ("+PHX::typeAsString<EvalT>()+")";
115  this->setName(n);
116 }
117 
118 //**********************************************************************
119 template<typename EvalT, typename TRAITS>
121 postRegistrationSetup(typename TRAITS::SetupData sd,
123 {
124  this->utils.setFieldData(dof_value,fm);
125  this->utils.setFieldData(dof_div,fm);
126 
127  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
128 }
129 
130 //**********************************************************************
131 template<typename EvalT, typename TRAITS>
133 evaluateFields(typename TRAITS::EvalData workset)
134 {
135  evaluateDiv_withSens<ScalarT>(workset.num_cells,basis_dimension,dof_div,dof_value,this->wda(workset).bases[basis_index]->div_basis);
136 }
137 
138 //**********************************************************************
139 
140 //**********************************************************************
141 // JACOBIAN EVALUATION TYPES
142 //**********************************************************************
143 
144 //**********************************************************************
145 template<typename TRAITS>
148  dof_value( p.get<std::string>("Name"),
149  p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->functional),
150  basis_name(p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->name())
151 {
153  = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
154 
155  // do you specialize because you know where the basis functions are and can
156  // skip a large number of AD calculations?
157  if(p.isType<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector")) {
158  offsets = *p.get<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector");
159  accelerate_jacobian = true; // short cut for identity matrix
160  }
161  else
162  accelerate_jacobian = false; // don't short cut for identity matrix
163  accelerate_jacobian = false; // don't short cut for identity matrix
164 
165  // Verify that this basis supports the div operation
166  TEUCHOS_TEST_FOR_EXCEPTION(!basis->supportsDiv(),std::logic_error,
167  "DOFDiv: Basis of type \"" << basis->name() << "\" does not support DIV");
169  "DOFDiv: Basis of type \"" << basis->name() << "\" in DOF Div should require orientations. So we are throwing.");
170 
171  // build dof_div
173  dof_div = PHX::MDField<ScalarT,Cell,IP>(p.get<std::string>("Div Name"),
175 
176  // add to evaluation graph
177  this->addEvaluatedField(dof_div);
178  this->addDependentField(dof_value);
179 
180  std::string n = "DOFDiv: " + dof_div.fieldTag().name() + " ("+PHX::typeAsString<panzer::Traits::Jacobian>()+")";
181  this->setName(n);
182 }
183 
184 //**********************************************************************
185 template<typename TRAITS>
187 postRegistrationSetup(typename TRAITS::SetupData sd,
189 {
190  this->utils.setFieldData(dof_value,fm);
191  this->utils.setFieldData(dof_div,fm);
192 
193  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
194 }
195 
196 template<typename TRAITS>
198 evaluateFields(typename TRAITS::EvalData workset)
199 {
200  if(!accelerate_jacobian) {
201  // do the case where we use the AD types to determine the derivatives
202  evaluateDiv_withSens(workset.num_cells,basis_dimension,dof_div,dof_value,this->wda(workset).bases[basis_index]->div_basis);
203  return;
204  }
205 
206  TEUCHOS_ASSERT(false);
207 }
208 
209 }
210 
211 #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::size_t basis_index
std::string name() const
A unique key that is the combination of the basis type and basis order.
std::string basis_name
DOFDiv(const Teuchos::ParameterList &p)
T & get(const std::string &name, T def_value)
bool supportsDiv() const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
int numPoints
int numFields
bool requiresOrientations() const
Teuchos::RCP< PHX::DataLayout > dl_scalar
Data layout for scalar fields.
PHX::MDField< ScalarT, Cell, Point > dof_value
bool isType(const std::string &name) const
int dimension() const
Returns the dimension of the basis from the topology.
PHX::MDField< ScalarT, Cell, Point > dof_value
Interpolates basis DOF values to IP DOF Gradient values.
void evaluateFields(typename TRAITS::EvalData d)
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
#define TEUCHOS_ASSERT(assertion_test)
std::string basis_name
Kokkos::View< const int *, PHX::Device > offsets
PHX::MDField< ScalarT, Cell, IP > dof_div