Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_STK_SquareTriMeshFactory.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#include <Teuchos_TimeMonitor.hpp>
45#include <PanzerAdaptersSTK_config.hpp>
46
47using Teuchos::RCP;
48using Teuchos::rcp;
49
50namespace panzer_stk {
51
53{
55}
56
59{
60}
61
63Teuchos::RCP<STK_Interface> SquareTriMeshFactory::buildMesh(stk::ParallelMachine parallelMach) const
64{
65 PANZER_FUNC_TIME_MONITOR("panzer::SquareTriMeshFactory::buildMesh()");
66
67 // build all meta data
68 RCP<STK_Interface> mesh = buildUncommitedMesh(parallelMach);
69
70 // commit meta data
71 mesh->initialize(parallelMach);
72
73 // build bulk data
74 completeMeshConstruction(*mesh,parallelMach);
75
76 return mesh;
77}
78
79Teuchos::RCP<STK_Interface> SquareTriMeshFactory::buildUncommitedMesh(stk::ParallelMachine parallelMach) const
80{
81 PANZER_FUNC_TIME_MONITOR("panzer::SquareTriMeshFactory::buildUncomittedMesh()");
82
83 RCP<STK_Interface> mesh = rcp(new STK_Interface(2));
84
85 machRank_ = stk::parallel_machine_rank(parallelMach);
86 machSize_ = stk::parallel_machine_size(parallelMach);
87
88 if (xProcs_ == -1 && yProcs_ == -1) {
89 // copied from galeri
90 xProcs_ = yProcs_ = Teuchos::as<int>(pow(Teuchos::as<double>(machSize_), 0.5));
91
92 if (xProcs_ * yProcs_ != Teuchos::as<int>(machSize_)) {
93 // Simple method to find a set of processor assignments
94 xProcs_ = yProcs_ = 1;
95
96 // This means that this works correctly up to about maxFactor^2
97 // processors.
98 const int maxFactor = 100;
99
100 int ProcTemp = machSize_;
101 int factors[maxFactor];
102 for (int jj = 0; jj < maxFactor; jj++) factors[jj] = 0;
103 for (int jj = 2; jj < maxFactor; jj++) {
104 bool flag = true;
105 while (flag) {
106 int temp = ProcTemp/jj;
107 if (temp*jj == ProcTemp) {
108 factors[jj]++;
109 ProcTemp = temp;
110
111 } else {
112 flag = false;
113 }
114 }
115 }
116 xProcs_ = ProcTemp;
117 for (int jj = maxFactor-1; jj > 0; jj--) {
118 while (factors[jj] != 0) {
119 if (xProcs_ <= yProcs_) xProcs_ = xProcs_*jj;
120 else yProcs_ = yProcs_*jj;
121 factors[jj]--;
122 }
123 }
124 }
125
126 } else if(xProcs_==-1) {
127 // default x only decomposition
129 yProcs_ = 1;
130 }
131 TEUCHOS_TEST_FOR_EXCEPTION(int(machSize_)!=xProcs_*yProcs_,std::logic_error,
132 "Cannot build SquareTriMeshFactory, the product of \"X Procs\" and \"Y Procs\""
133 " must equal the number of processors.");
135
136 // build meta information: blocks and side set setups
137 buildMetaData(parallelMach,*mesh);
138
139 mesh->addPeriodicBCs(periodicBCVec_);
140 mesh->setBoundingBoxSearchFlag(useBBoxSearch_);
141
142 return mesh;
143}
144
145void SquareTriMeshFactory::completeMeshConstruction(STK_Interface & mesh,stk::ParallelMachine parallelMach) const
146{
147 PANZER_FUNC_TIME_MONITOR("panzer::SquareTriMeshFactory::completeMeshConstruction()");
148
149 if(not mesh.isInitialized())
150 mesh.initialize(parallelMach);
151
152 // add node and element information
153 buildElements(parallelMach,mesh);
154
155 // finish up the edges
156 mesh.buildSubcells();
159 mesh.buildLocalEdgeIDs();
160 }
161
162 // now that edges are built, sidets can be added
163 addSideSets(mesh);
164
165 // add nodesets
166 addNodeSets(mesh);
167
169 addEdgeBlocks(mesh);
170 }
171
172 // calls Stk_MeshFactory::rebalance
173 this->rebalance(mesh);
174}
175
177void SquareTriMeshFactory::setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & paramList)
178{
179 paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
180
181 setMyParamList(paramList);
182
183 x0_ = paramList->get<double>("X0");
184 y0_ = paramList->get<double>("Y0");
185
186 xf_ = paramList->get<double>("Xf");
187 yf_ = paramList->get<double>("Yf");
188
189 xBlocks_ = paramList->get<int>("X Blocks");
190 yBlocks_ = paramList->get<int>("Y Blocks");
191
192 nXElems_ = paramList->get<int>("X Elements");
193 nYElems_ = paramList->get<int>("Y Elements");
194
195 xProcs_ = paramList->get<int>("X Procs");
196 yProcs_ = paramList->get<int>("Y Procs");
197
198 createEdgeBlocks_ = paramList->get<bool>("Create Edge Blocks");
199
200 // read in periodic boundary conditions
201 parsePeriodicBCList(Teuchos::rcpFromRef(paramList->sublist("Periodic BCs")),periodicBCVec_,useBBoxSearch_);
202}
203
205Teuchos::RCP<const Teuchos::ParameterList> SquareTriMeshFactory::getValidParameters() const
206{
207 static RCP<Teuchos::ParameterList> defaultParams;
208
209 // fill with default values
210 if(defaultParams == Teuchos::null) {
211 defaultParams = rcp(new Teuchos::ParameterList);
212
213 defaultParams->set<double>("X0",0.0);
214 defaultParams->set<double>("Y0",0.0);
215
216 defaultParams->set<double>("Xf",1.0);
217 defaultParams->set<double>("Yf",1.0);
218
219 defaultParams->set<int>("X Blocks",1);
220 defaultParams->set<int>("Y Blocks",1);
221
222 defaultParams->set<int>("X Procs",-1);
223 defaultParams->set<int>("Y Procs",1);
224
225 defaultParams->set<int>("X Elements",5);
226 defaultParams->set<int>("Y Elements",5);
227
228 // default to false for backward compatibility
229 defaultParams->set<bool>("Create Edge Blocks",false,"Create edge blocks in the mesh");
230
231 Teuchos::ParameterList & bcs = defaultParams->sublist("Periodic BCs");
232 bcs.set<int>("Count",0); // no default periodic boundary conditions
233 }
234
235 return defaultParams;
236}
237
239{
240 // get valid parameters
241 RCP<Teuchos::ParameterList> validParams = rcp(new Teuchos::ParameterList(*getValidParameters()));
242
243 // set that parameter list
244 setParameterList(validParams);
245
246 /* This is a tri mesh factory so all elements in all element blocks
247 * will be tri3. This means that all the edges will be line2.
248 * The edge block name is hard coded to reflect this.
249 */
251}
252
253void SquareTriMeshFactory::buildMetaData(stk::ParallelMachine /* parallelMach */, STK_Interface & mesh) const
254{
255 typedef shards::Triangle<> TriTopo;
256 const CellTopologyData * ctd = shards::getCellTopologyData<TriTopo>();
257 const CellTopologyData * side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(1,0);
258 const CellTopologyData * edge_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(1,0);
259
260 // build meta data
261 //mesh.setDimension(2);
262 for(int bx=0;bx<xBlocks_;bx++) {
263 for(int by=0;by<yBlocks_;by++) {
264
265 // add this element block
266 {
267 std::stringstream ebPostfix;
268 ebPostfix << "-" << bx << "_" << by;
269
270 // add element blocks
271 mesh.addElementBlock("eblock"+ebPostfix.str(),ctd);
273 mesh.addEdgeBlock("eblock"+ebPostfix.str(),
275 edge_ctd);
276 }
277 }
278
279 }
280 }
281
282 // add sidesets
283 mesh.addSideset("left",side_ctd);
284 mesh.addSideset("right",side_ctd);
285 mesh.addSideset("top",side_ctd);
286 mesh.addSideset("bottom",side_ctd);
287
288 // add nodesets
289 mesh.addNodeset("origin");
290}
291
292void SquareTriMeshFactory::buildElements(stk::ParallelMachine parallelMach,STK_Interface & mesh) const
293{
294 mesh.beginModification();
295 // build each block
296 for(int xBlock=0;xBlock<xBlocks_;xBlock++) {
297 for(int yBlock=0;yBlock<yBlocks_;yBlock++) {
298 buildBlock(parallelMach,xBlock,yBlock,mesh);
299 }
300 }
301 mesh.endModification();
302}
303
304void SquareTriMeshFactory::buildBlock(stk::ParallelMachine /* parallelMach */, int xBlock, int yBlock, STK_Interface& mesh) const
305{
306 // grab this processors rank and machine size
307 std::pair<int,int> sizeAndStartX = determineXElemSizeAndStart(xBlock,xProcs_,machRank_);
308 std::pair<int,int> sizeAndStartY = determineYElemSizeAndStart(yBlock,yProcs_,machRank_);
309
310 int myXElems_start = sizeAndStartX.first;
311 int myXElems_end = myXElems_start+sizeAndStartX.second;
312 int myYElems_start = sizeAndStartY.first;
313 int myYElems_end = myYElems_start+sizeAndStartY.second;
314 int totalXElems = nXElems_*xBlocks_;
315 int totalYElems = nYElems_*yBlocks_;
316
317 double deltaX = (xf_-x0_)/double(totalXElems);
318 double deltaY = (yf_-y0_)/double(totalYElems);
319
320 std::vector<double> coord(2,0.0);
321
322 // build the nodes
323 for(int nx=myXElems_start;nx<myXElems_end+1;++nx) {
324 coord[0] = this->getMeshCoord(nx, deltaX, x0_);
325 for(int ny=myYElems_start;ny<myYElems_end+1;++ny) {
326 coord[1] = this->getMeshCoord(ny, deltaY, y0_);
327
328 mesh.addNode(ny*(totalXElems+1)+nx+1,coord);
329 }
330 }
331
332 std::stringstream blockName;
333 blockName << "eblock-" << xBlock << "_" << yBlock;
334 stk::mesh::Part * block = mesh.getElementBlockPart(blockName.str());
335
336 // build the elements
337 for(int nx=myXElems_start;nx<myXElems_end;++nx) {
338 for(int ny=myYElems_start;ny<myYElems_end;++ny) {
339 stk::mesh::EntityId gid_a = 2*(totalXElems*ny+nx+1)-1;
340 stk::mesh::EntityId gid_b = gid_a+1;
341 std::vector<stk::mesh::EntityId> nodes(3);
342 stk::mesh::EntityId sw,se,ne,nw;
343 sw = nx+1+ny*(totalXElems+1);
344 se = sw+1;
345 ne = se+(totalXElems+1);
346 nw = ne-1;
347
348 nodes[0] = sw;
349 nodes[1] = se;
350 nodes[2] = ne;
351 mesh.addElement(rcp(new ElementDescriptor(gid_a,nodes)),block);
352
353 nodes[0] = sw;
354 nodes[1] = ne;
355 nodes[2] = nw;
356 mesh.addElement(rcp(new ElementDescriptor(gid_b,nodes)),block);
357 }
358 }
359}
360
361std::pair<int,int> SquareTriMeshFactory::determineXElemSizeAndStart(int xBlock,unsigned int size,unsigned int /* rank */) const
362{
363 std::size_t xProcLoc = procTuple_[0];
364 unsigned int minElements = nXElems_/size;
365 unsigned int extra = nXElems_ - minElements*size;
366
367 TEUCHOS_ASSERT(minElements>0);
368
369 // first "extra" elements get an extra column of elements
370 // this determines the starting X index and number of elements
371 int nume=0, start=0;
372 if(xProcLoc<extra) {
373 nume = minElements+1;
374 start = xProcLoc*(minElements+1);
375 }
376 else {
377 nume = minElements;
378 start = extra*(minElements+1)+(xProcLoc-extra)*minElements;
379 }
380
381 return std::make_pair(start+nXElems_*xBlock,nume);
382}
383
384std::pair<int,int> SquareTriMeshFactory::determineYElemSizeAndStart(int yBlock,unsigned int size,unsigned int /* rank */) const
385{
386 std::size_t yProcLoc = procTuple_[1];
387 unsigned int minElements = nYElems_/size;
388 unsigned int extra = nYElems_ - minElements*size;
389
390 TEUCHOS_ASSERT(minElements>0);
391
392 // first "extra" elements get an extra column of elements
393 // this determines the starting X index and number of elements
394 int nume=0, start=0;
395 if(yProcLoc<extra) {
396 nume = minElements+1;
397 start = yProcLoc*(minElements+1);
398 }
399 else {
400 nume = minElements;
401 start = extra*(minElements+1)+(yProcLoc-extra)*minElements;
402 }
403
404 return std::make_pair(start+nYElems_*yBlock,nume);
405}
406
408{
409 mesh.beginModification();
410
411 std::size_t totalXElems = nXElems_*xBlocks_;
412 std::size_t totalYElems = nYElems_*yBlocks_;
413
414 // get all part vectors
415 stk::mesh::Part * left = mesh.getSideset("left");
416 stk::mesh::Part * right = mesh.getSideset("right");
417 stk::mesh::Part * top = mesh.getSideset("top");
418 stk::mesh::Part * bottom = mesh.getSideset("bottom");
419
420 std::vector<stk::mesh::Entity> localElmts;
421 mesh.getMyElements(localElmts);
422
423 // loop over elements adding edges to sidesets
424 std::vector<stk::mesh::Entity>::const_iterator itr;
425 for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
426 stk::mesh::Entity element = (*itr);
427 stk::mesh::EntityId gid = mesh.elementGlobalId(element);
428
429 bool lower = (gid%2 != 0);
430 std::size_t block = lower ? (gid+1)/2 : gid/2;
431 std::size_t nx,ny;
432 ny = (block-1) / totalXElems;
433 nx = block-ny*totalXElems-1;
434
435 // vertical boundaries
437
438 if(nx+1==totalXElems && lower) {
439 stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 1);
440
441 // on the right
442 if(mesh.entityOwnerRank(edge)==machRank_)
443 mesh.addEntityToSideset(edge,right);
444 }
445
446 if(nx==0 && !lower) {
447 stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 2);
448
449 // on the left
450 if(mesh.entityOwnerRank(edge)==machRank_)
451 mesh.addEntityToSideset(edge,left);
452 }
453
454 // horizontal boundaries
456
457 if(ny==0 && lower) {
458 stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 0);
459
460 // on the bottom
461 if(mesh.entityOwnerRank(edge)==machRank_)
462 mesh.addEntityToSideset(edge,bottom);
463 }
464
465 if(ny+1==totalYElems && !lower) {
466 stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 1);
467
468 // on the top
469 if(mesh.entityOwnerRank(edge)==machRank_)
470 mesh.addEntityToSideset(edge,top);
471 }
472 }
473
474 mesh.endModification();
475}
476
478{
479 mesh.beginModification();
480
481 // get all part vectors
482 stk::mesh::Part * origin = mesh.getNodeset("origin");
483
484 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
485 if(machRank_==0)
486 {
487 stk::mesh::Entity node = bulkData->get_entity(mesh.getNodeRank(),1);
488
489 // add zero node to origin node set
490 mesh.addEntityToNodeset(node,origin);
491 }
492
493 mesh.endModification();
494}
495
497{
498 mesh.beginModification();
499
500 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
501 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
502
503 stk::mesh::Part * edge_block = mesh.getEdgeBlock(edgeBlockName_);
504
505 stk::mesh::Selector owned_block = metaData->locally_owned_part();
506
507 std::vector<stk::mesh::Entity> edges;
508 bulkData->get_entities(mesh.getEdgeRank(), owned_block, edges);
509 mesh.addEntitiesToEdgeBlock(edges, edge_block);
510
511 mesh.endModification();
512}
513
515Teuchos::Tuple<std::size_t,2> SquareTriMeshFactory::procRankToProcTuple(std::size_t procRank) const
516{
517 std::size_t i=0,j=0;
518
519 j = procRank/xProcs_;
520 procRank = procRank % xProcs_;
521 i = procRank;
522
523 return Teuchos::tuple(i,j);
524}
525
526} // end panzer_stk
stk::mesh::Entity findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true, const bool buildRefinementSupport=false)
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block part
static const std::string edgeBlockString
stk::mesh::EntityId elementGlobalId(std::size_t lid) const
stk::mesh::Part * getEdgeBlock(const std::string &name) const
get the block part
bool isInitialized() const
Has initialize been called on this mesh object?
stk::mesh::EntityRank getNodeRank() const
void addEdgeBlock(const std::string &elemBlockName, const std::string &edgeBlockName, const stk::topology &topology)
void buildSubcells()
force the mesh to build subcells: edges and faces
void addEntityToNodeset(stk::mesh::Entity entity, stk::mesh::Part *nodeset)
void addElement(const Teuchos::RCP< ElementDescriptor > &ed, stk::mesh::Part *block)
void addNode(stk::mesh::EntityId gid, const std::vector< double > &coord)
void addNodeset(const std::string &name)
void addSideset(const std::string &name, const CellTopologyData *ctData)
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
stk::mesh::Part * getNodeset(const std::string &name) const
unsigned entityOwnerRank(stk::mesh::Entity entity) const
void addEntitiesToEdgeBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *edgeblock)
void addEntityToSideset(stk::mesh::Entity entity, stk::mesh::Part *sideset)
void getMyElements(std::vector< stk::mesh::Entity > &elements) const
void addElementBlock(const std::string &name, const CellTopologyData *ctData)
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
stk::mesh::EntityRank getEdgeRank() const
stk::mesh::Part * getSideset(const std::string &name) const
void rebalance(STK_Interface &mesh) const
static void parsePeriodicBCList(const Teuchos::RCP< Teuchos::ParameterList > &pl, std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > &periodicBC, bool &useBBoxSearch)
double getMeshCoord(const int nx, const double deltaX, const double x0) const
std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > periodicBCVec_
void buildElements(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
From ParameterListAcceptor.
Teuchos::RCP< STK_Interface > buildMesh(stk::ParallelMachine parallelMach) const
Build the mesh object.
Teuchos::Tuple< std::size_t, 2 > procRankToProcTuple(std::size_t procRank) const
what is the 2D tuple describe this processor distribution
void buildMetaData(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
std::pair< int, int > determineYElemSizeAndStart(int yBlock, unsigned int size, unsigned int rank) const
void buildBlock(stk::ParallelMachine machRank, int xBlock, int yBlock, STK_Interface &mesh) const
std::pair< int, int > determineXElemSizeAndStart(int xBlock, unsigned int size, unsigned int rank) const
virtual void completeMeshConstruction(STK_Interface &mesh, stk::ParallelMachine parallelMach) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
From ParameterListAcceptor.
virtual Teuchos::RCP< STK_Interface > buildUncommitedMesh(stk::ParallelMachine parallelMach) const