Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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
145  iv->setupArrays(ir);
146  iv->evaluateValues(cell_vertex_coordinates,num_cells);
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
157 
158  // Create and store basis values
160  bv->setupArrays(b_layout);
162  bv->evaluateValues(iv->ref_ip_coordinates,
163  iv->jac,
164  iv->jac_det,
165  iv->jac_inv,
166  iv->weighted_measure,
168  true,
169  num_cells);
173  bv->evaluateValuesCV(iv->ref_ip_coordinates,
174  iv->jac,
175  iv->jac_det,
176  iv->jac_inv);
177  } else {
178  bv->evaluateValues(iv->cub_points,
179  iv->jac,
180  iv->jac_det,
181  iv->jac_inv,
182  iv->weighted_measure,
184  true,
185  num_cells);
186  }
187  _basis_map[basis_description.getKey()][integration_description.getKey()] = bv;
188  }
189 
190  }
191 
192  // Create the point terms and basis-integrator pairs
193  for(const panzer::PointDescriptor & point_description : point_descriptors){
194 
195  // get the generaotr, and build some points from a topology
196  auto points = point_description.getGenerator().getPoints(*cell_topology);
197 
198  // Create and store integration rule
200  cell_topology,
201  num_cells));
202  _point_rule_map[point_description.getKey()] = pr;
203 
204  // Create and store integration values
206  pv->setupArrays(pr);
207  pv->evaluateValues(cell_vertex_coordinates,points);
208 
209  _point_map[point_description.getKey()] = pv;
210 
211  // We need to generate a integration rule - basis pair for each basis
212  for(const panzer::BasisDescriptor & basis_description : basis_descriptors){
213 
214  // Grab the basis that was pre-calculated
215  const Teuchos::RCP<const panzer::PureBasis> & basis = _pure_basis_map[basis_description.getKey()];
216 
217  // Create a basis ir layout for this pair of integrator and basis
219 
220  // Create and store basis values
222  bv->setupArrays(b_layout);
223 
224  bv->evaluateValues(pv->coords_ref,
225  pv->jac,
226  pv->jac_det,
227  pv->jac_inv);
228 
229  _basis_map[basis_description.getKey()][point_description.getKey()] = bv;
230  }
231 
232  }
233 
234 }
235 
238 {
239  TEUCHOS_ASSERT(_face_connectivity != Teuchos::null);
240  return *_face_connectivity;
241 }
242 
245 {
246  const auto itr = _integrator_map.find(description.getKey());
247  TEUCHOS_ASSERT(itr != _integrator_map.end());
248  return *(itr->second);
249 }
250 
253 {
254  const auto itr = _integration_rule_map.find(description.getKey());
256  return *(itr->second);
257 }
258 
261  const panzer::IntegrationDescriptor & integration_description)
262 {
263  const auto itr = _basis_map.find(basis_description.getKey());
265  "Workset::getBasisValues: Can't find basis \"" + basis_description.getType() + "\" "
266  "of order " + std::to_string(basis_description.getOrder()));
267  const auto & integration_map = itr->second;
268  const auto itr2 = integration_map.find(integration_description.getKey());
269  TEUCHOS_TEST_FOR_EXCEPT_MSG(itr2 == integration_map.end(),
270  "Workset::getBasisValues: Can't find integration " + std::to_string(integration_description.getType()) + " "
271  "of order " + std::to_string(integration_description.getOrder()));
272  return *(itr2->second);
273 }
274 
277  const panzer::PointDescriptor & point_description) const
278 {
279  const auto itr = _basis_map.find(basis_description.getKey());
281  "Workset::getBasisValues: Can't find basis \"" + basis_description.getType() + "\" "
282  "of order " + std::to_string(basis_description.getOrder()));
283  const auto & point_map = itr->second;
284  const auto itr2 = point_map.find(point_description.getKey());
285  TEUCHOS_TEST_FOR_EXCEPT_MSG(itr2 == point_map.end(),
286  "Workset::getBasisValues: Can't find point values \"" + point_description.getType() + "\"");
287  return *(itr2->second);
288 }
289 
292  const panzer::IntegrationDescriptor & integration_description) const
293 {
294  const auto itr = _basis_map.find(basis_description.getKey());
296  "Workset::getBasisValues: Can't find basis \"" + basis_description.getType() + "\" "
297  "of order " + std::to_string(basis_description.getOrder()));
298  const auto & integration_map = itr->second;
299  const auto itr2 = integration_map.find(integration_description.getKey());
300  TEUCHOS_TEST_FOR_EXCEPT_MSG(itr2 == integration_map.end(),
301  "Workset::getBasisValues: Can't find integration " + std::to_string(integration_description.getType()) + " "
302  "of order " + std::to_string(integration_description.getOrder()));
303  return *(itr2->second);
304 }
305 
308 {
309  const auto itr = _point_map.find(point_description.getKey());
311  "Workset::getPointValues: Can't find point values \"" + point_description.getType() + "\"");
312  return *(itr->second);
313 }
314 
315 const panzer::PureBasis &
317 {
318  const auto itr = _pure_basis_map.find(description.getKey());
319  TEUCHOS_ASSERT(itr != _pure_basis_map.end());
320  return *(itr->second);
321 }
322 
323  std::ostream& operator<<(std::ostream& os, const panzer::Workset& w)
324  {
325  using std::endl;
326 
327  os << "Workset" << endl;
328  os << " block_id=" << w.block_id << endl;
329  os << " num_cells:" << w.num_cells << endl;
330  os << " cell_local_ids (size=" << w.cell_local_ids.size() << ")" << endl;
331  os << " subcell_dim = " << w.subcell_dim << endl;
332  os << " subcell_index = " << w.subcell_index << endl;
333 
334  os << " ir_degrees: " << endl;
335  for (std::vector<int>::const_iterator ir = w.ir_degrees->begin();
336  ir != w.ir_degrees->end(); ++ir)
337  os << " " << *ir << std::endl;
338 
339  std::vector<int>::const_iterator ir = w.ir_degrees->begin();
340  for (std::vector<Teuchos::RCP<panzer::IntegrationValues2<double> > >::const_iterator irv = w.int_rules.begin();
341  irv != w.int_rules.end(); ++irv,++ir) {
342 
343  os << " IR Values (Degree=" << *ir << "):" << endl;
344 
345  os << " cub_points:" << endl;
346  os << (*irv)->cub_points << endl;
347 
348  os << " side_cub_points:" << endl;
349  os << (*irv)->side_cub_points << endl;
350 
351  os << " cub_weights:" << endl;
352  os << (*irv)->cub_weights << endl;
353 
354  os << " node_coordinates:" << endl;
355  os << (*irv)->node_coordinates << endl;
356 
357  os << " jac:" << endl;
358  os << (*irv)->jac << endl;
359 
360  os << " jac_inv:" << endl;
361  os << (*irv)->jac_inv << endl;
362 
363  os << " jac_det:" << endl;
364  os << (*irv)->jac_det << endl;
365 
366  os << " weighted_measure:" << endl;
367  os << (*irv)->weighted_measure << endl;
368 
369  os << " covarient:" << endl;
370  os << (*irv)->covarient << endl;
371 
372  os << " contravarient:" << endl;
373  os << (*irv)->contravarient << endl;
374 
375  os << " norm_contravarient:" << endl;
376  os << (*irv)->norm_contravarient << endl;
377 
378  os << " ip_coordinates:" << endl;
379  os << (*irv)->ip_coordinates << endl;
380 
381  os << " int_rule->getName():" << (*irv)->int_rule->getName() << endl;
382  }
383 
384 
385  os << " basis_names: " << endl;
386  for (std::vector<std::string>::const_iterator b = w.basis_names->begin();
387  b != w.basis_names->end(); ++b)
388  os << " " << *b << std::endl;
389 
390  std::vector<std::string>::const_iterator b = w.basis_names->begin();
391 
392  for (std::vector<Teuchos::RCP< panzer::BasisValues2<double> > >::const_iterator bv = w.bases.begin(); bv != w.bases.end(); ++bv,++b) {
393 
394  os << " Basis Values (basis_name=" << *b << "):" << endl;
395 
396 /*
397  os << " basis_ref:" << endl;
398  os << (*bv)->basis_ref << endl;
399 
400  os << " basis:" << endl;
401  os << (*bv)->basis_scalar << endl;
402 
403  os << " grad_basis_ref:" << endl;
404  os << (*bv)->grad_basis_ref << endl;
405 
406  os << " grad_basis:" << endl;
407  os << (*bv)->grad_basis << endl;
408 
409  os << " curl_basis_ref:" << endl;
410  os << (*bv)->curl_basis_ref_vector << endl;
411 
412  os << " curl_basis:" << endl;
413  os << (*bv)->curl_basis_vector << endl;
414 
415  os << " basis_coordinates_ref:" << endl;
416  os << (*bv)->basis_coordinates_ref << endl;
417 
418  os << " basis_coordinates:" << endl;
419  os << (*bv)->basis_coordinates << endl;
420 */
421 
422  os << " basis_layout->name():" << (*bv)->basis_layout->name() << endl;
423  }
424 
425 
426 
427  return os;
428  }
429 
430 }
std::map< size_t, std::map< size_t, Teuchos::RCP< panzer::BasisValues2< double > > > > _basis_map
std::map< size_t, Teuchos::RCP< const panzer::IntegrationRule > > _integration_rule_map
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::size_t getKey() const
Get unique key associated with integrator of this order and type The key is used to sort through a ma...
Integral over a specific side of cells (side must be set)
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 ...
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.
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 int & getType() const
Get type of integrator.
CellCoordArray cell_vertex_coordinates
Kokkos::View< LO * > local_cells
const std::vector< panzer::PointDescriptor > & getPoints() const
Get a list of points being requested.
const std::vector< panzer::BasisDescriptor > & getBases() const
Get a list of bases being requested.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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 ...
const std::string & getType() const
Get type of basis.
#define TEUCHOS_TEST_FOR_EXCEPT_MSG(throw_exception_test, msg)
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)
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...
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
int getOrder() const
Get order of basis.
KOKKOS_INLINE_FUNCTION int numSubcells() const
Gives number of subcells (e.g. faces) in connectivity.
const panzer::IntegrationRule & getIntegrationRule(const panzer::IntegrationDescriptor &description) const
Grab the integration rule (contains data layouts) for a given integration description (throws error i...
const std::vector< panzer::IntegrationDescriptor > & getIntegrators() const
Get a list of integrators being requested.
Kokkos::View< const int *, PHX::Device > cell_local_ids_k
const std::string & getType() const
Get unique string associated with the type of point descriptor. This will be used generate a hash to ...
std::ostream & operator<<(std::ostream &os, const AssemblyEngineInArgs &in)
std::map< size_t, Teuchos::RCP< const panzer::PointValues2< double > > > _point_map
Description and data layouts associated with a particular basis.
#define TEUCHOS_ASSERT(assertion_test)
void populateValueArrays(std::size_t num_cells, bool isSide, const WorksetNeeds &needs, WorksetDetails &details, const Teuchos::RCP< WorksetDetails > other_details)
const panzer::SubcellConnectivity & getFaceConnectivity() const
Grab the face connectivity for this workset.
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::vector< GO > cell_local_ids
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.