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