Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_STKConnManager.cpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Panzer: A partial differential equation assembly
4 // engine for strongly coupled complex multiphysics systems
5 //
6 // Copyright 2011 NTESS and the Panzer contributors.
7 // SPDX-License-Identifier: BSD-3-Clause
8 // *****************************************************************************
9 // @HEADER
10 
12 
13 #include <vector>
14 
15 // Teuchos includes
16 #include "Teuchos_RCP.hpp"
17 
18 #include "Kokkos_DynRankView.hpp"
19 
24 
25 #include "Teuchos_FancyOStream.hpp"
26 
27 namespace panzer_stk {
28 
29 using Teuchos::RCP;
30 using Teuchos::rcp;
31 
32 // Initialize statics
34 std::vector<STKConnManager::CachedEntry> STKConnManager::cached_conn_managers_{};
36 
37 // Object describing how to sort a vector of elements using
38 // local ID as the key
39 class LocalIdCompare {
40 public:
41  LocalIdCompare(const RCP<const STK_Interface> & mesh) : mesh_(mesh) {}
42 
43  // Compares two stk mesh entities based on local ID
44  bool operator() (stk::mesh::Entity a,stk::mesh::Entity b)
45  { return mesh_->elementLocalId(a) < mesh_->elementLocalId(b);}
46 
47 private:
48  RCP<const STK_Interface> mesh_;
49 };
50 
52  : stkMeshDB_(stkMeshDB), ownedElementCount_(0)
53 {
54 }
55 
57  : stkMeshDB_(cm.stkMeshDB_),
58  elements_(cm.elements_),
59  elementBlocks_(cm.elementBlocks_),
60  neighborElementBlocks_(cm.neighborElementBlocks_),
61  blockIdToIndex_(cm.blockIdToIndex_),
62  elmtLidToConn_(cm.elmtLidToConn_),
63  connSize_(cm.connSize_),
64  connectivity_(cm.connectivity_),
65  ownedElementCount_(cm.ownedElementCount_),
66  sidesetsToAssociate_(cm.sidesetsToAssociate_),
67  sidesetYieldedAssociations_(cm.sidesetYieldedAssociations_)
68 {
70  for (size_t i=0; i < elmtToAssociatedElmts_.size(); ++i)
72 }
73 
76 {
78 }
79 
81 {
82  elements_ = Teuchos::null;
83 
84  elementBlocks_.clear();
85  elmtLidToConn_.clear();
86  connSize_.clear();
87  elmtToAssociatedElmts_.clear();
88 }
89 
91 {
92  clearLocalElementMapping(); // forget the past
93 
94  // build element block information
96  elements_ = Teuchos::rcp(new std::vector<stk::mesh::Entity>);
97 
98  // defines ordering of blocks
99  std::vector<std::string> blockIds;
100  stkMeshDB_->getElementBlockNames(blockIds);
101 
102  std::size_t blockIndex=0;
103  for(std::vector<std::string>::const_iterator idItr=blockIds.begin();
104  idItr!=blockIds.end();++idItr,++blockIndex) {
105  std::string blockId = *idItr;
106 
107  // grab elements on this block
108  std::vector<stk::mesh::Entity> blockElmts;
109  stkMeshDB_->getMyElements(blockId,blockElmts);
110 
111  // concatenate them into element LID lookup table
112  elements_->insert(elements_->end(),blockElmts.begin(),blockElmts.end());
113 
114  // build block to LID map
115  elementBlocks_[blockId] = Teuchos::rcp(new std::vector<LocalOrdinal>);
116  for(std::size_t i=0;i<blockElmts.size();i++)
117  elementBlocks_[blockId]->push_back(stkMeshDB_->elementLocalId(blockElmts[i]));
118  }
119 
120  ownedElementCount_ = elements_->size();
121 
122  blockIndex=0;
123  for(std::vector<std::string>::const_iterator idItr=blockIds.begin();
124  idItr!=blockIds.end();++idItr,++blockIndex) {
125  std::string blockId = *idItr;
126 
127  // grab elements on this block
128  std::vector<stk::mesh::Entity> blockElmts;
129  stkMeshDB_->getNeighborElements(blockId,blockElmts);
130 
131  // concatenate them into element LID lookup table
132  elements_->insert(elements_->end(),blockElmts.begin(),blockElmts.end());
133 
134  // build block to LID map
135  neighborElementBlocks_[blockId] = Teuchos::rcp(new std::vector<LocalOrdinal>);
136  for(std::size_t i=0;i<blockElmts.size();i++)
137  neighborElementBlocks_[blockId]->push_back(stkMeshDB_->elementLocalId(blockElmts[i]));
138  }
139 
140  // this expensive operation gurantees ordering of local IDs
141  std::sort(elements_->begin(),elements_->end(),LocalIdCompare(stkMeshDB_));
142 
143  // allocate space for element LID to Connectivty map
144  // connectivity size
145  elmtLidToConn_.clear();
146  elmtLidToConn_.resize(elements_->size(),0);
147 
148  connSize_.clear();
149  connSize_.resize(elements_->size(),0);
150 }
151 
152 void
154  LocalOrdinal & nodeIdCnt, LocalOrdinal & edgeIdCnt,
155  LocalOrdinal & faceIdCnt, LocalOrdinal & cellIdCnt,
156  GlobalOrdinal & nodeOffset, GlobalOrdinal & edgeOffset,
157  GlobalOrdinal & faceOffset, GlobalOrdinal & cellOffset) const
158 {
159  // get the global counts for all the nodes, faces, edges and cells
160  GlobalOrdinal maxNodeId = stkMeshDB_->getMaxEntityId(stkMeshDB_->getNodeRank());
161  GlobalOrdinal maxEdgeId = stkMeshDB_->getMaxEntityId(stkMeshDB_->getEdgeRank());
162  GlobalOrdinal maxFaceId = stkMeshDB_->getMaxEntityId(stkMeshDB_->getFaceRank());
163 
164  // compute ID counts for each sub cell type
165  int patternDim = fp.getDimension();
166  switch(patternDim) {
167  case 3:
168  faceIdCnt = fp.getSubcellIndices(2,0).size();
169  // Intentional fall-through.
170  case 2:
171  edgeIdCnt = fp.getSubcellIndices(1,0).size();
172  // Intentional fall-through.
173  case 1:
174  nodeIdCnt = fp.getSubcellIndices(0,0).size();
175  cellIdCnt = fp.getSubcellIndices(patternDim,0).size();
176  break;
177  case 0:
178  default:
179  TEUCHOS_ASSERT(false);
180  };
181 
182  // compute offsets for each sub cell type
183  nodeOffset = 0;
184  edgeOffset = nodeOffset+(maxNodeId+1)*nodeIdCnt;
185  faceOffset = edgeOffset+(maxEdgeId+1)*edgeIdCnt;
186  cellOffset = faceOffset+(maxFaceId+1)*faceIdCnt;
187 
188  // sanity check
189  TEUCHOS_ASSERT(nodeOffset <= edgeOffset
190  && edgeOffset <= faceOffset
191  && faceOffset <= cellOffset);
192 }
193 
196  unsigned subcellRank,
197  LocalOrdinal idCnt,
198  GlobalOrdinal offset,
199  const unsigned maxIds)
200 {
201  if(idCnt<=0)
202  return 0 ;
203 
204  // loop over all relations of specified type, unless requested
205  LocalOrdinal numIds = 0;
206  stk::mesh::BulkData& bulkData = *stkMeshDB_->getBulkData();
207  const stk::mesh::EntityRank rank = static_cast<stk::mesh::EntityRank>(subcellRank);
208 
209 #ifdef PANZER_DEBUG
210  TEUCHOS_TEST_FOR_EXCEPTION(maxIds > bulkData.num_connectivity(element, rank),
211  std::runtime_error,
212  "ERROR: The maxIds (num vertices from basis cell topology) is greater than the num vertices in the stk mesh topology!");
213 #endif
214 
215  const size_t num_rels = (maxIds > 0) ? maxIds : bulkData.num_connectivity(element, rank);
216  stk::mesh::Entity const* relations = bulkData.begin(element, rank);
217  for(std::size_t sc=0; sc<num_rels; ++sc) {
218  stk::mesh::Entity subcell = relations[sc];
219 
220  // add connectivities: adjust for STK indexing craziness
221  for(LocalOrdinal i=0;i<idCnt;i++) {
222  connectivity_.push_back(offset+idCnt*(bulkData.identifier(subcell)-1)+i);
223  }
224 
225  numIds += idCnt;
226  }
227  return numIds;
228 }
229 
230 void
232  unsigned subcellRank,unsigned subcellId,GlobalOrdinal newId,
233  GlobalOrdinal offset)
234 {
235  LocalOrdinal elmtLID = stkMeshDB_->elementLocalId(element);
236  auto * conn = this->getConnectivity(elmtLID);
237  const std::vector<int> & subCellIndices = fp.getSubcellIndices(subcellRank,subcellId);
238 
239  // add connectivities: adjust for STK indexing craziness
240  for(std::size_t i=0;i<subCellIndices.size();i++) {
241  conn[subCellIndices[i]] = offset+subCellIndices.size()*(newId-1)+i;
242  }
243 }
244 
246 {
247  PANZER_FUNC_TIME_MONITOR_DIFF("panzer_stk::STKConnManager::buildConnectivity", build_connectivity);
248 
249  if (cache_connectivity_) {
250  auto fp_rcp = fp.clone();
251  auto search = std::find_if(cached_conn_managers_.begin(),
252  cached_conn_managers_.end(),
253  FieldPatternCompare(fp_rcp));
254  if (search != cached_conn_managers_.end()) {
255  PANZER_FUNC_TIME_MONITOR_DIFF("panzer_stk::STKConnManager::copyingCachedConnectivity", copy_cached_connectivity);
256  {
257  STKConnManager& cm = *(search->second);
258  elements_ = cm.elements_;
263  connSize_ = cm.connSize_;
268  {
270  for (size_t i=0; i < elmtToAssociatedElmts_.size(); ++i)
272  }
274  }
275  return;
276  }
277  }
278 
279  stk::mesh::BulkData& bulkData = *stkMeshDB_->getBulkData();
280 
281  // get element info from STK_Interface
282  // object and build a local element mapping.
284 
285  // Build sub cell ID counts and offsets
286  // ID counts = How many IDs belong on each subcell (number of mesh DOF used)
287  // Offset = What is starting index for subcell ID type?
288  // Global numbering goes like [node ids, edge ids, face ids, cell ids]
289  LocalOrdinal nodeIdCnt=0, edgeIdCnt=0, faceIdCnt=0, cellIdCnt=0;
290  GlobalOrdinal nodeOffset=0, edgeOffset=0, faceOffset=0, cellOffset=0;
291  buildOffsetsAndIdCounts(fp, nodeIdCnt, edgeIdCnt, faceIdCnt, cellIdCnt,
292  nodeOffset, edgeOffset, faceOffset, cellOffset);
293 
294  // std::cout << "node: count = " << nodeIdCnt << ", offset = " << nodeOffset << std::endl;
295  // std::cout << "edge: count = " << edgeIdCnt << ", offset = " << edgeOffset << std::endl;
296  // std::cout << "face: count = " << faceIdCnt << ", offset = " << faceOffset << std::endl;
297  // std::cout << "cell: count = " << cellIdCnt << ", offset = " << cellOffset << std::endl;
298 
299  // Connectivity only requires lowest order mesh node information
300  // With the notion of extended topologies used by shards, it is
301  // sufficient to take the first num_vertices nodes for connectivity purposes.
302  const unsigned num_vertices = fp.getCellTopology().getVertexCount();
303 
304  // loop over elements and build global connectivity
305  for(std::size_t elmtLid=0;elmtLid!=elements_->size();++elmtLid) {
306  GlobalOrdinal numIds = 0;
307  stk::mesh::Entity element = (*elements_)[elmtLid];
308 
309  // get index into connectivity array
310  elmtLidToConn_[elmtLid] = connectivity_.size();
311 
312  // add connectivities for sub cells
313  // Second order or higher mesh nodes are only needed downstream by the FE bases
314  // The connection manager only expects first order nodes (vertices), so we subselect if necessary
315  // All edge and face IDs are stored
316  numIds += addSubcellConnectivities(element,stkMeshDB_->getNodeRank(),nodeIdCnt,nodeOffset,num_vertices);
317  numIds += addSubcellConnectivities(element,stkMeshDB_->getEdgeRank(),edgeIdCnt,edgeOffset);
318  numIds += addSubcellConnectivities(element,stkMeshDB_->getFaceRank(),faceIdCnt,faceOffset);
319 
320  // add connectivity for parent cells
321  if(cellIdCnt>0) {
322  // add connectivities: adjust for STK indexing craziness
323  for(LocalOrdinal i=0;i<cellIdCnt;i++)
324  connectivity_.push_back(cellOffset+cellIdCnt*(bulkData.identifier(element)-1));
325 
326  numIds += cellIdCnt;
327  }
328 
329  connSize_[elmtLid] = numIds;
330  }
331 
332  applyPeriodicBCs( fp, nodeOffset, edgeOffset, faceOffset, cellOffset);
333 
334  // This method does not modify connectivity_. But it should be called here
335  // because the data it initializes should be available at the same time as
336  // connectivity_.
339 
340  if (cache_connectivity_) {
341  auto fp_rcp = fp.clone();
342  cached_conn_managers_.push_back(std::make_pair(fp_rcp,Teuchos::rcp(new panzer_stk::STKConnManager(*this))));
343  }
344 }
345 
347 {
348  // walk through the element blocks and figure out which this ID belongs to
349  stk::mesh::Entity element = (*elements_)[localElmtId];
350 
351  return stkMeshDB_->containingBlockId(element);
352 }
353 
355  GlobalOrdinal faceOffset, GlobalOrdinal /* cellOffset */)
356 {
357  using Teuchos::RCP;
358  using Teuchos::rcp;
359 
360  PANZER_FUNC_TIME_MONITOR_DIFF("panzer_stk::STKConnManager::applyPeriodicBCs", apply_periodic_bcs);
361 
362  std::pair<Teuchos::RCP<std::vector<std::pair<std::size_t,std::size_t> > >, Teuchos::RCP<std::vector<unsigned int> > > matchedValues
363  = stkMeshDB_->getPeriodicNodePairing();
364 
366  = matchedValues.first;
368  = matchedValues.second;
369 
370  // no matchedNodes means nothing to do!
371  if(matchedNodes==Teuchos::null) return;
372 
373  for(std::size_t m=0;m<matchedNodes->size();m++) {
374  stk::mesh::EntityId oldNodeId = (*matchedNodes)[m].first;
375  std::size_t newNodeId = (*matchedNodes)[m].second;
376 
377  std::vector<stk::mesh::Entity> elements;
378  std::vector<int> localIds;
379 
380  GlobalOrdinal offset0 = 0; // to make numbering consistent with that in PeriodicBC_Matcher
381  GlobalOrdinal offset1 = 0; // offset for dof indexing
382  if((*matchTypes)[m] == 0)
383  offset1 = nodeOffset-offset0;
384  else if((*matchTypes)[m] == 1){
385  offset0 = stkMeshDB_->getMaxEntityId(stkMeshDB_->getNodeRank());
386  offset1 = edgeOffset-offset0;
387  } else if((*matchTypes)[m] == 2){
388  offset0 = stkMeshDB_->getMaxEntityId(stkMeshDB_->getNodeRank())+stkMeshDB_->getMaxEntityId(stkMeshDB_->getEdgeRank());
389  offset1 = faceOffset-offset0;
390  } else
391  TEUCHOS_ASSERT(false);
392 
393  // get relevent elements and node IDs
394  stkMeshDB_->getOwnedElementsSharingNode(stk::mesh::EntityId(oldNodeId-offset0),elements,localIds,(*matchTypes)[m]);
395 
396  // modify global numbering already built for each element
397  for(std::size_t e=0;e<elements.size();e++){
398  modifySubcellConnectivities(fp,elements[e],(*matchTypes)[m],localIds[e],newNodeId,offset1);
399  }
400 
401  }
402 }
403 
406 void STKConnManager::getDofCoords(const std::string & blockId,
407  const panzer::Intrepid2FieldPattern & coordProvider,
408  std::vector<std::size_t> & localCellIds,
409  Kokkos::DynRankView<double,PHX::Device> & points) const
410 {
411  int dim = coordProvider.getDimension();
412  int numIds = coordProvider.numberIds();
413 
414  // grab element nodes
415  Kokkos::DynRankView<double,PHX::Device> nodes;
416  workset_utils::getIdsAndNodes(*stkMeshDB_,blockId,localCellIds,nodes);
417 
418  // setup output array
419  points = Kokkos::DynRankView<double,PHX::Device>("points",localCellIds.size(),numIds,dim);
420  coordProvider.getInterpolatoryCoordinates(nodes,points,stkMeshDB_->getCellTopology(blockId));
421 }
422 
424 {
425  return ! sidesetsToAssociate_.empty();
426 }
427 
428 void STKConnManager::associateElementsInSideset(const std::string sideset_id)
429 {
430  sidesetsToAssociate_.push_back(sideset_id);
431  sidesetYieldedAssociations_.push_back(false);
432 }
433 
434 inline std::size_t
435 getElementIdx(const std::vector<stk::mesh::Entity>& elements,
436  stk::mesh::Entity const e)
437 {
438  return static_cast<std::size_t>(
439  std::distance(elements.begin(), std::find(elements.begin(), elements.end(), e)));
440 }
441 
443 {
444  stk::mesh::BulkData& bulkData = *stkMeshDB_->getBulkData();
445  elmtToAssociatedElmts_.resize(elements_->size());
446  for (std::size_t i = 0; i < sidesetsToAssociate_.size(); ++i) {
447  std::vector<stk::mesh::Entity> sides;
448  stkMeshDB_->getAllSides(sidesetsToAssociate_[i], sides);
449  sidesetYieldedAssociations_[i] = ! sides.empty();
450  for (std::vector<stk::mesh::Entity>::const_iterator si = sides.begin();
451  si != sides.end(); ++si) {
452  stk::mesh::Entity side = *si;
453  const size_t num_elements = bulkData.num_elements(side);
454  stk::mesh::Entity const* elements = bulkData.begin_elements(side);
455  if (num_elements != 2) {
456  // If relations.size() != 2 for one side in the sideset, then it's true
457  // for all, including the first.
458  TEUCHOS_ASSERT(si == sides.begin());
459  sidesetYieldedAssociations_[i] = false;
460  break;
461  }
462  const std::size_t ea_id = getElementIdx(*elements_, elements[0]),
463  eb_id = getElementIdx(*elements_, elements[1]);
464  elmtToAssociatedElmts_[ea_id].push_back(eb_id);
465  elmtToAssociatedElmts_[eb_id].push_back(ea_id);
466  }
467  }
468 }
469 
470 std::vector<std::string> STKConnManager::
472 {
473  std::vector<std::string> sidesets;
474  for (std::size_t i = 0; i < sidesetYieldedAssociations_.size(); ++i) {
475  int sya, my_sya = sidesetYieldedAssociations_[i] ? 1 : 0;
476  Teuchos::reduceAll(comm, Teuchos::REDUCE_MAX, 1, &my_sya, &sya);
477  if (sya == 0)
478  sidesets.push_back(sidesetsToAssociate_[i]);
479  }
480  return sidesets;
481 }
482 
483 const std::vector<STKConnManager::LocalOrdinal>&
485 {
486  return elmtToAssociatedElmts_[el];
487 }
488 
489 }
virtual Teuchos::RCP< panzer::FieldPattern > clone() const =0
virtual int getDimension() const =0
std::map< std::string, Teuchos::RCP< std::vector< LocalOrdinal > > > elementBlocks_
virtual const std::vector< LocalOrdinal > & getAssociatedNeighbors(const LocalOrdinal &el) const
virtual std::string getBlockId(LocalOrdinal localElmtId) const
panzer::ConnManager::GlobalOrdinal GlobalOrdinal
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
std::vector< bool > sidesetYieldedAssociations_
std::map< std::string, Teuchos::RCP< std::vector< LocalOrdinal > > > neighborElementBlocks_
Teuchos::RCP< const STK_Interface > stkMeshDB_
std::map< std::string, GlobalOrdinal > blockIdToIndex_
std::vector< GlobalOrdinal > connectivity_
void getIdsAndNodes(const panzer_stk::STK_Interface &mesh, std::string blockId, std::vector< std::size_t > &localIds, ArrayT &nodes)
virtual bool hasAssociatedNeighbors() const
std::vector< std::string > checkAssociateElementsInSidesets(const Teuchos::Comm< int > &comm) const
std::vector< std::vector< LocalOrdinal > > elmtToAssociatedElmts_
virtual int numberIds() const
virtual const panzer::GlobalOrdinal * getConnectivity(LocalOrdinal localElmtId) const
void buildOffsetsAndIdCounts(const panzer::FieldPattern &fp, LocalOrdinal &nodeIdCnt, LocalOrdinal &edgeIdCnt, LocalOrdinal &faceIdCnt, LocalOrdinal &cellIdCnt, GlobalOrdinal &nodeOffset, GlobalOrdinal &edgeOffset, GlobalOrdinal &faceOffset, GlobalOrdinal &cellOffset) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual Teuchos::RCP< panzer::ConnManager > noConnectivityClone() const
static std::vector< CachedEntry > cached_conn_managers_
void modifySubcellConnectivities(const panzer::FieldPattern &fp, stk::mesh::Entity element, unsigned subcellRank, unsigned subcellId, GlobalOrdinal newId, GlobalOrdinal offset)
void applyPeriodicBCs(const panzer::FieldPattern &fp, GlobalOrdinal nodeOffset, GlobalOrdinal edgeOffset, GlobalOrdinal faceOffset, GlobalOrdinal cellOffset)
panzer::ConnManager::LocalOrdinal LocalOrdinal
LocalOrdinal addSubcellConnectivities(stk::mesh::Entity element, unsigned subcellRank, LocalOrdinal idCnt, GlobalOrdinal offset, const unsigned maxIds=0)
Loops over relations of a given rank for a specified element and adds a unique ID to the connectivity...
virtual void getDofCoords(const std::string &blockId, const panzer::Intrepid2FieldPattern &coordProvider, std::vector< std::size_t > &localCellIds, Kokkos::DynRankView< double, PHX::Device > &points) const
std::vector< LocalOrdinal > elmtLidToConn_
std::vector< std::string > sidesetsToAssociate_
virtual shards::CellTopology getCellTopology() const =0
void getInterpolatoryCoordinates(Kokkos::DynRankView< double, PHX::Device > &coords) const
virtual const std::vector< int > & getSubcellIndices(int dim, int cellIndex) const =0
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< std::vector< stk::mesh::Entity > > elements_
std::size_t getElementIdx(const std::vector< stk::mesh::Entity > &elements, stk::mesh::Entity const e)
virtual void buildConnectivity(const panzer::FieldPattern &fp)
STKConnManager(const Teuchos::RCP< const STK_Interface > &stkMeshDB)
void associateElementsInSideset(const std::string sideset_id)
std::vector< LocalOrdinal > connSize_