Panzer  Version of the Day
Panzer_Workset.cpp
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 #include "Panzer_Workset.hpp"
44 
45 #include "Phalanx_DataLayout.hpp"
46 #include "Phalanx_DataLayout_MDALayout.hpp"
47 
49 #include "Panzer_Workset_Builder.hpp"
50 #include "Panzer_WorksetNeeds.hpp"
51 #include "Panzer_Dimension.hpp"
52 #include "Panzer_LocalMeshInfo.hpp"
54 #include "Panzer_PointValues2.hpp"
55 
57 
58 #include "Shards_CellTopology.hpp"
59 
60 namespace panzer {
61 
62 void
64  const panzer::WorksetNeeds & needs)
65 {
66 
67 
68  const size_t num_cells = partition.local_cells.extent(0);
69 
73 
74  subcell_index = -1;
75  block_id = partition.element_block_name;
76 
77  Kokkos::View<int*, PHX::Device> cell_ids = Kokkos::View<int*, PHX::Device>("cell_ids",num_cells);
78  Kokkos::deep_copy(cell_ids, partition.local_cells);
79  cell_local_ids_k = cell_ids;
80 
81  cell_local_ids.resize(num_cells,-1);
82  for(size_t cell=0;cell<num_cells;++cell){
83  const int local_cell = partition.local_cells(cell);
84  cell_local_ids[cell] = local_cell;
85  }
86 
87  auto fc = Teuchos::rcp(new panzer::FaceConnectivity());
88  fc->setup(partition);
89  _face_connectivity = fc;
90 
91  setupNeeds(partition.cell_topology, partition.cell_vertices, needs);
92 }
93 
94 void WorksetDetails::setupNeeds(Teuchos::RCP<const shards::CellTopology> cell_topology,
95  const Kokkos::View<double***,PHX::Device> & cell_vertices,
96  const panzer::WorksetNeeds & needs)
97 {
98 
99  const size_t num_cells = cell_vertices.extent(0);
100  const size_t num_vertices_per_cell = cell_vertices.extent(1);
101  const size_t num_dims_per_vertex = cell_vertices.extent(2);
102 
103  // Set cell vertices
104  {
105 
106  MDFieldArrayFactory af("",true);
107 
108  cell_vertex_coordinates = af.template buildStaticArray<double, Cell, NODE, Dim>("cell_vertices",num_cells, num_vertices_per_cell, num_dims_per_vertex);
109 
110  for(size_t i=0;i<num_cells;++i)
111  for(size_t j=0;j<num_vertices_per_cell;++j)
112  for(size_t k=0;k<num_dims_per_vertex;++k)
113  cell_vertex_coordinates(i,j,k) = cell_vertices(i,j,k);
114 
115  }
116 
117  // DEPRECATED - makes sure deprecated arrays are filled with something - this will probably segfault or throw an error
118  panzer::populateValueArrays(num_cells, false, needs, *this);
119 
120  const std::vector<panzer::BasisDescriptor> & basis_descriptors = needs.getBases();
121  const std::vector<panzer::IntegrationDescriptor> & integration_descriptors = needs.getIntegrators();
122  const std::vector<panzer::PointDescriptor> & point_descriptors = needs.getPoints();
123 
124  // Create the pure basis
125  for(const panzer::BasisDescriptor & basis_description : basis_descriptors){
126  // Create and store integration rule
127  Teuchos::RCP<panzer::PureBasis> basis = Teuchos::rcp(new panzer::PureBasis(basis_description, cell_topology, num_cells));
128  _pure_basis_map[basis_description.getKey()] = basis;
129  }
130 
131  // Create the integration terms and basis-integrator pairs
132  for(const panzer::IntegrationDescriptor & integration_description : integration_descriptors){
133 
134  int num_faces = -1;
135  if(integration_description.getType() == panzer::IntegrationRule::SURFACE){
136  num_faces = getFaceConnectivity().numSubcells();
137  }
138 
139  // Create and store integration rule
140  Teuchos::RCP<panzer::IntegrationRule> ir = Teuchos::rcp(new panzer::IntegrationRule(integration_description, cell_topology, num_cells, num_faces));
141  _integration_rule_map[integration_description.getKey()] = ir;
142 
143  // Create and store integration values
144  Teuchos::RCP<panzer::IntegrationValues2<double> > iv = Teuchos::rcp(new panzer::IntegrationValues2<double>("",true));
145  iv->setupArrays(ir);
146  iv->evaluateValues(cell_vertex_coordinates);
147  _integrator_map[integration_description.getKey()] = iv;
148 
149  // We need to generate a integration rule - basis pair for each basis
150  for(const panzer::BasisDescriptor & basis_description : basis_descriptors){
151 
152  // Grab the basis that was pre-calculated
153  const Teuchos::RCP<const panzer::PureBasis> & basis = _pure_basis_map[basis_description.getKey()];
154 
155  // Create a basis ir layout for this pair of integrator and basis
156  Teuchos::RCP<panzer::BasisIRLayout> b_layout = Teuchos::rcp(new panzer::BasisIRLayout(basis,*ir));
157 
158  // Create and store basis values
159  Teuchos::RCP<panzer::BasisValues2<double> > bv = Teuchos::rcp(new panzer::BasisValues2<double>("",true,true));
160  bv->setupArrays(b_layout);
161  if(ir->getType() == panzer::IntegrationDescriptor::SURFACE){
162  bv->evaluateValues(iv->ref_ip_coordinates,
163  iv->jac,
164  iv->jac_det,
165  iv->jac_inv,
166  iv->weighted_measure,
168  } else if((ir->getType() == panzer::IntegrationDescriptor::CV_VOLUME)
169  or (ir->getType() == panzer::IntegrationDescriptor::CV_SIDE)
170  or (ir->getType() == panzer::IntegrationDescriptor::CV_BOUNDARY)){
171  bv->evaluateValuesCV(iv->ref_ip_coordinates,
172  iv->jac,
173  iv->jac_det,
174  iv->jac_inv);
175  } else {
176  bv->evaluateValues(iv->cub_points,
177  iv->jac,
178  iv->jac_det,
179  iv->jac_inv,
180  iv->weighted_measure,
182  }
183  _basis_map[basis_description.getKey()][integration_description.getKey()] = bv;
184  }
185 
186  }
187 
188  // Create the point terms and basis-integrator pairs
189  for(const panzer::PointDescriptor & point_description : point_descriptors){
190 
191  // get the generaotr, and build some points from a topology
192  auto points = point_description.getGenerator().getPoints(*cell_topology);
193 
194  // Create and store integration rule
195  Teuchos::RCP<panzer::PointRule> pr = Teuchos::rcp(new panzer::PointRule(point_description,
196  cell_topology,
197  num_cells));
198  _point_rule_map[point_description.getKey()] = pr;
199 
200  // Create and store integration values
201  Teuchos::RCP<panzer::PointValues2<double> > pv = Teuchos::rcp(new panzer::PointValues2<double>("",true));
202  pv->setupArrays(pr);
203  pv->evaluateValues(cell_vertex_coordinates,points);
204 
205  _point_map[point_description.getKey()] = pv;
206 
207  // We need to generate a integration rule - basis pair for each basis
208  for(const panzer::BasisDescriptor & basis_description : basis_descriptors){
209 
210  // Grab the basis that was pre-calculated
211  const Teuchos::RCP<const panzer::PureBasis> & basis = _pure_basis_map[basis_description.getKey()];
212 
213  // Create a basis ir layout for this pair of integrator and basis
214  Teuchos::RCP<panzer::BasisIRLayout> b_layout = Teuchos::rcp(new panzer::BasisIRLayout(basis,*pr));
215 
216  // Create and store basis values
217  Teuchos::RCP<panzer::BasisValues2<double> > bv = Teuchos::rcp(new panzer::BasisValues2<double>("",true,true));
218  bv->setupArrays(b_layout);
219 
220  bv->evaluateValues(pv->coords_ref,
221  pv->jac,
222  pv->jac_det,
223  pv->jac_inv);
224 
225  _basis_map[basis_description.getKey()][point_description.getKey()] = bv;
226  }
227 
228  }
229 
230 }
231 
234 {
235  TEUCHOS_ASSERT(_face_connectivity != Teuchos::null);
236  return *_face_connectivity;
237 }
238 
241 {
242  const auto itr = _integrator_map.find(description.getKey());
243  TEUCHOS_ASSERT(itr != _integrator_map.end());
244  return *(itr->second);
245 }
246 
249 {
250  const auto itr = _integration_rule_map.find(description.getKey());
251  TEUCHOS_ASSERT(itr != _integration_rule_map.end());
252  return *(itr->second);
253 }
254 
257  const panzer::IntegrationDescriptor & integration_description)
258 {
259  const auto itr = _basis_map.find(basis_description.getKey());
260  TEUCHOS_TEST_FOR_EXCEPT_MSG(itr == _basis_map.end(),
261  "Workset::getBasisValues: Can't find basis \"" + basis_description.getType() + "\" "
262  "of order " + std::to_string(basis_description.getOrder()));
263  const auto & integration_map = itr->second;
264  const auto itr2 = integration_map.find(integration_description.getKey());
265  TEUCHOS_TEST_FOR_EXCEPT_MSG(itr2 == integration_map.end(),
266  "Workset::getBasisValues: Can't find integration " + std::to_string(integration_description.getType()) + " "
267  "of order " + std::to_string(integration_description.getOrder()));
268  return *(itr2->second);
269 }
270 
273  const panzer::PointDescriptor & point_description) const
274 {
275  const auto itr = _basis_map.find(basis_description.getKey());
276  TEUCHOS_TEST_FOR_EXCEPT_MSG(itr == _basis_map.end(),
277  "Workset::getBasisValues: Can't find basis \"" + basis_description.getType() + "\" "
278  "of order " + std::to_string(basis_description.getOrder()));
279  const auto & point_map = itr->second;
280  const auto itr2 = point_map.find(point_description.getKey());
281  TEUCHOS_TEST_FOR_EXCEPT_MSG(itr2 == point_map.end(),
282  "Workset::getBasisValues: Can't find point values \"" + point_description.getType() + "\"");
283  return *(itr2->second);
284 }
285 
288  const panzer::IntegrationDescriptor & integration_description) const
289 {
290  const auto itr = _basis_map.find(basis_description.getKey());
291  TEUCHOS_TEST_FOR_EXCEPT_MSG(itr == _basis_map.end(),
292  "Workset::getBasisValues: Can't find basis \"" + basis_description.getType() + "\" "
293  "of order " + std::to_string(basis_description.getOrder()));
294  const auto & integration_map = itr->second;
295  const auto itr2 = integration_map.find(integration_description.getKey());
296  TEUCHOS_TEST_FOR_EXCEPT_MSG(itr2 == integration_map.end(),
297  "Workset::getBasisValues: Can't find integration " + std::to_string(integration_description.getType()) + " "
298  "of order " + std::to_string(integration_description.getOrder()));
299  return *(itr2->second);
300 }
301 
304 {
305  const auto itr = _point_map.find(point_description.getKey());
306  TEUCHOS_TEST_FOR_EXCEPT_MSG(itr == _point_map.end(),
307  "Workset::getPointValues: Can't find point values \"" + point_description.getType() + "\"");
308  return *(itr->second);
309 }
310 
311 const panzer::PureBasis &
313 {
314  const auto itr = _pure_basis_map.find(description.getKey());
315  TEUCHOS_ASSERT(itr != _pure_basis_map.end());
316  return *(itr->second);
317 }
318 
319  std::ostream& operator<<(std::ostream& os, const panzer::Workset& w)
320  {
321  using std::endl;
322 
323  os << "Workset" << endl;
324  os << " block_id=" << w.block_id << endl;
325  os << " num_cells:" << w.num_cells << endl;
326  os << " cell_local_ids (size=" << w.cell_local_ids.size() << ")" << endl;
327  os << " subcell_dim = " << w.subcell_dim << endl;
328  os << " subcell_index = " << w.subcell_index << endl;
329 
330  os << " ir_degrees: " << endl;
331  for (std::vector<int>::const_iterator ir = w.ir_degrees->begin();
332  ir != w.ir_degrees->end(); ++ir)
333  os << " " << *ir << std::endl;
334 
335  std::vector<int>::const_iterator ir = w.ir_degrees->begin();
336  for (std::vector<Teuchos::RCP<panzer::IntegrationValues2<double> > >::const_iterator irv = w.int_rules.begin();
337  irv != w.int_rules.end(); ++irv,++ir) {
338 
339  os << " IR Values (Degree=" << *ir << "):" << endl;
340 
341  os << " cub_points:" << endl;
342  os << (*irv)->cub_points << endl;
343 
344  os << " side_cub_points:" << endl;
345  os << (*irv)->side_cub_points << endl;
346 
347  os << " cub_weights:" << endl;
348  os << (*irv)->cub_weights << endl;
349 
350  os << " node_coordinates:" << endl;
351  os << (*irv)->node_coordinates << endl;
352 
353  os << " jac:" << endl;
354  os << (*irv)->jac << endl;
355 
356  os << " jac_inv:" << endl;
357  os << (*irv)->jac_inv << endl;
358 
359  os << " jac_det:" << endl;
360  os << (*irv)->jac_det << endl;
361 
362  os << " weighted_measure:" << endl;
363  os << (*irv)->weighted_measure << endl;
364 
365  os << " covarient:" << endl;
366  os << (*irv)->covarient << endl;
367 
368  os << " contravarient:" << endl;
369  os << (*irv)->contravarient << endl;
370 
371  os << " norm_contravarient:" << endl;
372  os << (*irv)->norm_contravarient << endl;
373 
374  os << " ip_coordinates:" << endl;
375  os << (*irv)->ip_coordinates << endl;
376 
377  os << " int_rule->getName():" << (*irv)->int_rule->getName() << endl;
378  }
379 
380 
381  os << " basis_names: " << endl;
382  for (std::vector<std::string>::const_iterator b = w.basis_names->begin();
383  b != w.basis_names->end(); ++b)
384  os << " " << *b << std::endl;
385 
386  std::vector<std::string>::const_iterator b = w.basis_names->begin();
387 
388  for (std::vector<Teuchos::RCP< panzer::BasisValues2<double> > >::const_iterator bv = w.bases.begin(); bv != w.bases.end(); ++bv,++b) {
389 
390  os << " Basis Values (basis_name=" << *b << "):" << endl;
391 
392 /*
393  os << " basis_ref:" << endl;
394  os << (*bv)->basis_ref << endl;
395 
396  os << " basis:" << endl;
397  os << (*bv)->basis_scalar << endl;
398 
399  os << " grad_basis_ref:" << endl;
400  os << (*bv)->grad_basis_ref << endl;
401 
402  os << " grad_basis:" << endl;
403  os << (*bv)->grad_basis << endl;
404 
405  os << " curl_basis_ref:" << endl;
406  os << (*bv)->curl_basis_ref_vector << endl;
407 
408  os << " curl_basis:" << endl;
409  os << (*bv)->curl_basis_vector << endl;
410 
411  os << " basis_coordinates_ref:" << endl;
412  os << (*bv)->basis_coordinates_ref << endl;
413 
414  os << " basis_coordinates:" << endl;
415  os << (*bv)->basis_coordinates << endl;
416 */
417 
418  os << " basis_layout->name():" << (*bv)->basis_layout->name() << endl;
419  }
420 
421 
422 
423  return os;
424  }
425 
426 }
std::size_t getKey() const
Get unique key associated with basis of this order and type The key is used to sort through a map of ...
std::map< size_t, std::map< size_t, Teuchos::RCP< panzer::BasisValues2< double > > > > _basis_map
const panzer::IntegrationRule & getIntegrationRule(const panzer::IntegrationDescriptor &description) const
Grab the integration rule (contains data layouts) for a given integration description (throws error i...
std::map< size_t, Teuchos::RCP< const panzer::IntegrationRule > > _integration_rule_map
const std::vector< panzer::BasisDescriptor > & getBases() const
Get a list of bases being requested.
const int & getType() const
Get type of integrator.
Integral over a specific side of cells (side must be set)
Teuchos::RCP< std::vector< int > > ir_degrees
If workset corresponds to a sub cell, what is the index?
panzer::BasisValues2< double > & getBasisValues(const panzer::BasisDescriptor &basis_description, const panzer::IntegrationDescriptor &integration_description)
Grab the basis values for a given basis description and integration description (throws error if it d...
std::vector< Teuchos::RCP< panzer::BasisValues2< double > > > bases
Static basis function data, key is basis name, value is index in the static_bases vector...
Teuchos::RCP< std::vector< std::string > > basis_names
Value corresponds to basis type. Use the offest for indexing.
const panzer::PointValues2< double > & getPointValues(const panzer::PointDescriptor &point_description) const
Grab the basis values for a given basis description and integration description (throws error if it d...
std::map< size_t, Teuchos::RCP< const panzer::PureBasis > > _pure_basis_map
Generates a SubcellConnectivity associated with faces and cells given a partition of the local mesh...
Teuchos::RCP< const shards::CellTopology > cell_topology
const std::string & getType() const
Get unique string associated with the type of point descriptor. This will be used generate a hash to ...
CellCoordArray cell_vertex_coordinates
Kokkos::View< LO * > local_cells
const std::vector< panzer::IntegrationDescriptor > & getIntegrators() const
Get a list of integrators being requested.
const std::vector< panzer::PointDescriptor > & getPoints() const
Get a list of points being requested.
const panzer::IntegrationValues2< double > & getIntegrationValues(const panzer::IntegrationDescriptor &description) const
Grab the integration values for a given integration description (throws error if integration doesn&#39;t ...
std::map< size_t, Teuchos::RCP< const panzer::PointRule > > _point_rule_map
void setupNeeds(Teuchos::RCP< const shards::CellTopology > cell_topology, const Kokkos::View< double ***, PHX::Device > &cell_vertices, const panzer::WorksetNeeds &needs)
Teuchos::RCP< panzer::SubcellConnectivity > _face_connectivity
std::map< size_t, Teuchos::RCP< const panzer::IntegrationValues2< double > > > _integrator_map
std::vector< Teuchos::RCP< panzer::IntegrationValues2< double > > > int_rules
Kokkos::View< double ***, PHX::Device > cell_vertices
Kokkos::View< const int *, PHX::Device > cell_local_ids_k
std::size_t getKey() const
Get unique key associated with integrator of this order and type The key is used to sort through a ma...
std::ostream & operator<<(std::ostream &os, const AssemblyEngineInArgs &in)
int numSubcells() const
Gives number of subcells (e.g. faces) in connectivity.
std::size_t getKey() const
Get unique key associated with integrator of this order and type The key is used to sort through a ma...
const panzer::SubcellConnectivity & getFaceConnectivity() const
Grab the face connectivity for this workset.
std::map< size_t, Teuchos::RCP< const panzer::PointValues2< double > > > _point_map
Description and data layouts associated with a particular basis.
void populateValueArrays(std::size_t num_cells, bool isSide, const WorksetNeeds &needs, WorksetDetails &details, const Teuchos::RCP< WorksetDetails > other_details)
int getOrder() const
Get order of basis.
std::vector< GO > cell_local_ids
const std::string & getType() const
Get type of basis.
const panzer::PureBasis & getBasis(const panzer::BasisDescriptor &description) const
Grab the pure basis (contains data layouts) for a given basis description (throws error if integratio...
const int & getOrder() const
Get order of integrator.
void setup(const panzer::LocalMeshPartition< int, panzer::Ordinal64 > &partition, const panzer::WorksetNeeds &needs)
Constructs the workset details from a given chunk of the mesh.