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"
56 
59 
60 #include "Shards_CellTopology.hpp"
61 
62 namespace panzer {
63 
64 namespace {
65 
66 void
67 applyBV2Orientations(const int num_cells,
68  BasisValues2<double> & basis_values,
69  const Kokkos::View<const panzer::LocalOrdinal*,PHX::Device> & local_cell_ids,
70  const Teuchos::RCP<const OrientationsInterface> & orientations_interface)
71 {
72  // This call exists because there is a middle man we have to go through.
73  // orientations_interface is a 'local' object, not a workset object so we need to map local cells to workset cells
74 
75  // If the object doesn't exist, feel free to skip applying orientations, they aren't needed in some cases (e.g. DG/FV)
76  if(orientations_interface.is_null())
77  return;
78 
79  // Ignore this operation if it has already been applied
80  if(basis_values.orientationsApplied())
81  return;
82 
83  const auto & local_orientations = *orientations_interface->getOrientations();
84  std::vector<Intrepid2::Orientation> workset_orientations(num_cells);
85 
86  auto local_cell_ids_h = Kokkos::create_mirror_view(local_cell_ids);
87  Kokkos::deep_copy(local_cell_ids_h, local_cell_ids);
88 
89  // We can only apply orientations to owned and ghost cells - virtual cells are ignored (no orientations available)
90  auto local_cell_ids_host = Kokkos::create_mirror_view(local_cell_ids);
91  Kokkos::deep_copy(local_cell_ids_host, local_cell_ids);
92  for(int i=0; i<num_cells; ++i)
93  workset_orientations[i] = local_orientations[local_cell_ids_host[i]];
94  basis_values.applyOrientations(workset_orientations,num_cells);
95 }
96 
97 void
98 correctVirtualNormals(const Teuchos::RCP<panzer::IntegrationValues2<double> > iv,
99  const int num_real_cells,
100  const int num_virtual_cells,
101  const Teuchos::RCP<const shards::CellTopology> cell_topology,
102  const SubcellConnectivity & face_connectivity)
103 {
104 
105  if (iv->int_rule->getType() != panzer::IntegrationDescriptor::SURFACE)
106  return;
107 
108  // IntegrationValues2 doesn't know anything about virtual cells, so it sets up incorrect normals for those.
109  // What we want is for the adjoining face of the virtual cell to have normals that are the negated real cell's normals.
110  // we correct the normals here:
111  const int space_dim = cell_topology->getDimension();
112  const int faces_per_cell = cell_topology->getSubcellCount(space_dim-1);
113  const int points = iv->surface_normals.extent_int(1);
114  const int points_per_face = points / faces_per_cell;
115 
116  auto surface_normals_view = PHX::as_view(iv->surface_normals);
117  auto surface_normals_h = Kokkos::create_mirror_view(surface_normals_view);
118  Kokkos::deep_copy(surface_normals_h, surface_normals_view);
119 
120  auto surface_rotation_matrices_view = PHX::as_view(iv->surface_rotation_matrices);
121  auto surface_rotation_matrices_h = Kokkos::create_mirror_view(surface_rotation_matrices_view);
122  Kokkos::deep_copy(surface_rotation_matrices_h, surface_rotation_matrices_view);
123 
124 
125  for (int virtual_cell_ordinal=0; virtual_cell_ordinal<num_virtual_cells; virtual_cell_ordinal++)
126  {
127  const panzer::LocalOrdinal virtual_cell = virtual_cell_ordinal+num_real_cells;
128  int virtual_local_face_id = -1; // the virtual cell face that adjoins the real cell
129  int face_ordinal = -1;
130  for (int local_face_id=0; local_face_id<faces_per_cell; local_face_id++)
131  {
132  face_ordinal = face_connectivity.subcellForCellHost(virtual_cell, local_face_id);
133  if (face_ordinal >= 0)
134  {
135  virtual_local_face_id = local_face_id;
136  break;
137  }
138  }
139  if (face_ordinal >= 0)
140  {
141  const int first_cell_for_face = face_connectivity.cellForSubcellHost(face_ordinal, 0);
142  const panzer::LocalOrdinal other_side = (first_cell_for_face == virtual_cell) ? 1 : 0;
143  const panzer::LocalOrdinal real_cell = face_connectivity.cellForSubcellHost(face_ordinal,other_side);
144  const panzer::LocalOrdinal face_in_real_cell = face_connectivity.localSubcellForSubcellHost(face_ordinal,other_side);
145  TEUCHOS_ASSERT(real_cell < num_real_cells);
146  for (int point_ordinal=0; point_ordinal<points_per_face; point_ordinal++)
147  {
148  int virtual_cell_point = points_per_face * virtual_local_face_id + point_ordinal;
149  int real_cell_point = points_per_face * face_in_real_cell + point_ordinal;
150  // the following arrays will be used to produce/store the rotation matrix below
151  double normal[3], transverse[3], binormal[3];
152  for(int i=0;i<3;i++)
153  {
154  normal[i]=0.;
155  transverse[i]=0.;
156  binormal[i]=0.;
157  }
158 
159  for (int d=0; d<space_dim; d++)
160  {
161  const auto n_d = surface_normals_h(real_cell,real_cell_point,d);
162  surface_normals_h(virtual_cell,virtual_cell_point,d) = -n_d;
163  normal[d] = -n_d;
164  }
165 
166  panzer::convertNormalToRotationMatrix(normal,transverse,binormal);
167 
168  for(int dim=0; dim<3; ++dim){
169  surface_rotation_matrices_h(virtual_cell,virtual_cell_point,0,dim) = normal[dim];
170  surface_rotation_matrices_h(virtual_cell,virtual_cell_point,1,dim) = transverse[dim];
171  surface_rotation_matrices_h(virtual_cell,virtual_cell_point,2,dim) = binormal[dim];
172  }
173  }
174  // clear the other normals and rotation matrices for the virtual cell:
175  for (int local_face_id=0; local_face_id<faces_per_cell; local_face_id++)
176  {
177  if (local_face_id == virtual_local_face_id) continue;
178  for (int point_ordinal=0; point_ordinal<points_per_face; point_ordinal++)
179  {
180  int point = local_face_id * points_per_face + point_ordinal;
181  for (int dim=0; dim<space_dim; dim++)
182  {
183  surface_normals_h(virtual_cell,point,dim) = 0.0;
184  }
185 
186  for(int dim1=0; dim1<3; ++dim1)
187  {
188  for(int dim2=0; dim2<3; ++dim2)
189  {
190  surface_rotation_matrices_h(virtual_cell,point,dim1,dim2) = 0;
191  }
192  }
193  }
194  }
195  }
196  }
197 
198  Kokkos::deep_copy(surface_normals_view, surface_normals_h);
199  Kokkos::deep_copy(surface_rotation_matrices_view, surface_rotation_matrices_h);
200 }
201 
202 }
203 
206  : num_cells(0)
207  , subcell_dim(-1)
208  , subcell_index(-1)
209  , ir_degrees(new std::vector<int>())
210  , basis_names(new std::vector<std::string>())
211  , setup_(false)
212  , num_owned_cells_(0)
213  , num_ghost_cells_(0)
214  , num_virtual_cells_(0)
215  , num_dimensions_(-1)
216 { }
217 
218 void
221  const WorksetOptions & options)
222 {
223 
224  num_cells = partition.local_cells.extent(0);
225  num_owned_cells_ = partition.num_owned_cells;
226  num_ghost_cells_ = partition.num_ghstd_cells;
228  options_ = options;
229 
231 
232  TEUCHOS_ASSERT(partition.cell_topology != Teuchos::null);
233  cell_topology_ = partition.cell_topology;
234 
235  num_dimensions_ = cell_topology_->getDimension();
236  subcell_dim = partition.subcell_dimension;
237  subcell_index = partition.subcell_index;
238  block_id = partition.element_block_name;
239  sideset_ = partition.sideset_name;
240 
241 
242  // Allocate and fill the local cell indexes for this workset
243  {
244  Kokkos::View<LocalOrdinal*, PHX::Device> cell_ids = Kokkos::View<LocalOrdinal*, PHX::Device>("cell_ids",num_cells);
245  Kokkos::deep_copy(cell_ids, partition.local_cells);
246  cell_local_ids_k = cell_ids;
247 
248  // DEPRECATED - only retain for backward compatability
249  auto local_cells_h = Kokkos::create_mirror_view(partition.local_cells);
250  Kokkos::deep_copy(local_cells_h, partition.local_cells);
251  cell_local_ids.resize(num_cells,-1);
252  for(int cell=0;cell<num_cells;++cell){
253  const int local_cell = local_cells_h(cell);
254  cell_local_ids[cell] = local_cell;
255  }
256  }
257 
258  // Allocate and fill the cell vertices
259  {
260  // Double check this
261  TEUCHOS_ASSERT(partition.cell_vertices.Rank == 3);
262 
263  // Grab the size of the cell vertices array
264  const int num_partition_cells = partition.cell_vertices.extent(0);
265  const int num_vertices_per_cell = partition.cell_vertices.extent(1);
266  const int num_dims_per_vertex = partition.cell_vertices.extent(2);
267 
268  // Make sure there isn't some strange problem going on
269  TEUCHOS_ASSERT(num_partition_cells == num_cells);
270  TEUCHOS_ASSERT(num_vertices_per_cell > 0);
271  TEUCHOS_ASSERT(num_dims_per_vertex > 0);
272 
273  // Allocate the worksets copy of the cell vertices
274  MDFieldArrayFactory af("",true);
275  cell_vertex_coordinates = af.buildStaticArray<double, Cell, NODE, Dim>("cell vertices", num_partition_cells, num_vertices_per_cell, num_dims_per_vertex);
276 
277  // Copy vertices over
278  const auto partition_vertices = partition.cell_vertices;
279  auto cvc = cell_vertex_coordinates.get_view();
280  Kokkos::parallel_for(num_cells, KOKKOS_LAMBDA (int i) {
281  for(int j=0;j<num_vertices_per_cell;++j)
282  for(int k=0;k<num_dims_per_vertex;++k)
283  cvc(i,j,k) = partition_vertices(i,j,k);
284  });
285  Kokkos::fence();
286  }
287 
288  // Add the subcell connectivity
289  if(partition.has_connectivity){
290  auto face_connectivity = Teuchos::rcp(new FaceConnectivity);
291  face_connectivity->setup(partition);
292  face_connectivity_ = face_connectivity;
293  }
294  // We have enough information to construct Basis/Point/Integration Values on the fly
295  setup_ = true;
296 
297 }
298 
299 bool
301 hasSubcellConnectivity(const unsigned int subcell_dimension) const
302 {
303  TEUCHOS_ASSERT(setup_);
304  return (subcell_dimension == (numDimensions() - 1)) and (not face_connectivity_.is_null());
305 }
306 
307 
308 const SubcellConnectivity &
310 getSubcellConnectivity(const unsigned int subcell_dimension) const
311 {
312  TEUCHOS_ASSERT(setup_);
313  TEUCHOS_TEST_FOR_EXCEPT_MSG(not hasSubcellConnectivity(subcell_dimension),
314  "Workset::getSubcellConnectivity : Requested subcell dimension "<<subcell_dimension<<" for a "<<num_dimensions_<<"D workset. This is not supported.");
315  // Right now we have only one option
316  return *face_connectivity_;
317 }
318 
321 {
322  TEUCHOS_ASSERT(face_connectivity_ != Teuchos::null);
323  return *face_connectivity_;
324 }
325 
329 {
330  TEUCHOS_ASSERT(setup_);
331 
332  // Check if exists
333  const auto itr = integration_values_map_.find(description.getKey());
334  if(itr != integration_values_map_.end())
335  return *(itr->second);
336 
337  // Since it doesn't exist, we need to create it
338  const unsigned int subcell_dimension = numDimensions()-1;
339  int num_faces = -1;
340  if(hasSubcellConnectivity(subcell_dimension))
341  num_faces = getSubcellConnectivity(subcell_dimension).numSubcells();
342 
343  // For now, we need to make sure the descriptor lines up with the workset
345  TEUCHOS_TEST_FOR_EXCEPT_MSG(description.getSide() != getSubcellIndex(),
346  "Workset::getIntegrationValues : Attempted to build integration values for side '"<<description.getSide()<<"', but workset is constructed for side '"<<getSubcellIndex()<<"'");
347  }
348 
349  auto ir = Teuchos::rcp(new IntegrationRule(description, cell_topology_, numCells(), num_faces));
350 
351  auto iv = Teuchos::rcp(new IntegrationValues2<double>("",true));
352  iv->setupArrays(ir);
353  iv->evaluateValues(getCellVertices(), numCells(), face_connectivity_);
354 
356 
357  // This is an advanced feature that requires changes to the workset construction
358  // Basically there needs to be a way to grab the side assembly for both "details" belonging to the same workset, which requires a refactor
359  // Alternatively, we can do this using a face connectivity object, but the refactor is more important atm.
360  TEUCHOS_ASSERT(not (options_.side_assembly_ and options_.align_side_points_));
361 
362  integration_values_map_[description.getKey()] = iv;
363  ir_degrees->push_back(iv->int_rule->cubature_degree);
364  int_rules.push_back(iv);
365 
366  return *iv;
367 
368 }
369 
373  const bool lazy_version) const
374 {
375  TEUCHOS_ASSERT(setup_);
376 
377  // Check if one exists (can be of any integration order)
378  const auto itr = basis_integration_values_map_.find(description.getKey());
379  if(itr != basis_integration_values_map_.end()){
380  for(const auto & pr : itr->second)
381  return *pr.second;
382  }
383 
384  // TODO: We currently overlap BasisIntegrationValues and BasisValues
385  // To counter this we create a fake integration rule at this point to ensure the basis values exist
386 
388 
389  // We have to have the right integrator if this is a side workset
391  TEUCHOS_ASSERT(getSubcellIndex() >= 0);
393  }
394 
395  // Now just use the other call
396  return getBasisValues(description, id, lazy_version);
397 
398 }
399 
400 
403 getBasisValues(const panzer::BasisDescriptor & basis_description,
404  const panzer::IntegrationDescriptor & integration_description,
405  const bool lazy_version) const
406 {
407  TEUCHOS_ASSERT(setup_);
408 
409  // Check if exists
410  const auto itr = basis_integration_values_map_.find(basis_description.getKey());
411  if(itr != basis_integration_values_map_.end()){
412  const auto & submap = itr->second;
413  const auto itr2 = submap.find(integration_description.getKey());
414  if(itr2 != submap.end())
415  return *(itr2->second);
416 
417  }
418 
419  // Get the integration values for this description
420  const auto & iv = getIntegrationValues(integration_description);
421  auto bir = Teuchos::rcp(new BasisIRLayout(basis_description.getType(), basis_description.getOrder(), *iv.int_rule));
422 
423  Teuchos::RCP<BasisValues2<double>> biv;
424 
425  if(lazy_version){
426 
427  // Initialized for lazy evaluation
428 
429  biv = Teuchos::rcp(new BasisValues2<double>());
430 
431  if(integration_description.getType() == IntegrationDescriptor::VOLUME)
432  biv->setupUniform(bir, iv.cub_points, iv.jac, iv.jac_det, iv.jac_inv);
433  else
434  biv->setup(bir, iv.ref_ip_coordinates, iv.jac, iv.jac_det, iv.jac_inv);
435 
436  biv->setOrientations(options_.orientations_, numOwnedCells()+numGhostCells());
437  biv->setWeightedMeasure(iv.weighted_measure);
438  biv->setCellVertexCoordinates(cell_vertex_coordinates);
439 
440  } else {
441 
442  // Standard, fully allocated version of BasisValues2
443 
444  biv = Teuchos::rcp(new BasisValues2<double>("", true, true));
445  biv->setupArrays(bir);
446  if((integration_description.getType() == IntegrationDescriptor::CV_BOUNDARY) or
447  (integration_description.getType() == IntegrationDescriptor::CV_SIDE) or
448  (integration_description.getType() == IntegrationDescriptor::CV_VOLUME)){
449 
450  biv->evaluateValuesCV(iv.ref_ip_coordinates,
451  iv.jac,
452  iv.jac_det,
453  iv.jac_inv,
454  getCellVertices(),
455  true,
456  numCells());
457  } else {
458 
459  if(integration_description.getType() == IntegrationDescriptor::VOLUME){
460 
461  // TODO: Eventually we will use the other call, however, that will be part of the BasisValues2 refactor
462  // The reason we don't do it now is that there are small differences (machine precision) that break EMPIRE testing
463  biv->evaluateValues(iv.cub_points,
464  iv.jac,
465  iv.jac_det,
466  iv.jac_inv,
467  iv.weighted_measure,
468  getCellVertices(),
469  true,
470  numCells());
471 
472  } else {
473 
474  biv->evaluateValues(iv.ref_ip_coordinates,
475  iv.jac,
476  iv.jac_det,
477  iv.jac_inv,
478  iv.weighted_measure,
479  getCellVertices(),
480  true,
481  numCells());
482  }
483  }
484 
485  applyBV2Orientations(numOwnedCells()+numGhostCells(),*biv,getLocalCellIDs(),options_.orientations_);
486 
487  }
488 
489  basis_integration_values_map_[basis_description.getKey()][integration_description.getKey()] = biv;
490  bases.push_back(biv);
491  basis_names->push_back(biv->basis_layout->name());
492 
493  return *biv;
494 
495 }
496 
497 
500 getPointValues(const panzer::PointDescriptor & description) const
501 {
502  TEUCHOS_ASSERT(setup_);
503 
504  // Check if exists
505  const auto itr = point_values_map_.find(description.getKey());
506  if(itr != point_values_map_.end())
507  return *(itr->second);
508 
509  TEUCHOS_TEST_FOR_EXCEPT_MSG(not description.hasGenerator(),
510  "Point Descriptor of type '"<<description.getType()<<"' does not have associated generator.");
511 
512  auto pr = Teuchos::rcp(new PointRule(description, cell_topology_, numCells()));
513 
514  auto pv = Teuchos::rcp(new PointValues2<double>("",true));
515 
516  pv->setupArrays(pr);
517 
518  // Point values are not necessarily set at the workset level, but can be set by evaluators
519  if(description.hasGenerator())
520  if(description.getGenerator().hasPoints(*cell_topology_))
521  pv->evaluateValues(getCellVertices(), description.getGenerator().getPoints(*cell_topology_),false, numCells());
522 
523  point_values_map_[description.getKey()] = pv;
524 
525  return *pv;
526 
527 }
528 
531 getBasisValues(const panzer::BasisDescriptor & basis_description,
532  const panzer::PointDescriptor & point_description,
533  const bool lazy_version) const
534 {
535  TEUCHOS_ASSERT(setup_);
536 
537  // Check if exists
538  const auto itr = basis_point_values_map_.find(basis_description.getKey());
539  if(itr != basis_point_values_map_.end()){
540  const auto & submap = itr->second;
541  const auto itr2 = submap.find(point_description.getKey());
542  if(itr2 != submap.end())
543  return *(itr2->second);
544 
545  }
546 
547  // Get the integration values for this description
548  const auto & pv = getPointValues(point_description);
549 
550  auto bir = Teuchos::rcp(new BasisIRLayout(basis_description.getType(), basis_description.getOrder(), *pv.point_rule));
551 
552  Teuchos::RCP<BasisValues2<double>> bpv;
553 
554  if(lazy_version){
555 
556  // Initialized for lazy evaluation
557 
558  bpv = Teuchos::rcp(new BasisValues2<double>());
559 
560  bpv->setupUniform(bir, pv.coords_ref, pv.jac, pv.jac_det, pv.jac_inv);
561 
562  bpv->setOrientations(options_.orientations_, numOwnedCells()+numGhostCells());
563  bpv->setCellVertexCoordinates(cell_vertex_coordinates);
564 
565  } else {
566 
567  // Standard fully allocated version
568 
569  bpv = Teuchos::rcp(new BasisValues2<double>("", true, false));
570  bpv->setupArrays(bir);
571  bpv->evaluateValues(pv.coords_ref,
572  pv.jac,
573  pv.jac_det,
574  pv.jac_inv,
575  numCells());
576 
577  // TODO: We call this separately due to how BasisValues2 is structured - needs to be streamlined
578  bpv->evaluateBasisCoordinates(getCellVertices(),numCells());
579 
580  applyBV2Orientations(numOwnedCells()+numGhostCells(),*bpv, getLocalCellIDs(), options_.orientations_);
581 
582  }
583 
584  basis_point_values_map_[basis_description.getKey()][point_description.getKey()] = bpv;
585 
586  return *bpv;
587 
588 }
589 
593 {
594  const auto itr = _integration_rule_map.find(description.getKey());
595  if(itr == _integration_rule_map.end()){
596 
597  // Must run setup or cell topology wont be set properly
598  TEUCHOS_ASSERT(setup_);
599 
600  // Since it doesn't exist, we need to create it
601  const unsigned int subcell_dimension = numDimensions()-1;
602  int num_faces = -1;
603  if(hasSubcellConnectivity(subcell_dimension))
604  num_faces = getSubcellConnectivity(subcell_dimension).numSubcells();
605 
606  // For now, we need to make sure the descriptor lines up with the workset
608  TEUCHOS_TEST_FOR_EXCEPT_MSG(description.getSide() != getSubcellIndex(),
609  "Workset::getIntegrationValues : Attempted to build integration values for side '"<<description.getSide()<<"', but workset is constructed for side '"<<getSubcellIndex()<<"'");
610  }
611 
612  auto ir = Teuchos::rcp(new IntegrationRule(description, cell_topology_, numCells(), num_faces));
613 
614  _integration_rule_map[description.getKey()] = ir;
615 
616  return *ir;
617  }
618  return *(itr->second);
619 }
620 
621 const panzer::PureBasis &
623 getBasis(const panzer::BasisDescriptor & description) const
624 {
625  const auto itr = _pure_basis_map.find(description.getKey());
626  if(itr == _pure_basis_map.end()){
627 
628  // Must run setup or cell topology wont be set properly
629  TEUCHOS_ASSERT(setup_);
630 
631  // Create and storethe pure basis
632  Teuchos::RCP<panzer::PureBasis> basis = Teuchos::rcp(new panzer::PureBasis(description, cell_topology_, numCells()));
633  _pure_basis_map[description.getKey()] = basis;
634  return *basis;
635  }
636  return *(itr->second);
637 }
638 
639 
640 void
642 setNumberOfCells(const int o_cells,
643  const int g_cells,
644  const int v_cells)
645 {
646  num_owned_cells_ = o_cells;
647  num_ghost_cells_ = g_cells;
648  num_virtual_cells_ = v_cells;
649  num_cells = o_cells + g_cells + v_cells;
650 }
651 
652 std::ostream&
653 operator<<(std::ostream& os,
654  const panzer::Workset& w)
655 {
656  using std::endl;
657 
658  os << "Workset" << endl;
659  os << " block_id=" << w.getElementBlock() << endl;
660  os << " num_cells:" << w.num_cells << endl;
661  os << " num_owned_cells:" << w.numOwnedCells() << endl;
662  os << " num_ghost_cells:" << w.numGhostCells() << endl;
663  os << " num_virtual_cells:" << w.numVirtualCells() << endl;
664  os << " cell_local_ids (size=" << w.getLocalCellIDs().size() << ")" << endl;
665  os << " subcell_dim = " << w.getSubcellDimension() << endl;
666  os << " subcell_index = " << w.getSubcellIndex() << endl;
667 
668  os << " ir_degrees: " << endl;
669  for (std::vector<int>::const_iterator ir = w.ir_degrees->begin();
670  ir != w.ir_degrees->end(); ++ir)
671  os << " " << *ir << std::endl;
672 
673  std::vector<int>::const_iterator ir = w.ir_degrees->begin();
674  for (std::vector<Teuchos::RCP<panzer::IntegrationValues2<double> > >::const_iterator irv = w.int_rules.begin();
675  irv != w.int_rules.end(); ++irv,++ir) {
676 
677  os << " IR Values (Degree=" << *ir << "):" << endl;
678 
679  os << " cub_points:" << endl;
680  os << (*irv)->cub_points << endl;
681 
682  os << " side_cub_points:" << endl;
683  os << (*irv)->side_cub_points << endl;
684 
685  os << " cub_weights:" << endl;
686  os << (*irv)->cub_weights << endl;
687 
688  os << " node_coordinates:" << endl;
689  os << (*irv)->node_coordinates << endl;
690 
691  os << " jac:" << endl;
692  os << (*irv)->jac << endl;
693 
694  os << " jac_inv:" << endl;
695  os << (*irv)->jac_inv << endl;
696 
697  os << " jac_det:" << endl;
698  os << (*irv)->jac_det << endl;
699 
700  os << " weighted_measure:" << endl;
701  os << (*irv)->weighted_measure << endl;
702 
703  os << " covarient:" << endl;
704  os << (*irv)->covarient << endl;
705 
706  os << " contravarient:" << endl;
707  os << (*irv)->contravarient << endl;
708 
709  os << " norm_contravarient:" << endl;
710  os << (*irv)->norm_contravarient << endl;
711 
712  os << " ip_coordinates:" << endl;
713  os << (*irv)->ip_coordinates << endl;
714 
715  os << " int_rule->getName():" << (*irv)->int_rule->getName() << endl;
716  }
717 
718 
719  os << " basis_names: " << endl;
720  for (std::vector<std::string>::const_iterator b = w.basis_names->begin();
721  b != w.basis_names->end(); ++b)
722  os << " " << *b << std::endl;
723 
724  std::vector<std::string>::const_iterator b = w.basis_names->begin();
725 
726  for (std::vector<Teuchos::RCP< panzer::BasisValues2<double> > >::const_iterator bv = w.bases.begin(); bv != w.bases.end(); ++bv,++b) {
727 
728  os << " Basis Values (basis_name=" << *b << "):" << endl;
729 
730 /*
731  os << " basis_ref:" << endl;
732  os << (*bv)->basis_ref << endl;
733 
734  os << " basis:" << endl;
735  os << (*bv)->basis_scalar << endl;
736 
737  os << " grad_basis_ref:" << endl;
738  os << (*bv)->grad_basis_ref << endl;
739 
740  os << " grad_basis:" << endl;
741  os << (*bv)->grad_basis << endl;
742 
743  os << " curl_basis_ref:" << endl;
744  os << (*bv)->curl_basis_ref_vector << endl;
745 
746  os << " curl_basis:" << endl;
747  os << (*bv)->curl_basis_vector << endl;
748 
749  os << " basis_coordinates_ref:" << endl;
750  os << (*bv)->basis_coordinates_ref << endl;
751 
752  os << " basis_coordinates:" << endl;
753  os << (*bv)->basis_coordinates << endl;
754 */
755 
756  os << " basis_layout->name():" << (*bv)->basis_layout->name() << endl;
757  }
758 
759 
760 
761  return os;
762 }
763 
764 }
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 ...
PHX::MDField< Scalar, T0 > buildStaticArray(const std::string &str, int d0) const
int num_cells
DEPRECATED - use: numCells()
Teuchos::RCP< const shards::CellTopology > cell_topology
const panzer::IntegrationRule & getIntegrationRule(const panzer::IntegrationDescriptor &description) const
Grab the integration rule for a given integration description (throws error if integration doesn&#39;t ex...
const SubcellConnectivity & getSubcellConnectivity(const unsigned int subcell_dimension) const
Get the subcell connectivity for the workset topology.
std::map< size_t, Teuchos::RCP< const panzer::IntegrationRule > > _integration_rule_map
const int & getSide() const
Get side associated with integration - this is for backward compatibility.
PHX::View< const int * > cell_local_ids_k
bool hasGenerator() const
Check if the point descriptor has a generator for generating point values.
virtual Kokkos::DynRankView< double > getPoints(const shards::CellTopology &topo) const =0
Get the points for a particular topology.
Teuchos::RCP< panzer::SubcellConnectivity > face_connectivity_
KOKKOS_INLINE_FUNCTION int numSubcells() const
Gives number of subcells (e.g. faces) in connectivity.
int numVirtualCells() const
Number of cells not owned by any workset - these are used for boundary conditions.
std::map< size_t, std::map< size_t, Teuchos::RCP< panzer::BasisValues2< double > > > > basis_point_values_map_
const int & getType() const
Get type of integrator.
std::vector< size_t > cell_local_ids
Integral over a specific side of cells (side must be set)
Used to define options for lazy evaluation of BasisValues and IntegrationValues objects.
PHX::View< panzer::LocalOrdinal * > local_cells
Teuchos::RCP< std::vector< int > > ir_degrees
If workset corresponds to a sub cell, what is the index?
std::vector< Teuchos::RCP< panzer::BasisValues2< double > > > bases
Static basis function data, key is basis name, value is index in the static_bases vector...
std::map< size_t, Teuchos::RCP< const panzer::IntegrationValues2< double > > > integration_values_map_
Teuchos::RCP< std::vector< std::string > > basis_names
Value corresponds to basis type. Use the offest for indexing.
PHX::View< double *** > cell_vertices
int subcell_dim
DEPRECATED - use: getSubcellDimension()
unsigned int numDimensions() const
Get the cell dimension for the mesh.
std::map< size_t, Teuchos::RCP< const panzer::PureBasis > > _pure_basis_map
const std::string & getElementBlock() const
Get the element block id.
int numCells() const
Number of total cells in workset (owned, ghost, and virtual)
No integral specified - default state.
Generates a SubcellConnectivity associated with faces and cells given a partition of the local mesh...
const std::string & getType() const
Get unique string associated with the type of point descriptor. This will be used generate a hash to ...
panzer::LocalOrdinal num_owned_cells
CellCoordArray cell_vertex_coordinates
DEPRECATED - use: getCellVertices()
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...
void setNumberOfCells(const int owned_cells, const int ghost_cells, const int virtual_cells)
Provides access to set numbers of cells (required for backwards compatibility)
panzer::LocalOrdinal num_virtual_cells
Teuchos::RCP< const shards::CellTopology > cell_topology_
const panzer::IntegrationValues2< double > & getIntegrationValues(const panzer::IntegrationDescriptor &description) const
Get the integration values for a given integration description.
Teuchos::RCP< const OrientationsInterface > orientations_
Must be set to apply orientations - if it is set, then orientations will be applied to basis values...
int getSubcellIndex() const
Get the subcell index (returns -1 if not a subcell)
bool side_assembly_
Build integration values for sides.
std::vector< Teuchos::RCP< panzer::IntegrationValues2< double > > > int_rules
std::map< size_t, Teuchos::RCP< panzer::PointValues2< double > > > point_values_map_
const PointGenerator & getGenerator() const
const panzer::BasisValues2< double > & getBasisValues(const panzer::BasisDescriptor &basis_description, const bool lazy_version=false) const
bool hasSubcellConnectivity(const unsigned int subcell_dimension) const
Check if subcell connectivity exists for a given dimension.
KOKKOS_INLINE_FUNCTION void convertNormalToRotationMatrix(const Scalar normal[3], Scalar transverse[3], Scalar binormal[3])
int getSubcellDimension() const
Get the subcell dimension.
panzer::LocalOrdinal num_ghstd_cells
Integral over all sides of cells (closed surface integral)
std::map< size_t, std::map< size_t, Teuchos::RCP< panzer::BasisValues2< double > > > > basis_integration_values_map_
Kokkos::View< const int *, PHX::Device > getLocalCellIDs() const
Get the local cell IDs for the workset.
std::string block_id
DEPRECATED - use: getElementBlock()
WorksetDetails()
Default constructor.
int subcell_index
DEPRECATED - use: getSubcellIndex()
int numOwnedCells() const
Number of cells owned by this workset.
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)
bool align_side_points_
If workset side integration values must align with another workset, there must be a unique order assi...
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
int numGhostCells() const
Number of cells owned by a different workset.
virtual bool hasPoints(const shards::CellTopology &topo) const =0
Check if the generator can generate points for the given topology.
Description and data layouts associated with a particular basis.
CellCoordArray getCellVertices() const
Get the vertices for the cells.
void setup(const LocalMeshPartition &partition, const WorksetOptions &options)
Constructs the workset details from a given chunk of the mesh.
int getOrder() const
Get order of basis.
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...