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 //
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 <vector>
46 
47 // Teuchos includes
48 #include "Teuchos_RCP.hpp"
49 
50 #include "Kokkos_DynRankView.hpp"
51 
56 
57 #include "Teuchos_FancyOStream.hpp"
58 
59 namespace panzer_stk {
60 
61 using Teuchos::RCP;
62 using Teuchos::rcp;
63 
64 // Object describing how to sort a vector of elements using
65 // local ID as the key
66 class LocalIdCompare {
67 public:
68  LocalIdCompare(const RCP<const STK_Interface> & mesh) : mesh_(mesh) {}
69 
70  // Compares two stk mesh entities based on local ID
71  bool operator() (stk::mesh::Entity a,stk::mesh::Entity b)
72  { return mesh_->elementLocalId(a) < mesh_->elementLocalId(b);}
73 
74 private:
75  RCP<const STK_Interface> mesh_;
76 };
77 
79  : stkMeshDB_(stkMeshDB), ownedElementCount_(0)
80 {
81 }
82 
85 {
87 }
88 
90 {
91  elements_ = Teuchos::null;
92 
93  elementBlocks_.clear();
94  elmtLidToConn_.clear();
95  connSize_.clear();
96  elmtToAssociatedElmts_.clear();
97 }
98 
100 {
101  clearLocalElementMapping(); // forget the past
102 
103  // build element block information
105  elements_ = Teuchos::rcp(new std::vector<stk::mesh::Entity>);
106 
107  // defines ordering of blocks
108  std::vector<std::string> blockIds;
109  stkMeshDB_->getElementBlockNames(blockIds);
110 
111  std::size_t blockIndex=0;
112  for(std::vector<std::string>::const_iterator idItr=blockIds.begin();
113  idItr!=blockIds.end();++idItr,++blockIndex) {
114  std::string blockId = *idItr;
115 
116  // grab elements on this block
117  std::vector<stk::mesh::Entity> blockElmts;
118  stkMeshDB_->getMyElements(blockId,blockElmts);
119 
120  // concatenate them into element LID lookup table
121  elements_->insert(elements_->end(),blockElmts.begin(),blockElmts.end());
122 
123  // build block to LID map
124  elementBlocks_[blockId] = Teuchos::rcp(new std::vector<LocalOrdinal>);
125  for(std::size_t i=0;i<blockElmts.size();i++)
126  elementBlocks_[blockId]->push_back(stkMeshDB_->elementLocalId(blockElmts[i]));
127  }
128 
129  ownedElementCount_ = elements_->size();
130 
131  blockIndex=0;
132  for(std::vector<std::string>::const_iterator idItr=blockIds.begin();
133  idItr!=blockIds.end();++idItr,++blockIndex) {
134  std::string blockId = *idItr;
135 
136  // grab elements on this block
137  std::vector<stk::mesh::Entity> blockElmts;
138  stkMeshDB_->getNeighborElements(blockId,blockElmts);
139 
140  // concatenate them into element LID lookup table
141  elements_->insert(elements_->end(),blockElmts.begin(),blockElmts.end());
142 
143  // build block to LID map
144  neighborElementBlocks_[blockId] = Teuchos::rcp(new std::vector<LocalOrdinal>);
145  for(std::size_t i=0;i<blockElmts.size();i++)
146  neighborElementBlocks_[blockId]->push_back(stkMeshDB_->elementLocalId(blockElmts[i]));
147  }
148 
149  // this expensive operation gurantees ordering of local IDs
150  std::sort(elements_->begin(),elements_->end(),LocalIdCompare(stkMeshDB_));
151 
152  // allocate space for element LID to Connectivty map
153  // connectivity size
154  elmtLidToConn_.clear();
155  elmtLidToConn_.resize(elements_->size(),0);
156 
157  connSize_.clear();
158  connSize_.resize(elements_->size(),0);
159 }
160 
161 void
163  LocalOrdinal & nodeIdCnt, LocalOrdinal & edgeIdCnt,
164  LocalOrdinal & faceIdCnt, LocalOrdinal & cellIdCnt,
165  GlobalOrdinal & nodeOffset, GlobalOrdinal & edgeOffset,
166  GlobalOrdinal & faceOffset, GlobalOrdinal & cellOffset) const
167 {
168  // get the global counts for all the nodes, faces, edges and cells
169  GlobalOrdinal maxNodeId = stkMeshDB_->getMaxEntityId(stkMeshDB_->getNodeRank());
170  GlobalOrdinal maxEdgeId = stkMeshDB_->getMaxEntityId(stkMeshDB_->getEdgeRank());
171  GlobalOrdinal maxFaceId = stkMeshDB_->getMaxEntityId(stkMeshDB_->getFaceRank());
172 
173  // compute ID counts for each sub cell type
174  int patternDim = fp.getDimension();
175  switch(patternDim) {
176  case 3:
177  faceIdCnt = fp.getSubcellIndices(2,0).size();
178  // Intentional fall-through.
179  case 2:
180  edgeIdCnt = fp.getSubcellIndices(1,0).size();
181  // Intentional fall-through.
182  case 1:
183  nodeIdCnt = fp.getSubcellIndices(0,0).size();
184  cellIdCnt = fp.getSubcellIndices(patternDim,0).size();
185  break;
186  case 0:
187  default:
188  TEUCHOS_ASSERT(false);
189  };
190 
191  // compute offsets for each sub cell type
192  nodeOffset = 0;
193  edgeOffset = nodeOffset+(maxNodeId+1)*nodeIdCnt;
194  faceOffset = edgeOffset+(maxEdgeId+1)*edgeIdCnt;
195  cellOffset = faceOffset+(maxFaceId+1)*faceIdCnt;
196 
197  // sanity check
198  TEUCHOS_ASSERT(nodeOffset <= edgeOffset
199  && edgeOffset <= faceOffset
200  && faceOffset <= cellOffset);
201 }
202 
205  unsigned subcellRank,
206  LocalOrdinal idCnt,
207  GlobalOrdinal offset)
208 {
209  if(idCnt<=0)
210  return 0 ;
211 
212  // loop over all relations of specified type
213  LocalOrdinal numIds = 0;
214  stk::mesh::BulkData& bulkData = *stkMeshDB_->getBulkData();
215  const stk::mesh::EntityRank rank = static_cast<stk::mesh::EntityRank>(subcellRank);
216  const size_t num_rels = bulkData.num_connectivity(element, rank);
217  stk::mesh::Entity const* relations = bulkData.begin(element, rank);
218  for(std::size_t sc=0; sc<num_rels; ++sc) {
219  stk::mesh::Entity subcell = relations[sc];
220 
221  // add connectivities: adjust for STK indexing craziness
222  for(LocalOrdinal i=0;i<idCnt;i++)
223  connectivity_.push_back(offset+idCnt*(bulkData.identifier(subcell)-1)+i);
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 #ifdef HAVE_EXTRA_TIMERS
248  using Teuchos::TimeMonitor;
249  RCP<Teuchos::TimeMonitor> tM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(std::string("panzer_stk::STKConnManager::buildConnectivity"))));
250 #endif
251 
252  stk::mesh::BulkData& bulkData = *stkMeshDB_->getBulkData();
253 
254  // get element info from STK_Interface
255  // object and build a local element mapping.
257 
258  // Build sub cell ID counts and offsets
259  // ID counts = How many IDs belong on each subcell (number of mesh DOF used)
260  // Offset = What is starting index for subcell ID type?
261  // Global numbering goes like [node ids, edge ids, face ids, cell ids]
262  LocalOrdinal nodeIdCnt=0, edgeIdCnt=0, faceIdCnt=0, cellIdCnt=0;
263  GlobalOrdinal nodeOffset=0, edgeOffset=0, faceOffset=0, cellOffset=0;
264  buildOffsetsAndIdCounts(fp, nodeIdCnt, edgeIdCnt, faceIdCnt, cellIdCnt,
265  nodeOffset, edgeOffset, faceOffset, cellOffset);
266 
267  // std::cout << "node: count = " << nodeIdCnt << ", offset = " << nodeOffset << std::endl;
268  // std::cout << "edge: count = " << edgeIdCnt << ", offset = " << edgeOffset << std::endl;
269  // std::cout << "face: count = " << faceIdCnt << ", offset = " << faceOffset << std::endl;
270  // std::cout << "cell: count = " << cellIdCnt << ", offset = " << cellOffset << std::endl;
271 
272  // loop over elements and build global connectivity
273  for(std::size_t elmtLid=0;elmtLid!=elements_->size();++elmtLid) {
274  GlobalOrdinal numIds = 0;
275  stk::mesh::Entity element = (*elements_)[elmtLid];
276 
277  // get index into connectivity array
278  elmtLidToConn_[elmtLid] = connectivity_.size();
279 
280  // add connecviities for sub cells
281  numIds += addSubcellConnectivities(element,stkMeshDB_->getNodeRank(),nodeIdCnt,nodeOffset);
282  numIds += addSubcellConnectivities(element,stkMeshDB_->getEdgeRank(),edgeIdCnt,edgeOffset);
283  numIds += addSubcellConnectivities(element,stkMeshDB_->getFaceRank(),faceIdCnt,faceOffset);
284 
285  // add connectivity for parent cells
286  if(cellIdCnt>0) {
287  // add connectivities: adjust for STK indexing craziness
288  for(LocalOrdinal i=0;i<cellIdCnt;i++)
289  connectivity_.push_back(cellOffset+cellIdCnt*(bulkData.identifier(element)-1));
290 
291  numIds += cellIdCnt;
292  }
293 
294  connSize_[elmtLid] = numIds;
295  }
296 
297  applyPeriodicBCs( fp, nodeOffset, edgeOffset, faceOffset, cellOffset);
298 
299  // This method does not modify connectivity_. But it should be called here
300  // because the data it initializes should be available at the same time as
301  // connectivity_.
304 }
305 
307 {
308  // walk through the element blocks and figure out which this ID belongs to
309  stk::mesh::Entity element = (*elements_)[localElmtId];
310 
311  return stkMeshDB_->containingBlockId(element);
312 }
313 
315  GlobalOrdinal faceOffset, GlobalOrdinal /* cellOffset */)
316 {
317  using Teuchos::RCP;
318  using Teuchos::rcp;
319 
320 #ifdef HAVE_EXTRA_TIMERS
321  using Teuchos::TimeMonitor;
322  RCP<Teuchos::TimeMonitor> tM = rcp(new TimeMonitor(*TimeMonitor::getNewTimer(std::string("panzer_stk::STKConnManager::applyPeriodicBCs"))));
323 #endif
324 
325  std::pair<Teuchos::RCP<std::vector<std::pair<std::size_t,std::size_t> > >, Teuchos::RCP<std::vector<unsigned int> > > matchedValues
326  = stkMeshDB_->getPeriodicNodePairing();
327 
329  = matchedValues.first;
331  = matchedValues.second;
332 
333  // no matchedNodes means nothing to do!
334  if(matchedNodes==Teuchos::null) return;
335 
336  for(std::size_t m=0;m<matchedNodes->size();m++) {
337  stk::mesh::EntityId oldNodeId = (*matchedNodes)[m].first;
338  std::size_t newNodeId = (*matchedNodes)[m].second;
339 
340  std::vector<stk::mesh::Entity> elements;
341  std::vector<int> localIds;
342 
343  GlobalOrdinal offset0 = 0; // to make numbering consistent with that in PeriodicBC_Matcher
344  GlobalOrdinal offset1 = 0; // offset for dof indexing
345  if((*matchTypes)[m] == 0)
346  offset1 = nodeOffset-offset0;
347  else if((*matchTypes)[m] == 1){
348  offset0 = stkMeshDB_->getMaxEntityId(stkMeshDB_->getNodeRank());
349  offset1 = edgeOffset-offset0;
350  } else if((*matchTypes)[m] == 2){
351  offset0 = stkMeshDB_->getMaxEntityId(stkMeshDB_->getNodeRank())+stkMeshDB_->getMaxEntityId(stkMeshDB_->getEdgeRank());
352  offset1 = faceOffset-offset0;
353  } else
354  TEUCHOS_ASSERT(false);
355 
356  // get relevent elements and node IDs
357  stkMeshDB_->getOwnedElementsSharingNode(oldNodeId-offset0,elements,localIds,(*matchTypes)[m]);
358 
359  // modify global numbering already built for each element
360  for(std::size_t e=0;e<elements.size();e++){
361  modifySubcellConnectivities(fp,elements[e],(*matchTypes)[m],localIds[e],newNodeId,offset1);
362  }
363 
364  }
365 }
366 
369 void STKConnManager::getDofCoords(const std::string & blockId,
370  const panzer::Intrepid2FieldPattern & coordProvider,
371  std::vector<std::size_t> & localCellIds,
372  Kokkos::DynRankView<double,PHX::Device> & points) const
373 {
374  int dim = coordProvider.getDimension();
375  int numIds = coordProvider.numberIds();
376 
377  // grab element vertices
378  Kokkos::DynRankView<double,PHX::Device> vertices;
379  workset_utils::getIdsAndVertices(*stkMeshDB_,blockId,localCellIds,vertices);
380 
381  // setup output array
382  points = Kokkos::DynRankView<double,PHX::Device>("points",localCellIds.size(),numIds,dim);
383  coordProvider.getInterpolatoryCoordinates(vertices,points);
384 }
385 
387 {
388  return ! sidesetsToAssociate_.empty();
389 }
390 
391 void STKConnManager::associateElementsInSideset(const std::string sideset_id)
392 {
393  sidesetsToAssociate_.push_back(sideset_id);
394  sidesetYieldedAssociations_.push_back(false);
395 }
396 
397 inline std::size_t
398 getElementIdx(const std::vector<stk::mesh::Entity>& elements,
399  stk::mesh::Entity const e)
400 {
401  return static_cast<std::size_t>(
402  std::distance(elements.begin(), std::find(elements.begin(), elements.end(), e)));
403 }
404 
406 {
407  stk::mesh::BulkData& bulkData = *stkMeshDB_->getBulkData();
408  elmtToAssociatedElmts_.resize(elements_->size());
409  for (std::size_t i = 0; i < sidesetsToAssociate_.size(); ++i) {
410  std::vector<stk::mesh::Entity> sides;
411  stkMeshDB_->getAllSides(sidesetsToAssociate_[i], sides);
412  sidesetYieldedAssociations_[i] = ! sides.empty();
413  for (std::vector<stk::mesh::Entity>::const_iterator si = sides.begin();
414  si != sides.end(); ++si) {
415  stk::mesh::Entity side = *si;
416  const size_t num_elements = bulkData.num_elements(side);
417  stk::mesh::Entity const* elements = bulkData.begin_elements(side);
418  if (num_elements != 2) {
419  // If relations.size() != 2 for one side in the sideset, then it's true
420  // for all, including the first.
421  TEUCHOS_ASSERT(si == sides.begin());
422  sidesetYieldedAssociations_[i] = false;
423  break;
424  }
425  const std::size_t ea_id = getElementIdx(*elements_, elements[0]),
426  eb_id = getElementIdx(*elements_, elements[1]);
427  elmtToAssociatedElmts_[ea_id].push_back(eb_id);
428  elmtToAssociatedElmts_[eb_id].push_back(ea_id);
429  }
430  }
431 }
432 
433 std::vector<std::string> STKConnManager::
435 {
436  std::vector<std::string> sidesets;
437  for (std::size_t i = 0; i < sidesetYieldedAssociations_.size(); ++i) {
438  int sya, my_sya = sidesetYieldedAssociations_[i] ? 1 : 0;
439  Teuchos::reduceAll(comm, Teuchos::REDUCE_MAX, 1, &my_sya, &sya);
440  if (sya == 0)
441  sidesets.push_back(sidesetsToAssociate_[i]);
442  }
443  return sidesets;
444 }
445 
446 const std::vector<STKConnManager::LocalOrdinal>&
448 {
449  return elmtToAssociatedElmts_[el];
450 }
451 
452 }
virtual int getDimension() const =0
void getIdsAndVertices(const panzer_stk::STK_Interface &mesh, std::string blockId, std::vector< std::size_t > &localIds, ArrayT &vertices)
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
std::vector< bool > sidesetYieldedAssociations_
std::map< std::string, Teuchos::RCP< std::vector< LocalOrdinal > > > neighborElementBlocks_
Teuchos::RCP< const STK_Interface > stkMeshDB_
std::vector< GlobalOrdinal > connectivity_
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
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
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_
void getInterpolatoryCoordinates(Kokkos::DynRankView< double, PHX::Device > &coords) const
LocalOrdinal addSubcellConnectivities(stk::mesh::Entity element, unsigned subcellRank, LocalOrdinal idCnt, GlobalOrdinal offset)
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_