Panzer  Version of the Day
Panzer_LocalPartitioningUtilities.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 
44 
45 #include "Teuchos_RCP.hpp"
46 #include "Teuchos_Comm.hpp"
47 #include "Teuchos_Assert.hpp"
48 
49 #include "Tpetra_CrsMatrix.hpp"
50 #include "Tpetra_RowMatrixTransposer.hpp"
51 
52 #include "Panzer_FaceToElement.hpp"
53 #include "Panzer_ConnManager.hpp"
54 #include "Panzer_NodeType.hpp"
55 #include "Panzer_FieldPattern.hpp"
60 
61 #include "Panzer_Workset_Builder.hpp"
63 
64 #include "Phalanx_KokkosDeviceTypes.hpp"
65 
66 #include <set>
67 #include <unordered_set>
68 #include <unordered_map>
69 
70 namespace panzer
71 {
72 
73 namespace
74 {
75 
79 void
80 buildCellGlobalIDs(panzer::ConnManager & conn,
81  PHX::View<panzer::GlobalOrdinal*> & globals)
82 {
83  // extract topologies, and build global connectivity...currently assuming only one topology
84  std::vector<shards::CellTopology> elementBlockTopologies;
85  conn.getElementBlockTopologies(elementBlockTopologies);
86  const shards::CellTopology & topology = elementBlockTopologies[0];
87 
88  // FIXME: We assume that all element blocks have the same topology.
89  for(const auto & other_topology : elementBlockTopologies){
90  TEUCHOS_ASSERT(other_topology.getKey() == topology.getKey());
91  }
92 
93  Teuchos::RCP<panzer::FieldPattern> cell_pattern;
94  if(topology.getDimension() == 1){
95  cell_pattern = Teuchos::rcp(new panzer::EdgeFieldPattern(topology));
96  } else if(topology.getDimension() == 2){
97  cell_pattern = Teuchos::rcp(new panzer::FaceFieldPattern(topology));
98  } else if(topology.getDimension() == 3){
99  cell_pattern = Teuchos::rcp(new panzer::ElemFieldPattern(topology));
100  }
101 
102 // panzer::EdgeFieldPattern cell_pattern(elementBlockTopologies[0]);
103  conn.buildConnectivity(*cell_pattern);
104 
105  // calculate total number of local cells
106  std::vector<std::string> block_ids;
107  conn.getElementBlockIds(block_ids);
108 
109  std::size_t totalSize = 0;
110  for (std::size_t which_blk=0;which_blk<block_ids.size();which_blk++) {
111  // get the elem to face mapping
112  const std::vector<int> & localIDs = conn.getElementBlock(block_ids[which_blk]);
113  totalSize += localIDs.size();
114  }
115  globals = PHX::View<panzer::GlobalOrdinal*>("global_cells",totalSize);
116  auto globals_h = Kokkos::create_mirror_view(globals);
117 
118  for (std::size_t id=0;id<totalSize; ++id) {
119  // sanity check
120  int n_conn = conn.getConnectivitySize(id);
121  TEUCHOS_ASSERT(n_conn==1);
122 
123  const panzer::GlobalOrdinal * connectivity = conn.getConnectivity(id);
124  globals_h(id) = connectivity[0];
125  }
126  Kokkos::deep_copy(globals, globals_h);
127 
128 // print_view_1D("buildCellGlobalIDs : globals",globals);
129 }
130 
134 void
135 buildCellToNodes(panzer::ConnManager & conn, PHX::View<panzer::GlobalOrdinal**> & globals)
136 {
137  // extract topologies, and build global connectivity...currently assuming only one topology
138  std::vector<shards::CellTopology> elementBlockTopologies;
139  conn.getElementBlockTopologies(elementBlockTopologies);
140  const shards::CellTopology & topology = elementBlockTopologies[0];
141 
142  // FIXME: We assume that all element blocks have the same topology.
143  for(const auto & other_topology : elementBlockTopologies){
144  TEUCHOS_ASSERT(other_topology.getKey() == topology.getKey());
145  }
146 
147  panzer::NodalFieldPattern pattern(topology);
148  conn.buildConnectivity(pattern);
149 
150  // calculate total number of local cells
151  std::vector<std::string> block_ids;
152  conn.getElementBlockIds(block_ids);
153 
154  // compute total cells and maximum nodes
155  std::size_t totalCells=0, maxNodes=0;
156  for (std::size_t which_blk=0;which_blk<block_ids.size();which_blk++) {
157  // get the elem to face mapping
158  const std::vector<int> & localIDs = conn.getElementBlock(block_ids[which_blk]);
159  if ( localIDs.size() == 0 )
160  continue;
161  int thisSize = conn.getConnectivitySize(localIDs[0]);
162 
163  totalCells += localIDs.size();
164  maxNodes = maxNodes<Teuchos::as<std::size_t>(thisSize) ? Teuchos::as<std::size_t>(thisSize) : maxNodes;
165  }
166  globals = PHX::View<panzer::GlobalOrdinal**>("cell_to_node",totalCells,maxNodes);
167  auto globals_h = Kokkos::create_mirror_view(globals);
168 
169  // build connectivity array
170  for (std::size_t id=0;id<totalCells; ++id) {
171  const panzer::GlobalOrdinal * connectivity = conn.getConnectivity(id);
172  int nodeCnt = conn.getConnectivitySize(id);
173 
174  for(int n=0;n<nodeCnt;n++)
175  globals_h(id,n) = connectivity[n];
176  }
177  Kokkos::deep_copy(globals, globals_h);
178 
179 // print_view("buildCellToNodes : globals",globals);
180 }
181 
182 Teuchos::RCP<const Tpetra::Map<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> >
183 buildNodeMap(const Teuchos::RCP<const Teuchos::Comm<int> > & comm,
184  PHX::View<const panzer::GlobalOrdinal**> cells_to_nodes)
185 {
186  using Teuchos::RCP;
187  using Teuchos::rcp;
188 
189  /*
190 
191  This function converts
192 
193  cells_to_nodes(local cell, local node) = global node index
194 
195  to a map describing which global nodes are found on this process
196 
197  Note that we have to ensure that that the global nodes found on this process are unique.
198 
199  */
200 
201  typedef Tpetra::Map<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> map_type;
202 
203  // get locally unique global ids
204  auto cells_to_nodes_h = Kokkos::create_mirror_view(cells_to_nodes);
205  Kokkos::deep_copy(cells_to_nodes_h, cells_to_nodes);
206  std::set<panzer::GlobalOrdinal> global_nodes;
207  for(unsigned int i=0;i<cells_to_nodes.extent(0);i++)
208  for(unsigned int j=0;j<cells_to_nodes.extent(1);j++)
209  global_nodes.insert(cells_to_nodes_h(i,j));
210 
211  // build local vector contribution
212  PHX::View<panzer::GlobalOrdinal*> node_ids("global_nodes",global_nodes.size());
213  auto node_ids_h = Kokkos::create_mirror_view(node_ids);
214  int i = 0;
215  for(auto itr=global_nodes.begin();itr!=global_nodes.end();++itr,++i)
216  node_ids_h(i) = *itr;
217  Kokkos::deep_copy(node_ids, node_ids_h);
218 
219 // print_view("buildNodeMap : cells_to_nodes",cells_to_nodes);
220 // print_view_1D("buildNodeMap : node_ids",node_ids);
221 
222  return rcp(new map_type(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),node_ids,0,comm));
223 }
224 
228 Teuchos::RCP<Tpetra::CrsMatrix<panzer::LocalOrdinal,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> >
229 buildNodeToCellMatrix(const Teuchos::RCP<const Teuchos::Comm<int> > & comm,
230  PHX::View<const panzer::GlobalOrdinal*> owned_cells,
231  PHX::View<const panzer::GlobalOrdinal**> owned_cells_to_nodes)
232 {
233  using Teuchos::RCP;
234  using Teuchos::rcp;
235 
236  typedef Tpetra::Map<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> map_type;
237  typedef Tpetra::CrsMatrix<panzer::LocalOrdinal,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> crs_type;
238  typedef Tpetra::Import<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> import_type;
239 
240 
241  PANZER_FUNC_TIME_MONITOR_DIFF("panzer_stk::buildNodeToCellMatrix",BNTCM);
242 
243  TEUCHOS_ASSERT(owned_cells.extent(0)==owned_cells_to_nodes.extent(0));
244 
245  // build a unque node map to use with fill complete
246 
247  // This map identifies all nodes linked to cells on this process
248  auto node_map = buildNodeMap(comm,owned_cells_to_nodes);
249 
250  // This map identifies nodes owned by this process
251  auto unique_node_map = Tpetra::createOneToOne(node_map);
252 
253  // This map identifies the cells owned by this process
254  RCP<map_type> cell_map = rcp(new map_type(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),owned_cells,0,comm));
255 
256  // Create a CRS matrix that stores a pointless value for every global node that belongs to a global cell
257  // This is essentially another way to store cells_to_nodes
258  RCP<crs_type> cell_to_node;
259  {
260  PANZER_FUNC_TIME_MONITOR_DIFF("Build matrix",BuildMatrix);
261 
262  // fill in the cell to node matrix
263  const unsigned int num_local_cells = owned_cells_to_nodes.extent(0);
264  const unsigned int num_nodes_per_cell = owned_cells_to_nodes.extent(1);
265 
266  // The matrix is indexed by (global cell, global node) = local node
267  cell_to_node = rcp(new crs_type(cell_map,num_nodes_per_cell,
268  Tpetra::StaticProfile));
269 
270  std::vector<panzer::LocalOrdinal> local_node_indexes(num_nodes_per_cell);
271  std::vector<panzer::GlobalOrdinal> global_node_indexes(num_nodes_per_cell);
272  auto owned_cells_h = Kokkos::create_mirror_view(owned_cells);
273  auto owned_cells_to_nodes_h = Kokkos::create_mirror_view(owned_cells_to_nodes);
274  Kokkos::deep_copy(owned_cells_h, owned_cells);
275  Kokkos::deep_copy(owned_cells_to_nodes_h, owned_cells_to_nodes);
276  for(unsigned int i=0;i<num_local_cells;i++) {
277  const panzer::GlobalOrdinal global_cell_index = owned_cells_h(i);
278  for(unsigned int j=0;j<num_nodes_per_cell;j++) {
279  local_node_indexes[j] = Teuchos::as<panzer::LocalOrdinal>(j);
280  global_node_indexes[j] = owned_cells_to_nodes_h(i,j);
281  }
282 
283  // cell_to_node_mat->insertGlobalValues(cells(i),cols,vals);
284  cell_to_node->insertGlobalValues(global_cell_index,global_node_indexes,local_node_indexes);
285  }
286  cell_to_node->fillComplete(unique_node_map,cell_map);
287 
288  }
289 
290  // Transpose the cell_to_node array to create the node_to_cell array
291  RCP<crs_type> node_to_cell;
292  {
293  PANZER_FUNC_TIME_MONITOR_DIFF("Tranpose matrix",TransposeMatrix);
294  // Create an object designed to transpose the (global cell, global node) matrix to give
295  // a (global node, global cell) matrix
296  Tpetra::RowMatrixTransposer<panzer::LocalOrdinal,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> transposer(cell_to_node);
297 
298  // Create the transpose crs matrix
299  auto trans = transposer.createTranspose();
300 
301  // We want to import the portion of the transposed matrix relating to all nodes on this process
302  // The importer must import nodes required by this process (node_map) from the unique nodes (nodes living on a process)
303  RCP<import_type> import = rcp(new import_type(unique_node_map,node_map));
304 
305  // Create the crs matrix to store (ghost node, global cell) array
306  // This CRS matrix will have overlapping rows for shared global nodes
307  node_to_cell = rcp(new crs_type(node_map,0));
308 
309  // Transfer data from the transpose array (associated with unique_node_map) to node_to_cell (associated with node_map)
310  node_to_cell->doImport(*trans,*import,Tpetra::INSERT);
311 
312  // Set the fill - basicially locks down the structured of the CRS matrix - required before doing some operations
313  //node_to_cell->fillComplete();
314  node_to_cell->fillComplete(cell_map,unique_node_map);
315  }
316 
317  // Return the node to cell array
318  return node_to_cell;
319 }
320 
323 PHX::View<panzer::GlobalOrdinal*>
324 buildGhostedCellOneRing(const Teuchos::RCP<const Teuchos::Comm<int> > & comm,
325  PHX::View<const panzer::GlobalOrdinal*> cells,
326  PHX::View<const panzer::GlobalOrdinal**> cells_to_nodes)
327 {
328 
329  PANZER_FUNC_TIME_MONITOR_DIFF("panzer_stk::buildGhostedCellOneRing",BGCOR);
330  typedef Tpetra::CrsMatrix<int,int,panzer::GlobalOrdinal,panzer::TpetraNodeType> crs_type;
331 
332  // cells : (local cell index) -> global cell index
333  // cells_to_nodes : (local cell index, local node_index) -> global node index
334 
335  /*
336 
337  This function creates a list of global indexes relating to ghost cells
338 
339  It is a little misleading how it does this, but the idea is to use the indexing of a CRS matrix to identify
340  what global cells are connected to which global nodes. The values in the CRS matrix are meaningless, however,
341  the fact that they are filled allows us to ping what index combinations exist.
342 
343  To do this we are going to use cell 'nodes' which could also be cell 'vertices'. It is unclear.
344 
345  First we construct an array that stores that global node indexes make up the connectivity for a given global cell index (order doesn't matter)
346 
347  cell_to_node : cell_to_node[global cell index][global node index] = some value (not important, just has to exist)
348 
349  We are then going to transpose that array
350 
351  node_to_cell : node_to_cell[global node index][global cell index] = some value (not important, just has to exist)
352 
353  The coloring of the node_to_cell array tells us what global cells are connected to a given global node.
354 
355 
356  */
357 
358  // the node to cell matrix: Row = Global Node Id, Cell = Global Cell Id, Value = Cell Local Node Id
359  Teuchos::RCP<crs_type> node_to_cell = buildNodeToCellMatrix(comm,cells,cells_to_nodes);
360 
361  // the set of cells already known
362  std::unordered_set<panzer::GlobalOrdinal> unique_cells;
363 
364  // mark all owned cells as already known, e.g. and not in the list of
365  // ghstd cells to be constructed
366  auto cells_h = Kokkos::create_mirror_view(cells);
367  Kokkos::deep_copy(cells_h, cells);
368  for(size_t i=0;i<cells.extent(0);i++) {
369  unique_cells.insert(cells_h(i));
370  }
371 
372  // The set of ghost cells that share a global node with an owned cell
373  std::set<panzer::GlobalOrdinal> ghstd_cells_set;
374 
375  // Get a list of global node indexes associated with the cells owned by this process
376 // auto node_map = node_to_cell->getRangeMap()->getMyGlobalIndices();
377  auto node_map = node_to_cell->getMap()->getMyGlobalIndices();
378 
379  // Iterate through the global node indexes associated with this process
380  for(size_t i=0;i<node_map.extent(0);i++) {
381  const panzer::GlobalOrdinal global_node_index = node_map(i);
382  size_t numEntries = node_to_cell->getNumEntriesInGlobalRow(node_map(i));
383  typename crs_type::nonconst_global_inds_host_view_type indices("indices", numEntries);
384  typename crs_type::nonconst_values_host_view_type values("values", numEntries);
385 
386  // Copy the row for a global node index into a local vector
387  node_to_cell->getGlobalRowCopy(global_node_index,indices,values,numEntries);
388 
389  for(size_t i=0; i<indices.extent(0); ++i) {
390  auto index = indices(i);
391  // if this is a new index (not owned, not previously found ghstd index
392  // add it to the list of ghstd cells
393  if(unique_cells.find(index)==unique_cells.end()) {
394  ghstd_cells_set.insert(index);
395  unique_cells.insert(index); // make sure you don't find it again
396  }
397  }
398  }
399 
400  // build an array containing only the ghstd cells
401  int indx = 0;
402  PHX::View<panzer::GlobalOrdinal*> ghstd_cells("ghstd_cells",ghstd_cells_set.size());
403  auto ghstd_cells_h = Kokkos::create_mirror_view(ghstd_cells);
404  for(auto global_cell_index : ghstd_cells_set) {
405  ghstd_cells_h(indx) = global_cell_index;
406  indx++;
407  }
408 
409 // print_view_1D("ghstd_cells",ghstd_cells);
410  Kokkos::deep_copy(ghstd_cells, ghstd_cells_h);
411  return ghstd_cells;
412 }
413 
414 }
415 
416 namespace partitioning_utilities
417 {
418 
419 
420 void
422  const std::vector<panzer::LocalOrdinal> & owned_parent_cells,
423  panzer::LocalMeshInfoBase & sub_info)
424 {
425  using GO = panzer::GlobalOrdinal;
426  using LO = panzer::LocalOrdinal;
427 
428  PANZER_FUNC_TIME_MONITOR_DIFF("panzer::partitioning_utilities::setupSubLocalMeshInfo",setupSLMI);
429  // The goal of this function is to fill a LocalMeshInfoBase (sub_info) with
430  // a subset of cells from a given parent LocalMeshInfoBase (parent_info)
431 
432  // Note: owned_parent_cells are the owned cells for sub_info in the parent_info's indexing scheme
433  // We need to generate sub_info's ghosts and figure out the virtual cells
434 
435  // Note: We will only handle a single ghost layer
436 
437  // Note: We assume owned_parent_cells are owned cells of the parent
438  // i.e. owned_parent_indexes cannot refer to ghost or virtual cells in parent_info
439 
440  // Note: This function works with inter-face connectivity. NOT node connectivity.
441 
442  const int num_owned_cells = owned_parent_cells.size();
443  TEUCHOS_TEST_FOR_EXCEPT_MSG(num_owned_cells == 0, "panzer::partitioning_utilities::setupSubLocalMeshInfo : Input parent subcells must exist (owned_parent_cells)");
444 
445  const int num_parent_owned_cells = parent_info.num_owned_cells;
446  TEUCHOS_TEST_FOR_EXCEPT_MSG(num_parent_owned_cells <= 0, "panzer::partitioning_utilities::setupSubLocalMeshInfo : Input parent info must contain owned cells");
447 
448  const int num_parent_ghstd_cells = parent_info.num_ghstd_cells;
449 
450  const int num_parent_total_cells = parent_info.num_owned_cells + parent_info.num_ghstd_cells + parent_info.num_virtual_cells;
451 
452  // Just as a precaution, make sure the parent_info is setup properly
453  TEUCHOS_ASSERT(static_cast<int>(parent_info.cell_to_faces.extent(0)) == num_parent_total_cells);
454  const int num_faces_per_cell = parent_info.cell_to_faces.extent(1);
455 
456  // The first thing to do is construct a vector containing the parent cell indexes of all
457  // owned, ghstd, and virtual cells
458  std::vector<LO> ghstd_parent_cells;
459  std::vector<LO> virtual_parent_cells;
460  {
461  PANZER_FUNC_TIME_MONITOR_DIFF("Construct parent cell vector",ParentCell);
462  // We grab all of the owned cells and put their global indexes into sub_info
463  // We also put all of the owned cell indexes in the parent's indexing scheme into a set to use for lookups
464  std::unordered_set<LO> owned_parent_cells_set(owned_parent_cells.begin(), owned_parent_cells.end());
465 
466  // We need to create a list of ghstd and virtual cells
467  // We do this by running through sub_cell_indexes
468  // and looking at the neighbors to find neighbors that are not owned
469 
470  // Virtual cells are defined as cells with indexes outside of the range of owned_cells and ghstd_cells
471  const int virtual_parent_cell_offset = num_parent_owned_cells + num_parent_ghstd_cells;
472 
473  std::unordered_set<LO> ghstd_parent_cells_set;
474  std::unordered_set<LO> virtual_parent_cells_set;
475  auto cell_to_faces_h = Kokkos::create_mirror_view(parent_info.cell_to_faces);
476  auto face_to_cells_h = Kokkos::create_mirror_view(parent_info.face_to_cells);
477  Kokkos::deep_copy(cell_to_faces_h, parent_info.cell_to_faces);
478  Kokkos::deep_copy(face_to_cells_h, parent_info.face_to_cells);
479  for(int i=0;i<num_owned_cells;++i){
480  const LO parent_cell_index = owned_parent_cells[i];
481  for(int local_face_index=0;local_face_index<num_faces_per_cell;++local_face_index){
482  const LO parent_face = cell_to_faces_h(parent_cell_index, local_face_index);
483 
484  // Sidesets can have owned cells that border the edge of the domain (i.e. parent_face == -1)
485  // If we are at the edge of the domain, we can ignore this face.
486  if(parent_face < 0)
487  continue;
488 
489  // Find the side index for neighbor cell with respect to the face
490  const LO neighbor_parent_side = (face_to_cells_h(parent_face,0) == parent_cell_index) ? 1 : 0;
491 
492  // Get the neighbor cell index in the parent's indexing scheme
493  const LO neighbor_parent_cell = face_to_cells_h(parent_face, neighbor_parent_side);
494 
495  // If the face exists, then the neighbor should exist
496  TEUCHOS_ASSERT(neighbor_parent_cell >= 0);
497 
498  // We can easily check if this is a virtual cell
499  if(neighbor_parent_cell >= virtual_parent_cell_offset){
500  virtual_parent_cells_set.insert(neighbor_parent_cell);
501  } else if(neighbor_parent_cell >= num_parent_owned_cells){
502  // This is a quick check for a ghost cell
503  // This branch only exists to cut down on the number of times the next branch (much slower) is called
504  ghstd_parent_cells_set.insert(neighbor_parent_cell);
505  } else {
506  // There is still potential for this to be a ghost cell with respect to 'our' cells
507  // The only way to check this is with a super slow lookup call
508  if(owned_parent_cells_set.find(neighbor_parent_cell) == owned_parent_cells_set.end()){
509  // The neighbor cell is not owned by 'us', therefore it is a ghost
510  ghstd_parent_cells_set.insert(neighbor_parent_cell);
511  }
512  }
513  }
514  }
515 
516  // We now have a list of the owned, ghstd, and virtual cells in the parent's indexing scheme.
517  // We will take the 'unordered_set's ordering for the the sub-indexing scheme
518 
519  ghstd_parent_cells.assign(ghstd_parent_cells_set.begin(), ghstd_parent_cells_set.end());
520  virtual_parent_cells.assign(virtual_parent_cells_set.begin(), virtual_parent_cells_set.end());
521 
522  }
523 
524  const int num_ghstd_cells = ghstd_parent_cells.size();
525  const int num_virtual_cells = virtual_parent_cells.size();
526  const int num_real_cells = num_owned_cells + num_ghstd_cells;
527  const int num_total_cells = num_real_cells + num_virtual_cells;
528 
529  std::vector<std::pair<LO, LO> > all_parent_cells(num_total_cells);
530  for (std::size_t i=0; i< owned_parent_cells.size(); ++i)
531  all_parent_cells[i] = std::pair<LO, LO>(owned_parent_cells[i], i);
532 
533  for (std::size_t i=0; i< ghstd_parent_cells.size(); ++i) {
534  LO insert = owned_parent_cells.size()+i;
535  all_parent_cells[insert] = std::pair<LO, LO>(ghstd_parent_cells[i], insert);
536  }
537 
538  for (std::size_t i=0; i< virtual_parent_cells.size(); ++i) {
539  LO insert = owned_parent_cells.size()+ ghstd_parent_cells.size()+i;
540  all_parent_cells[insert] = std::pair<LO, LO>(virtual_parent_cells[i], insert);
541  }
542 
543  sub_info.num_owned_cells = owned_parent_cells.size();
544  sub_info.num_ghstd_cells = ghstd_parent_cells.size();
545  sub_info.num_virtual_cells = virtual_parent_cells.size();
546 
547  // We now have the indexing order for our sub_info
548 
549  // Just as a precaution, make sure the parent_info is setup properly
550  TEUCHOS_ASSERT(static_cast<int>(parent_info.cell_vertices.extent(0)) == num_parent_total_cells);
551  TEUCHOS_ASSERT(static_cast<int>(parent_info.local_cells.extent(0)) == num_parent_total_cells);
552  TEUCHOS_ASSERT(static_cast<int>(parent_info.global_cells.extent(0)) == num_parent_total_cells);
553 
554  const int num_vertices_per_cell = parent_info.cell_vertices.extent(1);
555  const int num_dims = parent_info.cell_vertices.extent(2);
556 
557  // Fill owned, ghstd, and virtual cells: global indexes, local indexes and vertices
558  sub_info.global_cells = PHX::View<GO*>("global_cells", num_total_cells);
559  sub_info.local_cells = PHX::View<LO*>("local_cells", num_total_cells);
560  sub_info.cell_vertices = PHX::View<double***>("cell_vertices", num_total_cells, num_vertices_per_cell, num_dims);
561  auto global_cells_h = Kokkos::create_mirror_view(sub_info.global_cells);
562  auto local_cells_h = Kokkos::create_mirror_view(sub_info.local_cells);
563  auto cell_vertices_h = Kokkos::create_mirror_view(sub_info.cell_vertices);
564  auto p_global_cells_h = Kokkos::create_mirror_view(parent_info.global_cells);
565  auto p_local_cells_h = Kokkos::create_mirror_view(parent_info.local_cells);
566  auto p_cell_vertices_h = Kokkos::create_mirror_view(parent_info.cell_vertices);
567  Kokkos::deep_copy(p_global_cells_h,parent_info.global_cells);
568  Kokkos::deep_copy(p_local_cells_h,parent_info.local_cells);
569  Kokkos::deep_copy(p_cell_vertices_h,parent_info.cell_vertices);
570 
571  for (int cell=0; cell<num_total_cells; ++cell) {
572  const LO parent_cell = all_parent_cells[cell].first;
573  global_cells_h(cell) = p_global_cells_h(parent_cell);
574  local_cells_h(cell) = p_local_cells_h(parent_cell);
575  for(int vertex=0;vertex<num_vertices_per_cell;++vertex){
576  for(int dim=0;dim<num_dims;++dim){
577  cell_vertices_h(cell,vertex,dim) = p_cell_vertices_h(parent_cell,vertex,dim);
578  }
579  }
580  }
581  Kokkos::deep_copy(sub_info.global_cells, global_cells_h);
582  Kokkos::deep_copy(sub_info.local_cells, local_cells_h);
583  Kokkos::deep_copy(sub_info.cell_vertices, cell_vertices_h);
584  // Now for the difficult part
585 
586  // We need to create a new face indexing scheme from the old face indexing scheme
587 
588  // Create an auxiliary list with all cells - note this preserves indexing
589 
590  struct face_t{
591  face_t(LO c0, LO c1, LO sc0, LO sc1)
592  {
593  cell_0=c0;
594  cell_1=c1;
595  subcell_index_0=sc0;
596  subcell_index_1=sc1;
597  }
598  LO cell_0;
599  LO cell_1;
600  LO subcell_index_0;
601  LO subcell_index_1;
602  };
603 
604 
605  // First create the faces
606  std::vector<face_t> faces;
607  {
608  PANZER_FUNC_TIME_MONITOR_DIFF("Create faces",CreateFaces);
609  // faces_set: cell_0, subcell_index_0, cell_1, subcell_index_1
610  std::unordered_map<LO,std::unordered_map<LO, std::pair<LO,LO> > > faces_set;
611 
612  std::sort(all_parent_cells.begin(), all_parent_cells.end());
613 
614  auto cell_to_faces_h = Kokkos::create_mirror_view(parent_info.cell_to_faces);
615  auto face_to_cells_h = Kokkos::create_mirror_view(parent_info.face_to_cells);
616  auto face_to_lidx_h = Kokkos::create_mirror_view(parent_info.face_to_lidx);
617  Kokkos::deep_copy(cell_to_faces_h, parent_info.cell_to_faces);
618  Kokkos::deep_copy(face_to_cells_h, parent_info.face_to_cells);
619  Kokkos::deep_copy(face_to_lidx_h, parent_info.face_to_lidx);
620  for(int owned_cell=0;owned_cell<num_owned_cells;++owned_cell){
621  const LO owned_parent_cell = owned_parent_cells[owned_cell];
622  for(int local_face=0;local_face<num_faces_per_cell;++local_face){
623  const LO parent_face = cell_to_faces_h(owned_parent_cell,local_face);
624 
625  // Skip faces at the edge of the domain
626  if(parent_face<0)
627  continue;
628 
629  // Get the cell on the other side of the face
630  const LO neighbor_side = (face_to_cells_h(parent_face,0) == owned_parent_cell) ? 1 : 0;
631 
632  const LO neighbor_parent_cell = face_to_cells_h(parent_face, neighbor_side);
633  const LO neighbor_subcell_index = face_to_lidx_h(parent_face, neighbor_side);
634 
635  // Convert parent cell index into sub cell index
636  std::pair<LO, LO> search_point(neighbor_parent_cell, 0);
637  auto itr = std::lower_bound(all_parent_cells.begin(), all_parent_cells.end(), search_point);
638 
639  TEUCHOS_TEST_FOR_EXCEPT_MSG(itr == all_parent_cells.end(), "panzer_stk::setupSubLocalMeshInfo : Neighbor cell was not found in owned, ghosted, or virtual cells");
640 
641  const LO neighbor_cell = itr->second;
642 
643  LO cell_0, cell_1, subcell_index_0, subcell_index_1;
644  if(owned_cell < neighbor_cell){
645  cell_0 = owned_cell;
646  subcell_index_0 = local_face;
647  cell_1 = neighbor_cell;
648  subcell_index_1 = neighbor_subcell_index;
649  } else {
650  cell_1 = owned_cell;
651  subcell_index_1 = local_face;
652  cell_0 = neighbor_cell;
653  subcell_index_0 = neighbor_subcell_index;
654  }
655 
656  // Add this interface to the set of faces - smaller cell index is 'left' (or '0') side of face
657  faces_set[cell_0][subcell_index_0] = std::pair<LO,LO>(cell_1, subcell_index_1);
658  }
659  }
660 
661  for(const auto & cell_pair : faces_set){
662  const LO cell_0 = cell_pair.first;
663  for(const auto & subcell_pair : cell_pair.second){
664  const LO subcell_index_0 = subcell_pair.first;
665  const LO cell_1 = subcell_pair.second.first;
666  const LO subcell_index_1 = subcell_pair.second.second;
667  faces.push_back(face_t(cell_0,cell_1,subcell_index_0,subcell_index_1));
668  }
669  }
670  }
671 
672  const int num_faces = faces.size();
673 
674  sub_info.face_to_cells = PHX::View<LO*[2]>("face_to_cells", num_faces);
675  sub_info.face_to_lidx = PHX::View<LO*[2]>("face_to_lidx", num_faces);
676  sub_info.cell_to_faces = PHX::View<LO**>("cell_to_faces", num_total_cells, num_faces_per_cell);
677  auto cell_to_faces_h = Kokkos::create_mirror_view(sub_info.cell_to_faces);
678  auto face_to_cells_h = Kokkos::create_mirror_view(sub_info.face_to_cells);
679  auto face_to_lidx_h = Kokkos::create_mirror_view(sub_info.face_to_lidx);
680 
681 
682  // Default the system with invalid cell index
683  Kokkos::deep_copy(cell_to_faces_h, -1);
684 
685  for(int face_index=0;face_index<num_faces;++face_index){
686  const face_t & face = faces[face_index];
687 
688  face_to_cells_h(face_index,0) = face.cell_0;
689  face_to_cells_h(face_index,1) = face.cell_1;
690 
691  cell_to_faces_h(face.cell_0,face.subcell_index_0) = face_index;
692  cell_to_faces_h(face.cell_1,face.subcell_index_1) = face_index;
693 
694  face_to_lidx_h(face_index,0) = face.subcell_index_0;
695  face_to_lidx_h(face_index,1) = face.subcell_index_1;
696 
697  }
698  Kokkos::deep_copy(sub_info.cell_to_faces, cell_to_faces_h);
699  Kokkos::deep_copy(sub_info.face_to_cells, face_to_cells_h);
700  Kokkos::deep_copy(sub_info.face_to_lidx, face_to_lidx_h);
701  // Complete.
702 
703 }
704 
705 void
707  const int splitting_size,
708  std::vector<panzer::LocalMeshPartition> & partitions)
709 {
710 
711  using LO = panzer::LocalOrdinal;
712 
713  // Make sure the splitting size makes sense
714  TEUCHOS_ASSERT((splitting_size > 0) or (splitting_size == WorksetSizeType::ALL_ELEMENTS));
715 
716  // Default partition size
717  const LO base_partition_size = std::min(mesh_info.num_owned_cells, (splitting_size > 0) ? splitting_size : mesh_info.num_owned_cells);
718 
719  // Cells to partition
720  std::vector<LO> partition_cells;
721  partition_cells.resize(base_partition_size);
722 
723  // Create the partitions
724  LO cell_count = 0;
725  while(cell_count < mesh_info.num_owned_cells){
726 
727  LO partition_size = base_partition_size;
728  if(cell_count + partition_size > mesh_info.num_owned_cells)
729  partition_size = mesh_info.num_owned_cells - cell_count;
730 
731  // Error check for a null partition - this should never happen by design
732  TEUCHOS_ASSERT(partition_size != 0);
733 
734  // In the final partition, we need to reduce the size of partition_cells
735  if(partition_size != base_partition_size)
736  partition_cells.resize(partition_size);
737 
738  // Set the partition indexes - not really a partition, just a chunk of cells
739  for(LO i=0; i<partition_size; ++i)
740  partition_cells[i] = cell_count+i;
741 
742  // Create an empty partition
743  partitions.push_back(panzer::LocalMeshPartition());
744 
745  // Fill the empty partition
746  partitioning_utilities::setupSubLocalMeshInfo(mesh_info,partition_cells,partitions.back());
747 
748  // Update the cell count
749  cell_count += partition_size;
750  }
751 
752 }
753 
754 }
755 
756 void
758  const panzer::WorksetDescriptor & description,
759  std::vector<panzer::LocalMeshPartition> & partitions)
760 {
761  // We have to make sure that the partitioning is possible
762  TEUCHOS_ASSERT(description.getWorksetSize() != panzer::WorksetSizeType::CLASSIC_MODE);
763  TEUCHOS_ASSERT(description.getWorksetSize() != 0);
764 
765  // This could just return, but it would be difficult to debug why no partitions were returned
766  TEUCHOS_ASSERT(description.requiresPartitioning());
767 
768  const std::string & element_block_name = description.getElementBlock();
769 
770  // We have two processes for in case this is a sideset or element block
771  if(description.useSideset()){
772 
773  // If the element block doesn't exist, there are no partitions to create
774  if(mesh_info.sidesets.find(element_block_name) == mesh_info.sidesets.end())
775  return;
776  const auto & sideset_map = mesh_info.sidesets.at(element_block_name);
777 
778  const std::string & sideset_name = description.getSideset();
779 
780  // If the sideset doesn't exist, there are no partitions to create
781  if(sideset_map.find(sideset_name) == sideset_map.end())
782  return;
783 
784  const panzer::LocalMeshSidesetInfo & sideset_info = sideset_map.at(sideset_name);
785 
786  // Partitioning is not important for sidesets
787  panzer::partitioning_utilities::splitMeshInfo(sideset_info, description.getWorksetSize(), partitions);
788 
789  for(auto & partition : partitions){
790  partition.sideset_name = sideset_name;
791  partition.element_block_name = element_block_name;
792  partition.cell_topology = sideset_info.cell_topology;
793  partition.has_connectivity = true;
794  }
795 
796  } else {
797 
798  // If the element block doesn't exist, there are no partitions to create
799  if(mesh_info.element_blocks.find(element_block_name) == mesh_info.element_blocks.end())
800  return;
801 
802  // Grab the element block we're interested in
803  const panzer::LocalMeshBlockInfo & block_info = mesh_info.element_blocks.at(element_block_name);
804 
806  // We only have one partition describing the entire local mesh
807  panzer::partitioning_utilities::splitMeshInfo(block_info, -1, partitions);
808  } else {
809  // We need to partition local mesh
810 
811  // FIXME: Until the above function is fixed, we will use this hack - this will lead to horrible partitions
812  panzer::partitioning_utilities::splitMeshInfo(block_info, description.getWorksetSize(), partitions);
813 
814  }
815 
816  for(auto & partition : partitions){
817  partition.element_block_name = element_block_name;
818  partition.cell_topology = block_info.cell_topology;
819  partition.has_connectivity = true;
820  }
821  }
822 
823 }
824 
825 
826 void
827 fillLocalCellIDs(const Teuchos::RCP<const Teuchos::Comm<int>> & comm,
828  panzer::ConnManager & conn,
829  PHX::View<panzer::GlobalOrdinal*> & owned_cells,
830  PHX::View<panzer::GlobalOrdinal*> & ghost_cells,
831  PHX::View<panzer::GlobalOrdinal*> & virtual_cells)
832 {
833 
834  using Teuchos::RCP;
835 
836  // build cell to node map
837  PHX::View<panzer::GlobalOrdinal**> owned_cell_to_nodes;
838  buildCellToNodes(conn, owned_cell_to_nodes);
839 
840  // Build the local to global cell ID map
841  buildCellGlobalIDs(conn, owned_cells);
842 
843  // Get ghost cells
844  ghost_cells = buildGhostedCellOneRing(comm,owned_cells,owned_cell_to_nodes);
845 
846  // Build virtual cells
847  // Note: virtual cells are currently defined by faces (only really used for FV/DG type discretizations)
848 
849  // this class comes from Mini-PIC and Matt B
850  auto faceToElement = Teuchos::rcp(new panzer::FaceToElement<panzer::LocalOrdinal,panzer::GlobalOrdinal>());
851  faceToElement->initialize(conn);
852  auto elems_by_face = faceToElement->getFaceToElementsMap();
853  auto face_to_lidx = faceToElement->getFaceToCellLocalIdxMap();
854 
855  const panzer::LocalOrdinal num_owned_cells = owned_cells.extent(0);
856 
857  // We also need to consider faces that connect to cells that do not exist, but are needed for boundary conditions
858  // We dub them virtual cell since there should be no geometry associated with them, or topology really
859  // They exist only for datastorage so that they are consistent with 'real' cells from an algorithm perspective
860 
861  // Each virtual face (face linked to a '-1' cell) requires a virtual cell (i.e. turn the '-1' into a virtual cell)
862  // Virtual cells are those that do not exist but are connected to an owned cell
863  // Note - in the future, ghosted cells will also need to connect to virtual cells at boundary conditions, but for the moment we will ignore this.
864 
865  // Iterate over all faces and identify the faces connected to a potential virtual cell
866  std::vector<int> all_boundary_faces;
867  const int num_faces = elems_by_face.extent(0);
868  auto elems_by_face_h = Kokkos::create_mirror_view(elems_by_face);
869  Kokkos::deep_copy(elems_by_face_h, elems_by_face);
870  for(int face=0;face<num_faces;++face)
871  if(elems_by_face_h(face,0) < 0 or elems_by_face_h(face,1) < 0)
872  all_boundary_faces.push_back(face);
873  const panzer::LocalOrdinal num_virtual_cells = all_boundary_faces.size();
874 
875  // Create some global indexes associated with the virtual cells
876  // Note: We are assuming that virtual cells belong to ranks and are not 'shared' - this will change later on
877  virtual_cells = PHX::View<panzer::GlobalOrdinal*>("virtual_cells",num_virtual_cells);
878  auto virtual_cells_h = Kokkos::create_mirror_view(virtual_cells);
879  {
880  PANZER_FUNC_TIME_MONITOR_DIFF("Initial global index creation",InitialGlobalIndexCreation);
881 
882  const int num_ranks = comm->getSize();
883  const int rank = comm->getRank();
884 
885  std::vector<panzer::GlobalOrdinal> owned_cell_distribution(num_ranks,0);
886  {
887  std::vector<panzer::GlobalOrdinal> my_owned_cell_distribution(num_ranks,0);
888  my_owned_cell_distribution[rank] = num_owned_cells;
889 
890  Teuchos::reduceAll(*comm,Teuchos::REDUCE_SUM, num_ranks, my_owned_cell_distribution.data(),owned_cell_distribution.data());
891  }
892 
893  std::vector<panzer::GlobalOrdinal> virtual_cell_distribution(num_ranks,0);
894  {
895  std::vector<panzer::GlobalOrdinal> my_virtual_cell_distribution(num_ranks,0);
896  my_virtual_cell_distribution[rank] = num_virtual_cells;
897 
898  Teuchos::reduceAll(*comm,Teuchos::REDUCE_SUM, num_ranks, my_virtual_cell_distribution.data(),virtual_cell_distribution.data());
899  }
900 
901  panzer::GlobalOrdinal num_global_real_cells=0;
902  for(int i=0;i<num_ranks;++i){
903  num_global_real_cells+=owned_cell_distribution[i];
904  }
905 
906  panzer::GlobalOrdinal global_virtual_start_idx = num_global_real_cells;
907  for(int i=0;i<rank;++i){
908  global_virtual_start_idx += virtual_cell_distribution[i];
909  }
910 
911  for(int i=0;i<num_virtual_cells;++i){
912  virtual_cells_h(i) = global_virtual_start_idx + panzer::GlobalOrdinal(i);
913  }
914  }
915  Kokkos::deep_copy(virtual_cells, virtual_cells_h);
916 }
917 
918 }
Backwards compatibility mode that ignores the worksetSize in the WorksetDescriptor.
std::map< std::string, std::map< std::string, LocalMeshSidesetInfo > > sidesets
Teuchos::RCP< const shards::CellTopology > cell_topology
PHX::View< panzer::LocalOrdinal * > local_cells
void fillLocalCellIDs(const Teuchos::RCP< const Teuchos::Comm< int >> &comm, panzer::ConnManager &conn, PHX::View< panzer::GlobalOrdinal *> &owned_cells, PHX::View< panzer::GlobalOrdinal *> &ghost_cells, PHX::View< panzer::GlobalOrdinal *> &virtual_cells)
Get the owned, ghost and virtual global cell ids for this process.
PHX::View< double *** > cell_vertices
void generateLocalMeshPartitions(const panzer::LocalMeshInfo &mesh_info, const panzer::WorksetDescriptor &description, std::vector< panzer::LocalMeshPartition > &partitions)
virtual LocalOrdinal getConnectivitySize(LocalOrdinal localElmtId) const =0
PHX::View< panzer::LocalOrdinal *[2]> face_to_cells
const std::string & getSideset(const int block=0) const
Get sideset name.
panzer::LocalOrdinal num_owned_cells
panzer::LocalOrdinal num_virtual_cells
Teuchos::RCP< const shards::CellTopology > cell_topology
virtual const std::vector< LocalOrdinal > & getElementBlock(const std::string &blockID) const =0
PHX::View< panzer::LocalOrdinal *[2]> face_to_lidx
bool requiresPartitioning() const
Do we need to partition the local mesh prior to generating worksets.
virtual void getElementBlockIds(std::vector< std::string > &elementBlockIds) const =0
Pure virtual base class for supplying mesh connectivity information to the DOF Manager.
virtual void getElementBlockTopologies(std::vector< shards::CellTopology > &elementBlockTopologies) const =0
int getWorksetSize() const
Get the requested workset size (default -2 (workset size is set elsewhere), -1 (largest possible work...
bool useSideset() const
This descriptor is for a side set.
panzer::LocalOrdinal num_ghstd_cells
std::map< std::string, LocalMeshBlockInfo > element_blocks
void splitMeshInfo(const panzer::LocalMeshInfoBase &mesh_info, const int splitting_size, std::vector< panzer::LocalMeshPartition > &partitions)
const std::string & getElementBlock(const int block=0) const
Get element block name.
void setupSubLocalMeshInfo(const panzer::LocalMeshInfoBase &parent_info, const std::vector< panzer::LocalOrdinal > &owned_parent_cells, panzer::LocalMeshInfoBase &sub_info)
PHX::View< panzer::GlobalOrdinal * > global_cells
virtual void buildConnectivity(const FieldPattern &fp)=0
Workset size is set to the total number of local elements in the MPI process.
PHX::View< panzer::LocalOrdinal ** > cell_to_faces
virtual const GlobalOrdinal * getConnectivity(LocalOrdinal localElmtId) const =0