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