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) {
89  // default x only decomposition
90  xProcs_ = machSize_;
91  yProcs_ = 1;
92  }
93  TEUCHOS_TEST_FOR_EXCEPTION(int(machSize_)!=xProcs_*yProcs_,std::logic_error,
94  "Cannot build SquareTriMeshFactory, the product of \"X Procs\" and \"Y Procs\""
95  " must equal the number of processors.");
97 
98  // build meta information: blocks and side set setups
99  buildMetaData(parallelMach,*mesh);
100 
101  mesh->addPeriodicBCs(periodicBCVec_);
102 
103  return mesh;
104 }
105 
106 void SquareTriMeshFactory::completeMeshConstruction(STK_Interface & mesh,stk::ParallelMachine parallelMach) const
107 {
108  PANZER_FUNC_TIME_MONITOR("panzer::SquareTriMeshFactory::completeMeshConstruction()");
109 
110  if(not mesh.isInitialized())
111  mesh.initialize(parallelMach);
112 
113  // add node and element information
114  buildElements(parallelMach,mesh);
115 
116  // finish up the edges
117  mesh.buildSubcells();
118  mesh.buildLocalElementIDs();
119 
120  // now that edges are built, sidets can be added
121  addSideSets(mesh);
122 
123  // calls Stk_MeshFactory::rebalance
124  this->rebalance(mesh);
125 }
126 
128 void SquareTriMeshFactory::setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & paramList)
129 {
130  paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
131 
132  setMyParamList(paramList);
133 
134  x0_ = paramList->get<double>("X0");
135  y0_ = paramList->get<double>("Y0");
136 
137  xf_ = paramList->get<double>("Xf");
138  yf_ = paramList->get<double>("Yf");
139 
140  xBlocks_ = paramList->get<int>("X Blocks");
141  yBlocks_ = paramList->get<int>("Y Blocks");
142 
143  nXElems_ = paramList->get<int>("X Elements");
144  nYElems_ = paramList->get<int>("Y Elements");
145 
146  xProcs_ = paramList->get<int>("X Procs");
147  yProcs_ = paramList->get<int>("Y Procs");
148 
149  // read in periodic boundary conditions
150  parsePeriodicBCList(Teuchos::rcpFromRef(paramList->sublist("Periodic BCs")),periodicBCVec_);
151 }
152 
154 Teuchos::RCP<const Teuchos::ParameterList> SquareTriMeshFactory::getValidParameters() const
155 {
156  static RCP<Teuchos::ParameterList> defaultParams;
157 
158  // fill with default values
159  if(defaultParams == Teuchos::null) {
160  defaultParams = rcp(new Teuchos::ParameterList);
161 
162  defaultParams->set<double>("X0",0.0);
163  defaultParams->set<double>("Y0",0.0);
164 
165  defaultParams->set<double>("Xf",1.0);
166  defaultParams->set<double>("Yf",1.0);
167 
168  defaultParams->set<int>("X Blocks",1);
169  defaultParams->set<int>("Y Blocks",1);
170 
171  defaultParams->set<int>("X Procs",-1);
172  defaultParams->set<int>("Y Procs",1);
173 
174  defaultParams->set<int>("X Elements",5);
175  defaultParams->set<int>("Y Elements",5);
176 
177  Teuchos::ParameterList & bcs = defaultParams->sublist("Periodic BCs");
178  bcs.set<int>("Count",0); // no default periodic boundary conditions
179  }
180 
181  return defaultParams;
182 }
183 
185 {
186  // get valid parameters
187  RCP<Teuchos::ParameterList> validParams = rcp(new Teuchos::ParameterList(*getValidParameters()));
188 
189  // set that parameter list
190  setParameterList(validParams);
191 }
192 
193 void SquareTriMeshFactory::buildMetaData(stk::ParallelMachine parallelMach, STK_Interface & mesh) const
194 {
195  typedef shards::Triangle<> TriTopo;
196  const CellTopologyData * ctd = shards::getCellTopologyData<TriTopo>();
197  const CellTopologyData * side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(1,0);
198 
199  // build meta data
200  //mesh.setDimension(2);
201  for(int bx=0;bx<xBlocks_;bx++) {
202  for(int by=0;by<yBlocks_;by++) {
203 
204  // add this element block
205  {
206  std::stringstream ebPostfix;
207  ebPostfix << "-" << bx << "_" << by;
208 
209  // add element blocks
210  mesh.addElementBlock("eblock"+ebPostfix.str(),ctd);
211  }
212 
213  }
214  }
215 
216  // add sidesets
217  mesh.addSideset("left",side_ctd);
218  mesh.addSideset("right",side_ctd);
219  mesh.addSideset("top",side_ctd);
220  mesh.addSideset("bottom",side_ctd);
221 }
222 
223 void SquareTriMeshFactory::buildElements(stk::ParallelMachine parallelMach,STK_Interface & mesh) const
224 {
225  mesh.beginModification();
226  // build each block
227  for(int xBlock=0;xBlock<xBlocks_;xBlock++) {
228  for(int yBlock=0;yBlock<yBlocks_;yBlock++) {
229  buildBlock(parallelMach,xBlock,yBlock,mesh);
230  }
231  }
232  mesh.endModification();
233 }
234 
235 void SquareTriMeshFactory::buildBlock(stk::ParallelMachine parallelMach,int xBlock,int yBlock,STK_Interface & mesh) const
236 {
237  // grab this processors rank and machine size
238  std::pair<int,int> sizeAndStartX = determineXElemSizeAndStart(xBlock,xProcs_,machRank_);
239  std::pair<int,int> sizeAndStartY = determineYElemSizeAndStart(yBlock,yProcs_,machRank_);
240 
241  int myXElems_start = sizeAndStartX.first;
242  int myXElems_end = myXElems_start+sizeAndStartX.second;
243  int myYElems_start = sizeAndStartY.first;
244  int myYElems_end = myYElems_start+sizeAndStartY.second;
245  int totalXElems = nXElems_*xBlocks_;
246  int totalYElems = nYElems_*yBlocks_;
247 
248  double deltaX = (xf_-x0_)/double(totalXElems);
249  double deltaY = (yf_-y0_)/double(totalYElems);
250 
251  std::vector<double> coord(2,0.0);
252 
253  // build the nodes
254  for(int nx=myXElems_start;nx<myXElems_end+1;++nx) {
255  coord[0] = double(nx)*deltaX+x0_;
256  for(int ny=myYElems_start;ny<myYElems_end+1;++ny) {
257  coord[1] = double(ny)*deltaY+y0_;
258 
259  mesh.addNode(ny*(totalXElems+1)+nx+1,coord);
260  }
261  }
262 
263  std::stringstream blockName;
264  blockName << "eblock-" << xBlock << "_" << yBlock;
265  stk::mesh::Part * block = mesh.getElementBlockPart(blockName.str());
266 
267  // build the elements
268  for(int nx=myXElems_start;nx<myXElems_end;++nx) {
269  for(int ny=myYElems_start;ny<myYElems_end;++ny) {
270  stk::mesh::EntityId gid_a = 2*(totalXElems*ny+nx+1)-1;
271  stk::mesh::EntityId gid_b = gid_a+1;
272  std::vector<stk::mesh::EntityId> nodes(3);
273  stk::mesh::EntityId sw,se,ne,nw;
274  sw = nx+1+ny*(totalXElems+1);
275  se = sw+1;
276  ne = se+(totalXElems+1);
277  nw = ne-1;
278 
279  nodes[0] = sw;
280  nodes[1] = se;
281  nodes[2] = ne;
282  mesh.addElement(rcp(new ElementDescriptor(gid_a,nodes)),block);
283 
284  nodes[0] = sw;
285  nodes[1] = ne;
286  nodes[2] = nw;
287  mesh.addElement(rcp(new ElementDescriptor(gid_b,nodes)),block);
288  }
289  }
290 }
291 
292 std::pair<int,int> SquareTriMeshFactory::determineXElemSizeAndStart(int xBlock,unsigned int size,unsigned int rank) const
293 {
294  std::size_t xProcLoc = procTuple_[0];
295  unsigned int minElements = nXElems_/size;
296  unsigned int extra = nXElems_ - minElements*size;
297 
298  TEUCHOS_ASSERT(minElements>0);
299 
300  // first "extra" elements get an extra column of elements
301  // this determines the starting X index and number of elements
302  int nume=0, start=0;
303  if(xProcLoc<extra) {
304  nume = minElements+1;
305  start = xProcLoc*(minElements+1);
306  }
307  else {
308  nume = minElements;
309  start = extra*(minElements+1)+(xProcLoc-extra)*minElements;
310  }
311 
312  return std::make_pair(start+nXElems_*xBlock,nume);
313 }
314 
315 std::pair<int,int> SquareTriMeshFactory::determineYElemSizeAndStart(int yBlock,unsigned int size,unsigned int rank) const
316 {
317  std::size_t yProcLoc = procTuple_[1];
318  unsigned int minElements = nYElems_/size;
319  unsigned int extra = nYElems_ - minElements*size;
320 
321  TEUCHOS_ASSERT(minElements>0);
322 
323  // first "extra" elements get an extra column of elements
324  // this determines the starting X index and number of elements
325  int nume=0, start=0;
326  if(yProcLoc<extra) {
327  nume = minElements+1;
328  start = yProcLoc*(minElements+1);
329  }
330  else {
331  nume = minElements;
332  start = extra*(minElements+1)+(yProcLoc-extra)*minElements;
333  }
334 
335  return std::make_pair(start+nYElems_*yBlock,nume);
336 }
337 
339 {
340  mesh.beginModification();
341 
342  std::size_t totalXElems = nXElems_*xBlocks_;
343  std::size_t totalYElems = nYElems_*yBlocks_;
344 
345  // get all part vectors
346  stk::mesh::Part * left = mesh.getSideset("left");
347  stk::mesh::Part * right = mesh.getSideset("right");
348  stk::mesh::Part * top = mesh.getSideset("top");
349  stk::mesh::Part * bottom = mesh.getSideset("bottom");
350 
351  std::vector<stk::mesh::Entity> localElmts;
352  mesh.getMyElements(localElmts);
353 
354  // loop over elements adding edges to sidesets
355  std::vector<stk::mesh::Entity>::const_iterator itr;
356  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
357  stk::mesh::Entity element = (*itr);
358  stk::mesh::EntityId gid = mesh.elementGlobalId(element);
359 
360  bool lower = (gid%2 != 0);
361  std::size_t block = lower ? (gid+1)/2 : gid/2;
362  std::size_t nx,ny;
363  ny = (block-1) / totalXElems;
364  nx = block-ny*totalXElems-1;
365 
366  // vertical boundaries
368 
369  if(nx+1==totalXElems && lower) {
370  stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 1);
371 
372  // on the right
373  if(mesh.entityOwnerRank(edge)==machRank_)
374  mesh.addEntityToSideset(edge,right);
375  }
376 
377  if(nx==0 && !lower) {
378  stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 2);
379 
380  // on the left
381  if(mesh.entityOwnerRank(edge)==machRank_)
382  mesh.addEntityToSideset(edge,left);
383  }
384 
385  // horizontal boundaries
387 
388  if(ny==0 && lower) {
389  stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 0);
390 
391  // on the bottom
392  if(mesh.entityOwnerRank(edge)==machRank_)
393  mesh.addEntityToSideset(edge,bottom);
394  }
395 
396  if(ny+1==totalYElems && !lower) {
397  stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 1);
398 
399  // on the top
400  if(mesh.entityOwnerRank(edge)==machRank_)
401  mesh.addEntityToSideset(edge,top);
402  }
403  }
404 
405  mesh.endModification();
406 }
407 
409 Teuchos::Tuple<std::size_t,2> SquareTriMeshFactory::procRankToProcTuple(std::size_t procRank) const
410 {
411  std::size_t i=0,j=0;
412 
413  j = procRank/xProcs_;
414  procRank = procRank % xProcs_;
415  i = procRank;
416 
417  return Teuchos::tuple(i,j);
418 }
419 
420 } // end panzer_stk
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block count
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)
void buildMetaData(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
Teuchos::Tuple< std::size_t, 2 > procRankToProcTuple(std::size_t procRank) const
what is the 2D tuple describe this processor distribution
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
void buildBlock(stk::ParallelMachine machRank, int xBlock, int yBlock, STK_Interface &mesh) 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
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
virtual void completeMeshConstruction(STK_Interface &mesh, stk::ParallelMachine parallelMach) const
void addNode(stk::mesh::EntityId gid, const std::vector< double > &coord)
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)
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
static void parsePeriodicBCList(const Teuchos::RCP< Teuchos::ParameterList > &pl, std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > &periodicBC)