Panzer  Version of the Day
Panzer_STK_LocalMeshUtilities.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_NodeType.hpp"
45 #include "Panzer_STK_Interface.hpp"
48 
49 #include "Panzer_HashUtils.hpp"
50 #include "Panzer_LocalMeshInfo.hpp"
52 #include "Panzer_FaceToElement.hpp"
53 
54 #include "Panzer_FieldPattern.hpp"
59 
60 #include "Panzer_ConnManager.hpp"
61 
62 #include "Phalanx_KokkosDeviceTypes.hpp"
63 
64 #include "Teuchos_Assert.hpp"
65 #include "Teuchos_OrdinalTraits.hpp"
66 
67 #include "Tpetra_Import.hpp"
68 
69 #include <string>
70 #include <map>
71 #include <vector>
72 #include <unordered_set>
73 
74 namespace panzer_stk
75 {
76 
77 // No external access
78 namespace
79 {
80 
81 
85 Kokkos::DynRankView<double,PHX::Device>
86 buildGhostedVertices(const Tpetra::Import<int,panzer::GlobalOrdinal,panzer::TpetraNodeType> & importer,
87  Kokkos::DynRankView<const double,PHX::Device> owned_vertices)
88 {
89  using Teuchos::RCP;
90  using Teuchos::rcp;
91 
92  typedef Tpetra::MultiVector<double,int,panzer::GlobalOrdinal,panzer::TpetraNodeType> mvec_type;
93  typedef typename mvec_type::dual_view_type dual_view_type;
94 
95  size_t owned_cell_cnt = importer.getSourceMap()->getNodeNumElements();
96  size_t ghstd_cell_cnt = importer.getTargetMap()->getNodeNumElements();
97  int vertices_per_cell = owned_vertices.extent(1);
98  int space_dim = owned_vertices.extent(2);
99 
100  TEUCHOS_ASSERT(owned_vertices.extent(0)==owned_cell_cnt);
101 
102  // build vertex multivector
103  RCP<mvec_type> owned_vertices_mv = rcp(new mvec_type(importer.getSourceMap(),vertices_per_cell*space_dim));
104  RCP<mvec_type> ghstd_vertices_mv = rcp(new mvec_type(importer.getTargetMap(),vertices_per_cell*space_dim));
105 
106  {
107  auto owned_vertices_view = owned_vertices_mv->getLocalViewDevice(Tpetra::Access::OverwriteAll);
108  Kokkos::parallel_for(owned_cell_cnt, KOKKOS_LAMBDA (size_t i) {
109  int l = 0;
110  for(int j=0;j<vertices_per_cell;j++)
111  for(int k=0;k<space_dim;k++,l++)
112  owned_vertices_view(i,l) = owned_vertices(i,j,k);
113  });
114  }
115 
116  // communicate ghstd vertices
117  ghstd_vertices_mv->doImport(*owned_vertices_mv,importer,Tpetra::INSERT);
118 
119  // copy multivector into ghstd vertex structure
120  Kokkos::DynRankView<double,PHX::Device> ghstd_vertices("ghstd_vertices",ghstd_cell_cnt,vertices_per_cell,space_dim);
121  {
122  auto ghstd_vertices_view = ghstd_vertices_mv->getLocalViewDevice(Tpetra::Access::ReadOnly);
123  Kokkos::parallel_for(ghstd_cell_cnt, KOKKOS_LAMBDA (size_t i) {
124  int l = 0;
125  for(int j=0;j<vertices_per_cell;j++)
126  for(int k=0;k<space_dim;k++,l++)
127  ghstd_vertices(i,j,k) = ghstd_vertices_view(i,l);
128  } );
129  Kokkos::fence();
130  }
131 
132  return ghstd_vertices;
133 } // end build ghstd vertices
134 
135 void
136 setupLocalMeshBlockInfo(const panzer_stk::STK_Interface & mesh,
137  panzer::ConnManager & conn,
138  const panzer::LocalMeshInfo & mesh_info,
139  const std::string & element_block_name,
140  panzer::LocalMeshBlockInfo & block_info)
141 {
142 
143  // This function identifies all cells in mesh_info that belong to element_block_name
144  // and creates a block_info from it.
145 
146  const int num_parent_owned_cells = mesh_info.num_owned_cells;
147 
148  // Make sure connectivity is setup for interfaces between cells
149  {
150  const shards::CellTopology & topology = *(mesh.getCellTopology(element_block_name));
151  Teuchos::RCP<panzer::FieldPattern> cell_pattern;
152  if(topology.getDimension() == 1){
153  cell_pattern = Teuchos::rcp(new panzer::EdgeFieldPattern(topology));
154  } else if(topology.getDimension() == 2){
155  cell_pattern = Teuchos::rcp(new panzer::FaceFieldPattern(topology));
156  } else if(topology.getDimension() == 3){
157  cell_pattern = Teuchos::rcp(new panzer::ElemFieldPattern(topology));
158  }
159 
160  {
161  PANZER_FUNC_TIME_MONITOR("Build connectivity");
162  conn.buildConnectivity(*cell_pattern);
163  }
164  }
165 
166  std::vector<panzer::LocalOrdinal> owned_block_cells;
167  auto local_cells_h = Kokkos::create_mirror_view(mesh_info.local_cells);
168  Kokkos::deep_copy(local_cells_h, mesh_info.local_cells);
169  for(int parent_owned_cell=0;parent_owned_cell<num_parent_owned_cells;++parent_owned_cell){
170  const panzer::LocalOrdinal local_cell = local_cells_h(parent_owned_cell);
171  const bool is_in_block = conn.getBlockId(local_cell) == element_block_name;
172 
173  if(is_in_block){
174  owned_block_cells.push_back(parent_owned_cell);
175  }
176 
177  }
178 
179  if ( owned_block_cells.size() == 0 )
180  return;
181  block_info.num_owned_cells = owned_block_cells.size();
182  block_info.element_block_name = element_block_name;
183  block_info.cell_topology = mesh.getCellTopology(element_block_name);
184  {
185  PANZER_FUNC_TIME_MONITOR("panzer::partitioning_utilities::setupSubLocalMeshInfo");
186  panzer::partitioning_utilities::setupSubLocalMeshInfo(mesh_info, owned_block_cells, block_info);
187  }
188 }
189 
190 
191 void
192 setupLocalMeshSidesetInfo(const panzer_stk::STK_Interface & mesh,
193  panzer::ConnManager& /* conn */,
194  const panzer::LocalMeshInfo & mesh_info,
195  const std::string & element_block_name,
196  const std::string & sideset_name,
197  panzer::LocalMeshSidesetInfo & sideset_info)
198 {
199 
200  // We use these typedefs to make the algorithm slightly more clear
201  using LocalOrdinal = panzer::LocalOrdinal;
202  using ParentOrdinal = int;
203  using SubcellOrdinal = short;
204 
205  // This function identifies all cells in mesh_info that belong to element_block_name
206  // and creates a block_info from it.
207 
208  // This is a list of all entities that make up the 'side'
209  std::vector<stk::mesh::Entity> side_entities;
210 
211  // Grab the side entities associated with this sideset on the mesh
212  // Note: Throws exception if element block or sideset doesn't exist
213  try{
214 
215  mesh.getAllSides(sideset_name, element_block_name, side_entities);
216 
217  } catch(STK_Interface::SidesetException & e) {
218  std::stringstream ss;
219  std::vector<std::string> sideset_names;
220  mesh.getSidesetNames(sideset_names);
221 
222  // build an error message
223  ss << e.what() << "\nChoose existing sideset:\n";
224  for(const auto & optional_sideset_name : sideset_names){
225  ss << "\t\"" << optional_sideset_name << "\"\n";
226  }
227 
228  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
229 
230  } catch(STK_Interface::ElementBlockException & e) {
231  std::stringstream ss;
232  std::vector<std::string> element_block_names;
233  mesh.getElementBlockNames(element_block_names);
234 
235  // build an error message
236  ss << e.what() << "\nChoose existing element block:\n";
237  for(const auto & optional_element_block_name : element_block_names){
238  ss << "\t\"" << optional_element_block_name << "\"\n";
239  }
240 
241  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
242 
243  } catch(std::logic_error & e) {
244  std::stringstream ss;
245  ss << e.what() << "\nUnrecognized logic error.\n";
246 
247  TEUCHOS_TEST_FOR_EXCEPTION_PURE_MSG(true,std::logic_error,ss.str());
248 
249  }
250 
251  // We now have a list of sideset entities, lets unwrap them and create the sideset_info!
252  std::map<ParentOrdinal,std::vector<SubcellOrdinal> > owned_parent_cell_to_subcell_indexes;
253  {
254 
255  // This is the subcell dimension we are trying to line up on the sideset
256  const size_t face_subcell_dimension = static_cast<size_t>(mesh.getCellTopology(element_block_name)->getDimension()-1);
257 
258  // List of local subcell indexes indexed by element:
259  // For example: a Tet (element) would have
260  // - 4 triangular faces (subcell_index 0-3, subcell_dimension=2)
261  // - 6 edges (subcell_index 0-5, subcell_dimension=1)
262  // - 4 vertices (subcell_index 0-3, subcell_dimension=0)
263  // Another example: a Line (element) would have
264  // - 2 vertices (subcell_index 0-1, subcell_dimension=0)
265  std::vector<stk::mesh::Entity> elements;
266  std::vector<size_t> subcell_indexes;
267  std::vector<size_t> subcell_dimensions;
268  panzer_stk::workset_utils::getSideElementCascade(mesh, element_block_name, side_entities, subcell_dimensions, subcell_indexes, elements);
269  const size_t num_elements = subcell_dimensions.size();
270 
271  // We need a fast lookup for maping local indexes to parent indexes
272  std::unordered_map<LocalOrdinal,ParentOrdinal> local_to_parent_lookup;
273  auto local_cells_h = Kokkos::create_mirror_view(mesh_info.local_cells);
274  Kokkos::deep_copy(local_cells_h, mesh_info.local_cells);
275  for(ParentOrdinal parent_index=0; parent_index<static_cast<ParentOrdinal>(mesh_info.local_cells.extent(0)); ++parent_index)
276  local_to_parent_lookup[local_cells_h(parent_index)] = parent_index;
277 
278  // Add the subcell indexes for each element on the sideset
279  // TODO: There is a lookup call here to map from local indexes to parent indexes which slows things down. Maybe there is a way around that
280  for(size_t element_index=0; element_index<num_elements; ++element_index) {
281  const size_t subcell_dimension = subcell_dimensions[element_index];
282 
283  // Add subcell to map
284  if(subcell_dimension == face_subcell_dimension){
285  const SubcellOrdinal subcell_index = static_cast<SubcellOrdinal>(subcell_indexes[element_index]);
286  const LocalOrdinal element_local_index = static_cast<LocalOrdinal>(mesh.elementLocalId(elements[element_index]));
287 
288  // Look up the parent cell index using the local cell index
289  const auto itr = local_to_parent_lookup.find(element_local_index);
290  TEUCHOS_ASSERT(itr!= local_to_parent_lookup.end());
291  const ParentOrdinal element_parent_index = itr->second;
292 
293  owned_parent_cell_to_subcell_indexes[element_parent_index].push_back(subcell_index);
294  }
295  }
296  }
297 
298  // We now know the mapping of parent cell indexes to subcell indexes touching the sideset
299 
300  const panzer::LocalOrdinal num_owned_cells = owned_parent_cell_to_subcell_indexes.size();
301 
302  sideset_info.element_block_name = element_block_name;
303  sideset_info.sideset_name = sideset_name;
304  sideset_info.cell_topology = mesh.getCellTopology(element_block_name);
305 
306  sideset_info.num_owned_cells = num_owned_cells;
307 
308  struct face_t{
309  face_t(const ParentOrdinal c0,
310  const ParentOrdinal c1,
311  const SubcellOrdinal sc0,
312  const SubcellOrdinal sc1)
313  {
314  cell_0=c0;
315  cell_1=c1;
316  subcell_index_0=sc0;
317  subcell_index_1=sc1;
318  }
319  ParentOrdinal cell_0;
320  ParentOrdinal cell_1;
321  SubcellOrdinal subcell_index_0;
322  SubcellOrdinal subcell_index_1;
323  };
324 
325 
326  // Figure out how many cells on the other side of the sideset are ghost or virtual
327  std::unordered_set<panzer::LocalOrdinal> owned_parent_cells_set, ghstd_parent_cells_set, virtual_parent_cells_set;
328  std::vector<face_t> faces;
329  {
330  auto cell_to_faces_h = Kokkos::create_mirror_view(mesh_info.cell_to_faces);
331  auto face_to_cells_h = Kokkos::create_mirror_view(mesh_info.face_to_cells);
332  auto face_to_lidx_h = Kokkos::create_mirror_view(mesh_info.face_to_lidx);
333  Kokkos::deep_copy(cell_to_faces_h, mesh_info.cell_to_faces);
334  Kokkos::deep_copy(face_to_cells_h, mesh_info.face_to_cells);
335  Kokkos::deep_copy(face_to_lidx_h, mesh_info.face_to_lidx);
336  panzer::LocalOrdinal parent_virtual_cell_offset = mesh_info.num_owned_cells + mesh_info.num_ghstd_cells;
337  for(const auto & local_cell_index_pair : owned_parent_cell_to_subcell_indexes){
338 
339  const ParentOrdinal parent_cell = local_cell_index_pair.first;
340  const auto & subcell_indexes = local_cell_index_pair.second;
341 
342  owned_parent_cells_set.insert(parent_cell);
343 
344  for(const SubcellOrdinal & subcell_index : subcell_indexes){
345 
346  const LocalOrdinal face = cell_to_faces_h(parent_cell, subcell_index);
347  const LocalOrdinal face_other_side = (face_to_cells_h(face,0) == parent_cell) ? 1 : 0;
348 
349  TEUCHOS_ASSERT(subcell_index == face_to_lidx_h(face, 1-face_other_side));
350 
351  const ParentOrdinal other_side_cell = face_to_cells_h(face, face_other_side);
352  const SubcellOrdinal other_side_subcell_index = face_to_lidx_h(face, face_other_side);
353 
354  faces.push_back(face_t(parent_cell, other_side_cell, subcell_index, other_side_subcell_index));
355 
356  if(other_side_cell >= parent_virtual_cell_offset){
357  virtual_parent_cells_set.insert(other_side_cell);
358  } else {
359  ghstd_parent_cells_set.insert(other_side_cell);
360  }
361  }
362  }
363  }
364 
365  std::vector<ParentOrdinal> all_cells;
366  all_cells.insert(all_cells.end(),owned_parent_cells_set.begin(),owned_parent_cells_set.end());
367  all_cells.insert(all_cells.end(),ghstd_parent_cells_set.begin(),ghstd_parent_cells_set.end());
368  all_cells.insert(all_cells.end(),virtual_parent_cells_set.begin(),virtual_parent_cells_set.end());
369 
370  sideset_info.num_ghstd_cells = ghstd_parent_cells_set.size();
371  sideset_info.num_virtual_cells = virtual_parent_cells_set.size();
372 
373  const LocalOrdinal num_real_cells = sideset_info.num_owned_cells + sideset_info.num_ghstd_cells;
374  const LocalOrdinal num_total_cells = num_real_cells + sideset_info.num_virtual_cells;
375  const LocalOrdinal num_vertices_per_cell = mesh_info.cell_vertices.extent(1);
376  const LocalOrdinal num_dims = mesh_info.cell_vertices.extent(2);
377 
378  // Copy local indexes, global indexes, and cell vertices to sideset info
379  {
380  sideset_info.global_cells = PHX::View<panzer::GlobalOrdinal*>("global_cells", num_total_cells);
381  sideset_info.local_cells = PHX::View<LocalOrdinal*>("local_cells", num_total_cells);
382  sideset_info.cell_vertices = PHX::View<double***>("cell_vertices", num_total_cells, num_vertices_per_cell, num_dims);
383  Kokkos::deep_copy(sideset_info.cell_vertices,0.);
384 
385  typename PHX::View<ParentOrdinal*>::HostMirror all_cells_h("all_cells_h", num_total_cells);
386  PHX::View<ParentOrdinal*> all_cells_d("all_cells_d", num_total_cells);
387  for(LocalOrdinal i=0; i<num_total_cells; ++i)
388  all_cells_h(i) = all_cells[i];
389  Kokkos::deep_copy(all_cells_d, all_cells_h);
390  Kokkos::parallel_for(num_total_cells, KOKKOS_LAMBDA (LocalOrdinal i) {
391  const ParentOrdinal parent_cell = all_cells_d(i);
392  sideset_info.local_cells(i) = mesh_info.local_cells(parent_cell);
393  sideset_info.global_cells(i) = mesh_info.global_cells(parent_cell);
394  for(LocalOrdinal j=0; j<num_vertices_per_cell; ++j)
395  for(LocalOrdinal k=0; k<num_dims; ++k)
396  sideset_info.cell_vertices(i,j,k) = mesh_info.cell_vertices(parent_cell,j,k);
397  });
398  }
399 
400  // Now we have to set the connectivity for the faces.
401 
402  const LocalOrdinal num_faces = faces.size();
403  const LocalOrdinal num_faces_per_cell = mesh_info.cell_to_faces.extent(1);
404 
405  sideset_info.face_to_cells = PHX::View<LocalOrdinal*[2]>("face_to_cells", num_faces);
406  sideset_info.face_to_lidx = PHX::View<LocalOrdinal*[2]>("face_to_lidx", num_faces);
407  sideset_info.cell_to_faces = PHX::View<LocalOrdinal**>("cell_to_faces", num_total_cells, num_faces_per_cell);
408  auto cell_to_faces_h = Kokkos::create_mirror_view(sideset_info.cell_to_faces);
409  auto face_to_cells_h = Kokkos::create_mirror_view(sideset_info.face_to_cells);
410  auto face_to_lidx_h = Kokkos::create_mirror_view(sideset_info.face_to_lidx);
411 
412  // Default the system with invalid cell index - this will be most of the entries
413  Kokkos::deep_copy(cell_to_faces_h, -1);
414 
415  // Lookup for sideset cell index given parent cell index
416  std::unordered_map<ParentOrdinal,ParentOrdinal> parent_to_sideset_lookup;
417  for(unsigned int i=0; i<all_cells.size(); ++i)
418  parent_to_sideset_lookup[all_cells[i]] = i;
419 
420  for(int face_index=0;face_index<num_faces;++face_index){
421  const face_t & face = faces[face_index];
422  const ParentOrdinal & parent_cell_0 = face.cell_0;
423  const ParentOrdinal & parent_cell_1 = face.cell_1;
424 
425  // Convert the parent cell index into a sideset cell index
426  const auto itr_0 = parent_to_sideset_lookup.find(parent_cell_0);
427  TEUCHOS_ASSERT(itr_0 != parent_to_sideset_lookup.end());
428  const ParentOrdinal sideset_cell_0 = itr_0->second;
429 
430  const auto itr_1 = parent_to_sideset_lookup.find(parent_cell_1);
431  TEUCHOS_ASSERT(itr_1 != parent_to_sideset_lookup.end());
432  const ParentOrdinal sideset_cell_1 = itr_1->second;
433 
434 // const ParentOrdinal sideset_cell_0 = std::distance(all_cells.begin(), std::find(all_cells.begin(), all_cells.end(), parent_cell_0));
435 // const ParentOrdinal sideset_cell_1 = std::distance(all_cells.begin(), std::find(all_cells.begin()+num_owned_cells, all_cells.end(), parent_cell_1));
436 
437  face_to_cells_h(face_index,0) = sideset_cell_0;
438  face_to_cells_h(face_index,1) = sideset_cell_1;
439 
440  face_to_lidx_h(face_index,0) = face.subcell_index_0;
441  face_to_lidx_h(face_index,1) = face.subcell_index_1;
442 
443  cell_to_faces_h(sideset_cell_0,face.subcell_index_0) = face_index;
444  cell_to_faces_h(sideset_cell_1,face.subcell_index_1) = face_index;
445 
446  }
447  Kokkos::deep_copy(sideset_info.cell_to_faces, cell_to_faces_h);
448  Kokkos::deep_copy(sideset_info.face_to_cells, face_to_cells_h);
449  Kokkos::deep_copy(sideset_info.face_to_lidx, face_to_lidx_h );
450 
451 }
452 
453 } // namespace
454 
455 Teuchos::RCP<panzer::LocalMeshInfo>
457 {
458  using Teuchos::RCP;
459  using Teuchos::rcp;
460 
461  //typedef Tpetra::CrsMatrix<int,panzer::LocalOrdinal,panzer::GlobalOrdinal> crs_type;
462  typedef Tpetra::Map<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> map_type;
463  typedef Tpetra::Import<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> import_type;
464  //typedef Tpetra::MultiVector<double,panzer::LocalOrdinal,panzer::GlobalOrdinal> mvec_type;
465  //typedef Tpetra::MultiVector<panzer::GlobalOrdinal,panzer::LocalOrdinal,panzer::GlobalOrdinal> ordmvec_type;
466 
467  auto mesh_info_rcp = Teuchos::rcp(new panzer::LocalMeshInfo);
468  auto & mesh_info = *mesh_info_rcp;
469 
470  // Make sure the STK interface is valid
471  TEUCHOS_ASSERT(mesh.isInitialized());
472 
473  // This is required by some of the STK stuff
474  TEUCHOS_ASSERT(typeid(panzer::LocalOrdinal) == typeid(int));
475 
476  Teuchos::RCP<const Teuchos::Comm<int> > comm = mesh.getComm();
477 
478  TEUCHOS_FUNC_TIME_MONITOR_DIFF("panzer_stk::generateLocalMeshInfo",GenerateLocalMeshInfo);
479 
480  // This horrible line of code is required since the connection manager only takes rcps of a mesh
481  RCP<const panzer_stk::STK_Interface> mesh_rcp = Teuchos::rcpFromRef(mesh);
482  // We're allowed to do this since the connection manager only exists in this scope... even though it is also an RCP...
483 
484  // extract topology handle
485  RCP<panzer::ConnManager> conn_rcp = rcp(new panzer_stk::STKConnManager(mesh_rcp));
486  panzer::ConnManager & conn = *conn_rcp;
487 
488  PHX::View<panzer::GlobalOrdinal*> owned_cells, ghost_cells, virtual_cells;
489  panzer::fillLocalCellIDs(comm, conn, owned_cells, ghost_cells, virtual_cells);
490 
491  // build cell maps
493 
494  RCP<map_type> owned_cell_map = rcp(new map_type(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),owned_cells,0,comm));
495  RCP<map_type> ghstd_cell_map = rcp(new map_type(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),ghost_cells,0,comm));
496 
497  // build importer: cell importer, owned to ghstd
498  RCP<import_type> cellimport_own2ghst = rcp(new import_type(owned_cell_map,ghstd_cell_map));
499 
500  // read all the vertices associated with these elements, get ghstd contributions
502 
503  // TODO: This all needs to be rewritten for when element blocks have different cell topologies
504  std::vector<std::string> element_block_names;
505  mesh.getElementBlockNames(element_block_names);
506 
507  const shards::CellTopology & cell_topology = *(mesh.getCellTopology(element_block_names[0]));
508 
509  const int space_dim = cell_topology.getDimension();
510  const int vertices_per_cell = cell_topology.getVertexCount();
511  const int faces_per_cell = cell_topology.getSubcellCount(space_dim-1);
512 
513  // PHX::View<double***> owned_vertices("owned_vertices",localCells.size(),vertices_per_cell,space_dim);
514  Kokkos::DynRankView<double,PHX::Device> owned_vertices("owned_vertices",owned_cells.extent(0),vertices_per_cell,space_dim);
515  {
516  std::vector<std::size_t> localCells(owned_cells.extent(0),Teuchos::OrdinalTraits<std::size_t>::invalid());
517  for(size_t i=0;i<localCells.size();i++)
518  localCells[i] = i;
519  mesh.getElementVerticesNoResize(localCells,owned_vertices);
520  }
521 
522  // this builds a ghstd vertex array
523  Kokkos::DynRankView<double,PHX::Device> ghstd_vertices = buildGhostedVertices(*cellimport_own2ghst,owned_vertices);
524 
525  // build edge to cell neighbor mapping
527 
528  auto owned_cells_h = Kokkos::create_mirror_view(owned_cells);
529  auto ghost_cells_h = Kokkos::create_mirror_view(ghost_cells);
530  Kokkos::deep_copy(owned_cells_h, owned_cells);
531  Kokkos::deep_copy(ghost_cells_h, ghost_cells);
532  std::unordered_map<panzer::GlobalOrdinal,int> global_to_local;
533  global_to_local[-1] = -1; // this is the "no neighbor" flag
534  for(size_t i=0;i<owned_cells.extent(0);i++)
535  global_to_local[owned_cells_h(i)] = i;
536  for(size_t i=0;i<ghost_cells.extent(0);i++)
537  global_to_local[ghost_cells_h(i)] = i+Teuchos::as<int>(owned_cells.extent(0));
538 
539  // this class comes from Mini-PIC and Matt B
540  RCP<panzer::FaceToElement<panzer::LocalOrdinal,panzer::GlobalOrdinal> > faceToElement = rcp(new panzer::FaceToElement<panzer::LocalOrdinal,panzer::GlobalOrdinal>());
541  faceToElement->initialize(conn);
542  auto elems_by_face = faceToElement->getFaceToElementsMap();
543  auto face_to_lidx = faceToElement->getFaceToCellLocalIdxMap();
544 
545  // We also need to consider faces that connect to cells that do not exist, but are needed for boundary conditions
546  // We dub them virtual cell since there should be no geometry associated with them, or topology really
547  // They exist only for datastorage so that they are consistent with 'real' cells from an algorithm perspective
548 
549  // Each virtual face (face linked to a '-1' cell) requires a virtual cell (i.e. turn the '-1' into a virtual cell)
550  // Virtual cells are those that do not exist but are connected to an owned cell
551  // 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.
552 
553  // Iterate over all faces and identify the faces connected to a potential virtual cell
554 
555  const panzer::LocalOrdinal num_owned_cells = owned_cells.extent(0);
556  const panzer::LocalOrdinal num_ghstd_cells = ghost_cells.extent(0);
557  const panzer::LocalOrdinal num_virtual_cells = virtual_cells.extent(0);
558 
559  // total cells and faces include owned, ghosted, and virtual
560  const panzer::LocalOrdinal num_real_cells = num_owned_cells + num_ghstd_cells;
561  const panzer::LocalOrdinal num_total_cells = num_real_cells + num_virtual_cells;
562  const panzer::LocalOrdinal num_total_faces = elems_by_face.extent(0);
563 
564  // Lookup cells connected to a face
565  PHX::View<panzer::LocalOrdinal*[2]> face_to_cells("face_to_cells",num_total_faces);
566 
567  // Lookup local face indexes given cell and left/right state (0/1)
568  PHX::View<panzer::LocalOrdinal*[2]> face_to_localidx("face_to_localidx",num_total_faces);
569 
570  // Lookup face index given a cell and local face index
571  PHX::View<panzer::LocalOrdinal**> cell_to_face("cell_to_face",num_total_cells,faces_per_cell);
572 
573  // initialize with negative one cells that are not associated with a face
574  Kokkos::deep_copy(cell_to_face,-1);
575 
576  // Transfer information from 'faceToElement' datasets to local arrays
577  {
578  PANZER_FUNC_TIME_MONITOR_DIFF("Transer faceToElement to local",TransferFaceToElementLocal);
579 
580  int virtual_cell_index = num_real_cells;
581  auto elems_by_face_h = Kokkos::create_mirror_view(elems_by_face);
582  auto face_to_lidx_h = Kokkos::create_mirror_view(face_to_lidx);
583  auto face_to_cells_h = Kokkos::create_mirror_view(face_to_cells);
584  auto face_to_localidx_h = Kokkos::create_mirror_view(face_to_localidx);
585  auto cell_to_face_h = Kokkos::create_mirror_view(cell_to_face);
586  Kokkos::deep_copy(elems_by_face_h, elems_by_face);
587  Kokkos::deep_copy(face_to_lidx_h, face_to_lidx);
588  for(size_t f=0;f<elems_by_face.extent(0);f++) {
589 
590  const panzer::GlobalOrdinal global_c0 = elems_by_face_h(f,0);
591  const panzer::GlobalOrdinal global_c1 = elems_by_face_h(f,1);
592 
593  // make sure that no bonus cells get in here
594  TEUCHOS_ASSERT(global_to_local.find(global_c0)!=global_to_local.end());
595  TEUCHOS_ASSERT(global_to_local.find(global_c1)!=global_to_local.end());
596 
597  auto c0 = global_to_local[global_c0];
598  auto lidx0 = face_to_lidx_h(f,0);
599  auto c1 = global_to_local[global_c1];
600  auto lidx1 = face_to_lidx_h(f,1);
601 
602  // Test for virtual cells
603 
604  // Left cell
605  if(c0 < 0){
606  // Virtual cell - create it!
607  c0 = virtual_cell_index++;
608 
609  // We need the subcell_index to line up between real and virtual cell
610  // This way the face has the same geometry... though the face normal
611  // will point in the wrong direction
612  lidx0 = lidx1;
613  }
614  cell_to_face_h(c0,lidx0) = f;
615 
616 
617  // Right cell
618  if(c1 < 0){
619  // Virtual cell - create it!
620  c1 = virtual_cell_index++;
621 
622  // We need the subcell_index to line up between real and virtual cell
623  // This way the face has the same geometry... though the face normal
624  // will point in the wrong direction
625  lidx1 = lidx0;
626  }
627  cell_to_face_h(c1,lidx1) = f;
628 
629  // Faces point from low cell index to high cell index
630  if(c0<c1){
631  face_to_cells_h(f,0) = c0;
632  face_to_localidx_h(f,0) = lidx0;
633  face_to_cells_h(f,1) = c1;
634  face_to_localidx_h(f,1) = lidx1;
635  } else {
636  face_to_cells_h(f,0) = c1;
637  face_to_localidx_h(f,0) = lidx1;
638  face_to_cells_h(f,1) = c0;
639  face_to_localidx_h(f,1) = lidx0;
640  }
641 
642  // We should avoid having two virtual cells linked together.
643  TEUCHOS_ASSERT(c0<num_real_cells or c1<num_real_cells)
644 
645  }
646  Kokkos::deep_copy(face_to_cells, face_to_cells_h);
647  Kokkos::deep_copy(face_to_localidx, face_to_localidx_h);
648  Kokkos::deep_copy(cell_to_face, cell_to_face_h);
649  }
650  // at this point all the data structures have been built, so now we can "do" DG.
651  // There are two of everything, an "owned" data structured corresponding to "owned"
652  // cells. And a "ghstd" data structure corresponding to ghosted cells
654  {
655  PANZER_FUNC_TIME_MONITOR_DIFF("Assign Indices",AssignIndices);
656  mesh_info.cell_to_faces = cell_to_face;
657  mesh_info.face_to_cells = face_to_cells; // faces
658  mesh_info.face_to_lidx = face_to_localidx;
659  mesh_info.subcell_dimension = space_dim;
660  mesh_info.subcell_index = -1;
661  mesh_info.has_connectivity = true;
662 
663  mesh_info.num_owned_cells = owned_cells.extent(0);
664  mesh_info.num_ghstd_cells = ghost_cells.extent(0);
665  mesh_info.num_virtual_cells = virtual_cells.extent(0);
666 
667  mesh_info.global_cells = PHX::View<panzer::GlobalOrdinal*>("global_cell_indices",num_total_cells);
668  mesh_info.local_cells = PHX::View<panzer::LocalOrdinal*>("local_cell_indices",num_total_cells);
669 
670  Kokkos::parallel_for(num_owned_cells,KOKKOS_LAMBDA (int i) {
671  mesh_info.global_cells(i) = owned_cells(i);
672  mesh_info.local_cells(i) = i;
673  });
674 
675  Kokkos::parallel_for(num_ghstd_cells,KOKKOS_LAMBDA (int i) {
676  mesh_info.global_cells(i+num_owned_cells) = ghost_cells(i);
677  mesh_info.local_cells(i+num_owned_cells) = i+num_owned_cells;
678  });
679 
680  Kokkos::parallel_for(num_virtual_cells,KOKKOS_LAMBDA (int i) {
681  mesh_info.global_cells(i+num_real_cells) = virtual_cells(i);
682  mesh_info.local_cells(i+num_real_cells) = i+num_real_cells;
683  });
684 
685  mesh_info.cell_vertices = PHX::View<double***>("cell_vertices",num_total_cells,vertices_per_cell,space_dim);
686 
687  // Initialize coordinates to zero
688  Kokkos::deep_copy(mesh_info.cell_vertices, 0.);
689 
690  Kokkos::parallel_for(num_owned_cells,KOKKOS_LAMBDA (int i) {
691  for(int j=0;j<vertices_per_cell;++j){
692  for(int k=0;k<space_dim;++k){
693  mesh_info.cell_vertices(i,j,k) = owned_vertices(i,j,k);
694  }
695  }
696  });
697 
698  Kokkos::parallel_for(num_ghstd_cells,KOKKOS_LAMBDA (int i) {
699  for(int j=0;j<vertices_per_cell;++j){
700  for(int k=0;k<space_dim;++k){
701  mesh_info.cell_vertices(i+num_owned_cells,j,k) = ghstd_vertices(i,j,k);
702  }
703  }
704  });
705 
706  // This will backfire at some point, but we're going to make the virtual cell have the same geometry as the cell it interfaces with
707  // This way we can define a virtual cell geometry without extruding the face outside of the domain
708  {
709  PANZER_FUNC_TIME_MONITOR_DIFF("Assign geometry traits",AssignGeometryTraits);
710  Kokkos::parallel_for(num_virtual_cells,KOKKOS_LAMBDA (int i) {
711  const panzer::LocalOrdinal virtual_cell = i+num_real_cells;
712  for(int local_face=0; local_face<faces_per_cell; ++local_face){
713  const panzer::LocalOrdinal face = cell_to_face(virtual_cell, local_face);
714  if(face >= 0){
715  const panzer::LocalOrdinal other_side = (face_to_cells(face, 0) == virtual_cell) ? 1 : 0;
716  const panzer::LocalOrdinal real_cell = face_to_cells(face,other_side);
717  for(int j=0;j<vertices_per_cell;++j){
718  for(int k=0;k<space_dim;++k){
719  mesh_info.cell_vertices(virtual_cell,j,k) = mesh_info.cell_vertices(real_cell,j,k);
720  }
721  }
722  break;
723  }
724  }
725  });
726 
727  }
728  }
729 
730  // Setup element blocks and sidesets
731  std::vector<std::string> sideset_names;
732  mesh.getSidesetNames(sideset_names);
733 
734  for(const std::string & element_block_name : element_block_names){
735  PANZER_FUNC_TIME_MONITOR_DIFF("Set up setupLocalMeshBlockInfo",SetupLocalMeshBlockInfo);
736  panzer::LocalMeshBlockInfo & block_info = mesh_info.element_blocks[element_block_name];
737  setupLocalMeshBlockInfo(mesh, conn, mesh_info, element_block_name, block_info);
738  block_info.subcell_dimension = space_dim;
739  block_info.subcell_index = -1;
740  block_info.has_connectivity = true;
741 
742  // Setup sidesets
743  for(const std::string & sideset_name : sideset_names){
744  PANZER_FUNC_TIME_MONITOR_DIFF("Setup LocalMeshSidesetInfo",SetupLocalMeshSidesetInfo);
745  panzer::LocalMeshSidesetInfo & sideset_info = mesh_info.sidesets[element_block_name][sideset_name];
746  setupLocalMeshSidesetInfo(mesh, conn, mesh_info, element_block_name, sideset_name, sideset_info);
747  sideset_info.subcell_dimension = space_dim;
748  sideset_info.subcell_index = -1;
749  sideset_info.has_connectivity = true;
750  }
751 
752  }
753 
754  return mesh_info_rcp;
755 
756 }
757 
758 }
std::map< std::string, std::map< std::string, LocalMeshSidesetInfo > > sidesets
std::size_t elementLocalId(stk::mesh::Entity elmt) const
Teuchos::RCP< panzer::LocalMeshInfo > generateLocalMeshInfo(const panzer_stk::STK_Interface &mesh)
Create a structure containing information about the local portion of a given element block...
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
PHX::View< panzer::LocalOrdinal *[2]> face_to_cells
panzer::LocalOrdinal num_owned_cells
void getElementBlockNames(std::vector< std::string > &names) const
panzer::LocalOrdinal num_virtual_cells
Teuchos::RCP< const shards::CellTopology > cell_topology
PHX::View< panzer::LocalOrdinal *[2]> face_to_lidx
Pure virtual base class for supplying mesh connectivity information to the DOF Manager.
bool isInitialized() const
Has initialize been called on this mesh object?
void getElementVerticesNoResize(const std::vector< std::size_t > &localIds, ArrayT &vertices) const
panzer::LocalOrdinal num_ghstd_cells
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
get the comm associated with this mesh
std::map< std::string, LocalMeshBlockInfo > element_blocks
Teuchos::RCP< const shards::CellTopology > getCellTopology(const std::string &eBlock) const
void getAllSides(const std::string &sideName, std::vector< stk::mesh::Entity > &sides) const
void getSidesetNames(std::vector< std::string > &name) const
void getSideElementCascade(const panzer_stk::STK_Interface &mesh, const std::string &blockId, const std::vector< stk::mesh::Entity > &sides, std::vector< std::size_t > &localSubcellDim, std::vector< std::size_t > &localSubcellIds, std::vector< stk::mesh::Entity > &elements)
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
virtual std::string getBlockId(LocalOrdinal localElmtId) const =0
PHX::View< panzer::LocalOrdinal ** > cell_to_faces