Panzer  Version of the Day
Panzer_Integrator_GradBasisDotVector_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_EVALUATOR_GRADBASISDOTVECTOR_IMPL_HPP
44 #define PANZER_EVALUATOR_GRADBASISDOTVECTOR_IMPL_HPP
45 
46 #include "Intrepid2_FunctionSpaceTools.hpp"
48 #include "Panzer_BasisIRLayout.hpp"
50 #include "Kokkos_ViewFactory.hpp"
51 
52 #define PANZER_USE_FAST_QUAD 1
53 // #define PANZER_USE_FAST_QUAD 0
54 
55 namespace panzer {
56 
57 //**********************************************************************
58 PHX_EVALUATOR_CTOR(Integrator_GradBasisDotVector,p) :
59  residual( p.get<std::string>("Residual Name"),
60  p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->functional),
61 // flux( p.get<std::string>("Flux Name"),
62 // p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_vector ),
63  basis_name(p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->name())
64 {
65  Teuchos::RCP<const PureBasis> basis
66  = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
67 
68 
69  // Default to this so we don't break anything
70  Teuchos::RCP<PHX::DataLayout> vector = p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_vector;
71  if(p.isType<Teuchos::RCP<PHX::DataLayout> >("Vector Data Layout")){
72  vector = p.get<Teuchos::RCP<PHX::DataLayout> >("Vector Data Layout");
73  }
74 
75  flux = PHX::MDField<const ScalarT,Cell,IP,Dim>( p.get<std::string>("Flux Name"), vector);
76 
77  // Number of dimensions has to be based on the spatial dimensions NOT THE DIMENSIONS OF THE VECTOR
78  num_dim = p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_vector->dimension(2);
79 
80  // Verify that this basis supports the gradient operation
81  TEUCHOS_TEST_FOR_EXCEPTION(!basis->supportsGrad(),std::logic_error,
82  "Integrator_GradBasisDotVector: Basis of type \"" << basis->name() << "\" does not support GRAD");
83 
84  // Make sure the Grad dimensions includes the vector
85  TEUCHOS_TEST_FOR_EXCEPTION(vector->dimension(2) < num_dim,std::logic_error,
86  "Integrator_GradBasisDotVector: Dimension of space exceeds dimension of vector.");
87 
88  this->addEvaluatedField(residual);
89  this->addDependentField(flux);
90 
91  multiplier = p.get<double>("Multiplier");
92 
93  if (p.isType<Teuchos::RCP<const std::vector<std::string> > >("Field Multipliers"))
94  {
95  const std::vector<std::string>& field_multiplier_names =
96  *(p.get<Teuchos::RCP<const std::vector<std::string> > >("Field Multipliers"));
97 
98  for (std::vector<std::string>::const_iterator name = field_multiplier_names.begin();
99  name != field_multiplier_names.end(); ++name)
100  {
101  PHX::MDField<const ScalarT,Cell,IP> tmp_field(*name, p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_scalar);
102  field_multipliers.push_back(tmp_field);
103  }
104  }
105 
106  for (auto & field : field_multipliers){
107  this->addDependentField(field);
108  }
109 
110  std::string n =
111  "Integrator_GradBasisDotVector: " + residual.fieldTag().name();
112 
113  this->setName(n);
114 }
115 
116 //**********************************************************************
117 PHX_POST_REGISTRATION_SETUP(Integrator_GradBasisDotVector,sd,fm)
118 {
119  this->utils.setFieldData(residual,fm);
120  this->utils.setFieldData(flux,fm);
121 
122  for (auto & field : field_multipliers)
123  this->utils.setFieldData(field,fm);
124 
125  num_nodes = residual.dimension(1);
126  num_qp = flux.dimension(1);
127 
128  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
129 
130  tmp = Kokkos::createDynRankView(residual.get_static_view(),"tmp",flux.dimension(0), num_qp, num_dim);
131 }
132 
133 //**********************************************************************
134 PHX_EVALUATE_FIELDS(Integrator_GradBasisDotVector,workset)
135 {
136  //for (int i=0; i < residual.size(); ++i)
137  // residual[i] = 0.0;
138 
139  Kokkos::deep_copy(residual.get_static_view(), ScalarT(0.0));
140 
141 #if PANZER_USE_FAST_QUAD
142  // do a scaled copy
143  for (int i=0; i < flux.extent_int(0); ++i)
144  for (int j=0; j < flux.extent_int(1); ++j)
145  for (int k=0; k < num_dim; ++k)
146  tmp(i,j,k) = multiplier * flux(i,j,k);
147 //Irina modified
148 // for (int i=0; i < flux.size(); ++i)
149 // tmp[i] = multiplier * flux[i];
150 
151  for (auto & field : field_multipliers) {
152  PHX::MDField<const ScalarT,Cell,IP> field_data = field;
153 
154  for (index_t cell = 0; cell < workset.num_cells; ++cell) {
155  for (std::size_t qp = 0; qp < num_qp; ++qp) {
156  ScalarT tmpVar = field_data(cell,qp);
157 
158  for (std::size_t dim = 0; dim < num_dim; ++dim)
159  tmp(cell,qp,dim) *= tmpVar;
160  }
161  }
162  }
163 
164  // const Kokkos::DynRankView<double,PHX::Device> & weighted_grad_basis = this->wda(workset).bases[basis_index]->weighted_grad_basis;
165  const BasisValues2<double> & bv = *this->wda(workset).bases[basis_index];
166 
167  // perform integration and vector dot product (at the same time! whoah!)
168  for (index_t cell = 0; cell < workset.num_cells; ++cell) {
169  for (std::size_t basis = 0; basis < num_nodes; ++basis) {
170  for (std::size_t qp = 0; qp < num_qp; ++qp) {
171  for (std::size_t dim = 0; dim < num_dim; ++dim) {
172  residual(cell,basis) += tmp(cell,qp,dim)*bv.weighted_grad_basis(cell,basis,qp,dim);
173 
174  /*
175  std::cout << residual(cell,basis) << " "
176  << tmp(cell,qp,dim) << " "
177  << bv.weighted_grad_basis(cell,basis,qp,dim) << std::endl;
178  */
179  }
180  }
181  }
182  }
183 #else
184  for (index_t cell = 0; cell < workset.num_cells; ++cell)
185  {
186  for (std::size_t qp = 0; qp < num_qp; ++qp)
187  {
188  ScalarT tmpVar = 1.0;
189  for (typename std::vector<PHX::MDField<const ScalarT,Cell,IP> >::iterator field = field_multipliers.begin();
190  field != field_multipliers.end(); ++field)
191  tmpVar = tmpVar * (*field)(cell,qp);
192 
193  for (std::size_t dim = 0; dim < num_dim; ++dim)
194  tmp(cell,qp,dim) = multiplier * tmpVar * flux(cell,qp,dim);
195  }
196  }
197 
198  if(workset.num_cells>0)
199  Intrepid2::FunctionSpaceTools::
200  integrate<ScalarT>(residual, tmp,
201  (this->wda(workset).bases[basis_index])->weighted_grad_basis,
202  Intrepid2::COMP_CPP);
203 #endif
204 }
205 
206 //**********************************************************************
207 
208 }
209 
210 #endif
PHX::MDField< ScalarT > residual
Evaluates a Dirichlet BC residual corresponding to a field value.
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 > tmp
std::size_t num_qp
PHX::MDField< ScalarT > flux
PHX::MDField< ScalarT > vector
Kokkos::DynRankView< typename InputArray::value_type, PHX::Device > createDynRankView(const InputArray &a, const std::string &name, const DimensionPack... dims)
Wrapper to simplify Panzer use of Sacado ViewFactory.
PHX_EVALUATOR_CTOR(BasisValues_Evaluator, p)
PHX::MDField< const ScalarT, Cell, IP > field
double multiplier
Array_CellBasisIPDim weighted_grad_basis
PHX_EVALUATE_FIELDS(BasisValues_Evaluator, workset)
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
std::vector< PHX::MDField< const ScalarT, Cell, IP > > field_multipliers
PHX_POST_REGISTRATION_SETUP(BasisValues_Evaluator, sd, fm)