Panzer  Version of the Day
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 
47 using Teuchos::RCP;
48 using Teuchos::rcp;
49 
50 namespace panzer_stk {
51 
53 {
55 }
56 
59 {
60 }
61 
63 Teuchos::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 
79 Teuchos::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
128  xProcs_ = machSize_;
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 
141  return mesh;
142 }
143 
144 void SquareTriMeshFactory::completeMeshConstruction(STK_Interface & mesh,stk::ParallelMachine parallelMach) const
145 {
146  PANZER_FUNC_TIME_MONITOR("panzer::SquareTriMeshFactory::completeMeshConstruction()");
147 
148  if(not mesh.isInitialized())
149  mesh.initialize(parallelMach);
150 
151  // add node and element information
152  buildElements(parallelMach,mesh);
153 
154  // finish up the edges
155  mesh.buildSubcells();
156  mesh.buildLocalElementIDs();
157  if(createEdgeBlocks_) {
158  mesh.buildLocalEdgeIDs();
159  }
160 
161  // now that edges are built, sidets can be added
162  addSideSets(mesh);
163 
164  // add nodesets
165  addNodeSets(mesh);
166 
167  if(createEdgeBlocks_) {
168  addEdgeBlocks(mesh);
169  }
170 
171  // calls Stk_MeshFactory::rebalance
172  this->rebalance(mesh);
173 }
174 
176 void SquareTriMeshFactory::setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & paramList)
177 {
178  paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
179 
180  setMyParamList(paramList);
181 
182  x0_ = paramList->get<double>("X0");
183  y0_ = paramList->get<double>("Y0");
184 
185  xf_ = paramList->get<double>("Xf");
186  yf_ = paramList->get<double>("Yf");
187 
188  xBlocks_ = paramList->get<int>("X Blocks");
189  yBlocks_ = paramList->get<int>("Y Blocks");
190 
191  nXElems_ = paramList->get<int>("X Elements");
192  nYElems_ = paramList->get<int>("Y Elements");
193 
194  xProcs_ = paramList->get<int>("X Procs");
195  yProcs_ = paramList->get<int>("Y Procs");
196 
197  createEdgeBlocks_ = paramList->get<bool>("Create Edge Blocks");
198 
199  // read in periodic boundary conditions
200  parsePeriodicBCList(Teuchos::rcpFromRef(paramList->sublist("Periodic BCs")),periodicBCVec_);
201 }
202 
204 Teuchos::RCP<const Teuchos::ParameterList> SquareTriMeshFactory::getValidParameters() const
205 {
206  static RCP<Teuchos::ParameterList> defaultParams;
207 
208  // fill with default values
209  if(defaultParams == Teuchos::null) {
210  defaultParams = rcp(new Teuchos::ParameterList);
211 
212  defaultParams->set<double>("X0",0.0);
213  defaultParams->set<double>("Y0",0.0);
214 
215  defaultParams->set<double>("Xf",1.0);
216  defaultParams->set<double>("Yf",1.0);
217 
218  defaultParams->set<int>("X Blocks",1);
219  defaultParams->set<int>("Y Blocks",1);
220 
221  defaultParams->set<int>("X Procs",-1);
222  defaultParams->set<int>("Y Procs",1);
223 
224  defaultParams->set<int>("X Elements",5);
225  defaultParams->set<int>("Y Elements",5);
226 
227  // default to false for backward compatibility
228  defaultParams->set<bool>("Create Edge Blocks",false,"Create edge blocks in the mesh");
229 
230  Teuchos::ParameterList & bcs = defaultParams->sublist("Periodic BCs");
231  bcs.set<int>("Count",0); // no default periodic boundary conditions
232  }
233 
234  return defaultParams;
235 }
236 
238 {
239  // get valid parameters
240  RCP<Teuchos::ParameterList> validParams = rcp(new Teuchos::ParameterList(*getValidParameters()));
241 
242  // set that parameter list
243  setParameterList(validParams);
244 }
245 
246 void SquareTriMeshFactory::buildMetaData(stk::ParallelMachine /* parallelMach */, STK_Interface & mesh) const
247 {
248  typedef shards::Triangle<> TriTopo;
249  const CellTopologyData * ctd = shards::getCellTopologyData<TriTopo>();
250  const CellTopologyData * side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(1,0);
251 
252  // build meta data
253  //mesh.setDimension(2);
254  for(int bx=0;bx<xBlocks_;bx++) {
255  for(int by=0;by<yBlocks_;by++) {
256 
257  // add this element block
258  {
259  std::stringstream ebPostfix;
260  ebPostfix << "-" << bx << "_" << by;
261 
262  // add element blocks
263  mesh.addElementBlock("eblock"+ebPostfix.str(),ctd);
264  }
265 
266  }
267  }
268 
269  // add sidesets
270  mesh.addSideset("left",side_ctd);
271  mesh.addSideset("right",side_ctd);
272  mesh.addSideset("top",side_ctd);
273  mesh.addSideset("bottom",side_ctd);
274 
275  // add nodesets
276  mesh.addNodeset("origin");
277 
278  if(createEdgeBlocks_) {
279  const CellTopologyData * edge_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(1,0);
281  }
282 }
283 
284 void SquareTriMeshFactory::buildElements(stk::ParallelMachine parallelMach,STK_Interface & mesh) const
285 {
286  mesh.beginModification();
287  // build each block
288  for(int xBlock=0;xBlock<xBlocks_;xBlock++) {
289  for(int yBlock=0;yBlock<yBlocks_;yBlock++) {
290  buildBlock(parallelMach,xBlock,yBlock,mesh);
291  }
292  }
293  mesh.endModification();
294 }
295 
296 void SquareTriMeshFactory::buildBlock(stk::ParallelMachine /* parallelMach */, int xBlock, int yBlock, STK_Interface& mesh) const
297 {
298  // grab this processors rank and machine size
299  std::pair<int,int> sizeAndStartX = determineXElemSizeAndStart(xBlock,xProcs_,machRank_);
300  std::pair<int,int> sizeAndStartY = determineYElemSizeAndStart(yBlock,yProcs_,machRank_);
301 
302  int myXElems_start = sizeAndStartX.first;
303  int myXElems_end = myXElems_start+sizeAndStartX.second;
304  int myYElems_start = sizeAndStartY.first;
305  int myYElems_end = myYElems_start+sizeAndStartY.second;
306  int totalXElems = nXElems_*xBlocks_;
307  int totalYElems = nYElems_*yBlocks_;
308 
309  double deltaX = (xf_-x0_)/double(totalXElems);
310  double deltaY = (yf_-y0_)/double(totalYElems);
311 
312  std::vector<double> coord(2,0.0);
313 
314  // build the nodes
315  for(int nx=myXElems_start;nx<myXElems_end+1;++nx) {
316  coord[0] = this->getMeshCoord(nx, deltaX, x0_);
317  for(int ny=myYElems_start;ny<myYElems_end+1;++ny) {
318  coord[1] = this->getMeshCoord(ny, deltaY, y0_);
319 
320  mesh.addNode(ny*(totalXElems+1)+nx+1,coord);
321  }
322  }
323 
324  std::stringstream blockName;
325  blockName << "eblock-" << xBlock << "_" << yBlock;
326  stk::mesh::Part * block = mesh.getElementBlockPart(blockName.str());
327 
328  // build the elements
329  for(int nx=myXElems_start;nx<myXElems_end;++nx) {
330  for(int ny=myYElems_start;ny<myYElems_end;++ny) {
331  stk::mesh::EntityId gid_a = 2*(totalXElems*ny+nx+1)-1;
332  stk::mesh::EntityId gid_b = gid_a+1;
333  std::vector<stk::mesh::EntityId> nodes(3);
334  stk::mesh::EntityId sw,se,ne,nw;
335  sw = nx+1+ny*(totalXElems+1);
336  se = sw+1;
337  ne = se+(totalXElems+1);
338  nw = ne-1;
339 
340  nodes[0] = sw;
341  nodes[1] = se;
342  nodes[2] = ne;
343  mesh.addElement(rcp(new ElementDescriptor(gid_a,nodes)),block);
344 
345  nodes[0] = sw;
346  nodes[1] = ne;
347  nodes[2] = nw;
348  mesh.addElement(rcp(new ElementDescriptor(gid_b,nodes)),block);
349  }
350  }
351 }
352 
353 std::pair<int,int> SquareTriMeshFactory::determineXElemSizeAndStart(int xBlock,unsigned int size,unsigned int /* rank */) const
354 {
355  std::size_t xProcLoc = procTuple_[0];
356  unsigned int minElements = nXElems_/size;
357  unsigned int extra = nXElems_ - minElements*size;
358 
359  TEUCHOS_ASSERT(minElements>0);
360 
361  // first "extra" elements get an extra column of elements
362  // this determines the starting X index and number of elements
363  int nume=0, start=0;
364  if(xProcLoc<extra) {
365  nume = minElements+1;
366  start = xProcLoc*(minElements+1);
367  }
368  else {
369  nume = minElements;
370  start = extra*(minElements+1)+(xProcLoc-extra)*minElements;
371  }
372 
373  return std::make_pair(start+nXElems_*xBlock,nume);
374 }
375 
376 std::pair<int,int> SquareTriMeshFactory::determineYElemSizeAndStart(int yBlock,unsigned int size,unsigned int /* rank */) const
377 {
378  std::size_t yProcLoc = procTuple_[1];
379  unsigned int minElements = nYElems_/size;
380  unsigned int extra = nYElems_ - minElements*size;
381 
382  TEUCHOS_ASSERT(minElements>0);
383 
384  // first "extra" elements get an extra column of elements
385  // this determines the starting X index and number of elements
386  int nume=0, start=0;
387  if(yProcLoc<extra) {
388  nume = minElements+1;
389  start = yProcLoc*(minElements+1);
390  }
391  else {
392  nume = minElements;
393  start = extra*(minElements+1)+(yProcLoc-extra)*minElements;
394  }
395 
396  return std::make_pair(start+nYElems_*yBlock,nume);
397 }
398 
400 {
401  mesh.beginModification();
402 
403  std::size_t totalXElems = nXElems_*xBlocks_;
404  std::size_t totalYElems = nYElems_*yBlocks_;
405 
406  // get all part vectors
407  stk::mesh::Part * left = mesh.getSideset("left");
408  stk::mesh::Part * right = mesh.getSideset("right");
409  stk::mesh::Part * top = mesh.getSideset("top");
410  stk::mesh::Part * bottom = mesh.getSideset("bottom");
411 
412  std::vector<stk::mesh::Entity> localElmts;
413  mesh.getMyElements(localElmts);
414 
415  // loop over elements adding edges to sidesets
416  std::vector<stk::mesh::Entity>::const_iterator itr;
417  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
418  stk::mesh::Entity element = (*itr);
419  stk::mesh::EntityId gid = mesh.elementGlobalId(element);
420 
421  bool lower = (gid%2 != 0);
422  std::size_t block = lower ? (gid+1)/2 : gid/2;
423  std::size_t nx,ny;
424  ny = (block-1) / totalXElems;
425  nx = block-ny*totalXElems-1;
426 
427  // vertical boundaries
429 
430  if(nx+1==totalXElems && lower) {
431  stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 1);
432 
433  // on the right
434  if(mesh.entityOwnerRank(edge)==machRank_)
435  mesh.addEntityToSideset(edge,right);
436  }
437 
438  if(nx==0 && !lower) {
439  stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 2);
440 
441  // on the left
442  if(mesh.entityOwnerRank(edge)==machRank_)
443  mesh.addEntityToSideset(edge,left);
444  }
445 
446  // horizontal boundaries
448 
449  if(ny==0 && lower) {
450  stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 0);
451 
452  // on the bottom
453  if(mesh.entityOwnerRank(edge)==machRank_)
454  mesh.addEntityToSideset(edge,bottom);
455  }
456 
457  if(ny+1==totalYElems && !lower) {
458  stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 1);
459 
460  // on the top
461  if(mesh.entityOwnerRank(edge)==machRank_)
462  mesh.addEntityToSideset(edge,top);
463  }
464  }
465 
466  mesh.endModification();
467 }
468 
470 {
471  mesh.beginModification();
472 
473  // get all part vectors
474  stk::mesh::Part * origin = mesh.getNodeset("origin");
475 
476  Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
477  if(machRank_==0)
478  {
479  stk::mesh::Entity node = bulkData->get_entity(mesh.getNodeRank(),1);
480 
481  // add zero node to origin node set
482  mesh.addEntityToNodeset(node,origin);
483  }
484 
485  mesh.endModification();
486 }
487 
489 {
490  mesh.beginModification();
491 
492  stk::mesh::Part * edge_block = mesh.getEdgeBlock(panzer_stk::STK_Interface::edgeBlockString);
493 
494  Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
495  Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
496 
497  std::vector<stk::mesh::Entity> edges;
498  bulkData->get_entities(mesh.getEdgeRank(),metaData->locally_owned_part(),edges);
499  mesh.addEntitiesToEdgeBlock(edges, edge_block);
500 
501  mesh.endModification();
502 }
503 
505 Teuchos::Tuple<std::size_t,2> SquareTriMeshFactory::procRankToProcTuple(std::size_t procRank) const
506 {
507  std::size_t i=0,j=0;
508 
509  j = procRank/xProcs_;
510  procRank = procRank % xProcs_;
511  i = procRank;
512 
513  return Teuchos::tuple(i,j);
514 }
515 
516 } // end panzer_stk
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block part
void addNodeset(const std::string &name)
void addEntityToNodeset(stk::mesh::Entity entity, stk::mesh::Part *nodeset)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
From ParameterListAcceptor.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
From ParameterListAcceptor.
void getMyElements(std::vector< stk::mesh::Entity > &elements) const
void addEntityToSideset(stk::mesh::Entity entity, stk::mesh::Part *sideset)
stk::mesh::EntityRank getNodeRank() const
void buildMetaData(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
Teuchos::Tuple< std::size_t, 2 > procRankToProcTuple(std::size_t procRank) const
what is the 2D tuple describe this processor distribution
stk::mesh::Entity findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const
void buildBlock(stk::ParallelMachine machRank, int xBlock, int yBlock, STK_Interface &mesh) const
stk::mesh::Part * getEdgeBlock(const std::string &name) const
get the block part
stk::mesh::EntityRank getEdgeRank() const
stk::mesh::Part * getSideset(const std::string &name) const
void buildElements(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
void addElement(const Teuchos::RCP< ElementDescriptor > &ed, stk::mesh::Part *block)
stk::mesh::EntityId elementGlobalId(std::size_t lid) const
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true, const bool buildRefinementSupport=false)
double getMeshCoord(const int nx, const double deltaX, const double x0) const
bool isInitialized() const
Has initialize been called on this mesh object?
static const std::string edgeBlockString
std::pair< int, int > determineYElemSizeAndStart(int yBlock, unsigned int size, unsigned int rank) const
virtual void completeMeshConstruction(STK_Interface &mesh, stk::ParallelMachine parallelMach) const
void addNode(stk::mesh::EntityId gid, const std::vector< double > &coord)
void addEntitiesToEdgeBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *edgeblock)
void buildSubcells()
force the mesh to build subcells: edges and faces
std::pair< int, int > determineXElemSizeAndStart(int xBlock, unsigned int size, unsigned int rank) const
void addSideset(const std::string &name, const CellTopologyData *ctData)
void addEdgeBlock(const std::string &name, const CellTopologyData *ctData)
Teuchos::RCP< STK_Interface > buildMesh(stk::ParallelMachine parallelMach) const
Build the mesh object.
std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > periodicBCVec_
unsigned entityOwnerRank(stk::mesh::Entity entity) const
void rebalance(STK_Interface &mesh) const
void addElementBlock(const std::string &name, const CellTopologyData *ctData)
virtual Teuchos::RCP< STK_Interface > buildUncommitedMesh(stk::ParallelMachine parallelMach) const
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)