Panzer Version of the Day
Loading...
Searching...
No Matches
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
59namespace panzer_stk {
60
61using Teuchos::RCP;
62using Teuchos::rcp;
63
64// Object describing how to sort a vector of elements using
65// local ID as the key
66class LocalIdCompare {
67public:
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
74private:
75 RCP<const STK_Interface> mesh_;
76};
77
78STKConnManager::STKConnManager(const Teuchos::RCP<const STK_Interface> & stkMeshDB)
79 : stkMeshDB_(stkMeshDB), ownedElementCount_(0)
80{
81}
82
83Teuchos::RCP<panzer::ConnManager>
85{
86 return Teuchos::rcp(new STKConnManager(stkMeshDB_));
87}
88
90{
91 elements_ = Teuchos::null;
92
93 elementBlocks_.clear();
94 elmtLidToConn_.clear();
95 connSize_.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
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
161void
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
230void
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
328 Teuchos::RCP<std::vector<std::pair<std::size_t,std::size_t> > > matchedNodes
329 = matchedValues.first;
330 Teuchos::RCP<std::vector<unsigned int> > matchTypes
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(stk::mesh::EntityId(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
369void 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
391void STKConnManager::associateElementsInSideset(const std::string sideset_id)
392{
393 sidesetsToAssociate_.push_back(sideset_id);
394 sidesetYieldedAssociations_.push_back(false);
395}
396
397inline std::size_t
398getElementIdx(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());
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
433std::vector<std::string> STKConnManager::
434checkAssociateElementsInSidesets(const Teuchos::Comm<int>& comm) const
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
446const std::vector<STKConnManager::LocalOrdinal>&
448{
449 return elmtToAssociatedElmts_[el];
450}
451
452}
virtual int getDimension() const =0
virtual int numberIds() const
virtual const std::vector< int > & getSubcellIndices(int dim, int cellIndex) const =0
void getInterpolatoryCoordinates(Kokkos::DynRankView< double, PHX::Device > &coords) const
std::vector< std::string > checkAssociateElementsInSidesets(const Teuchos::Comm< int > &comm) const
void applyPeriodicBCs(const panzer::FieldPattern &fp, GlobalOrdinal nodeOffset, GlobalOrdinal edgeOffset, GlobalOrdinal faceOffset, GlobalOrdinal cellOffset)
Teuchos::RCP< std::vector< stk::mesh::Entity > > elements_
virtual Teuchos::RCP< panzer::ConnManager > noConnectivityClone() const
virtual const panzer::GlobalOrdinal * getConnectivity(LocalOrdinal localElmtId) const
std::map< std::string, Teuchos::RCP< std::vector< LocalOrdinal > > > neighborElementBlocks_
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< std::vector< LocalOrdinal > > elmtToAssociatedElmts_
std::vector< LocalOrdinal > connSize_
void buildOffsetsAndIdCounts(const panzer::FieldPattern &fp, LocalOrdinal &nodeIdCnt, LocalOrdinal &edgeIdCnt, LocalOrdinal &faceIdCnt, LocalOrdinal &cellIdCnt, GlobalOrdinal &nodeOffset, GlobalOrdinal &edgeOffset, GlobalOrdinal &faceOffset, GlobalOrdinal &cellOffset) const
STKConnManager(const Teuchos::RCP< const STK_Interface > &stkMeshDB)
panzer::ConnManager::GlobalOrdinal GlobalOrdinal
virtual void buildConnectivity(const panzer::FieldPattern &fp)
std::vector< bool > sidesetYieldedAssociations_
virtual bool hasAssociatedNeighbors() const
std::vector< std::string > sidesetsToAssociate_
virtual const std::vector< LocalOrdinal > & getAssociatedNeighbors(const LocalOrdinal &el) const
std::map< std::string, Teuchos::RCP< std::vector< LocalOrdinal > > > elementBlocks_
Teuchos::RCP< const STK_Interface > stkMeshDB_
virtual std::string getBlockId(LocalOrdinal localElmtId) const
std::vector< LocalOrdinal > elmtLidToConn_
void associateElementsInSideset(const std::string sideset_id)
void modifySubcellConnectivities(const panzer::FieldPattern &fp, stk::mesh::Entity element, unsigned subcellRank, unsigned subcellId, GlobalOrdinal newId, GlobalOrdinal offset)
std::vector< GlobalOrdinal > connectivity_
LocalOrdinal addSubcellConnectivities(stk::mesh::Entity element, unsigned subcellRank, LocalOrdinal idCnt, GlobalOrdinal offset)
void getIdsAndVertices(const panzer_stk::STK_Interface &mesh, std::string blockId, std::vector< std::size_t > &localIds, ArrayT &vertices)
std::size_t getElementIdx(const std::vector< stk::mesh::Entity > &elements, stk::mesh::Entity const e)