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