Panzer  Version of the Day
Panzer_STK_CubeHexMeshFactory.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 #include <stk_mesh/base/FEMHelpers.hpp>
47 
48 using Teuchos::RCP;
49 using Teuchos::rcp;
50 
51 namespace panzer_stk {
52 
54 {
56 }
57 
60 {
61 }
62 
64 Teuchos::RCP<STK_Interface> CubeHexMeshFactory::buildMesh(stk::ParallelMachine parallelMach) const
65 {
66  PANZER_FUNC_TIME_MONITOR("panzer::CubeHexMeshFactory::buildMesh()");
67 
68  // build all meta data
69  RCP<STK_Interface> mesh = buildUncommitedMesh(parallelMach);
70 
71  // commit meta data
72  mesh->initialize(parallelMach);
73 
74  // build bulk data
75  completeMeshConstruction(*mesh,parallelMach);
76 
77  return mesh;
78 }
79 
80 Teuchos::RCP<STK_Interface> CubeHexMeshFactory::buildUncommitedMesh(stk::ParallelMachine parallelMach) const
81 {
82  PANZER_FUNC_TIME_MONITOR("panzer::CubeHexMeshFactory::buildUncomittedMesh()");
83 
84  RCP<STK_Interface> mesh = rcp(new STK_Interface(3));
85 
86  machRank_ = stk::parallel_machine_rank(parallelMach);
87  machSize_ = stk::parallel_machine_size(parallelMach);
88 
89  if (xProcs_ == -1 && yProcs_ == -1 && zProcs_ == -1) {
90  // copied from galeri
91  xProcs_ = yProcs_ = zProcs_ = Teuchos::as<int>(pow(Teuchos::as<double>(machSize_), 0.333334));
92 
93  if (xProcs_ * yProcs_ * zProcs_ != Teuchos::as<int>(machSize_)) {
94  // Simple method to find a set of processor assignments
95  xProcs_ = yProcs_ = zProcs_ = 1;
96 
97  // This means that this works correctly up to about maxFactor^3
98  // processors.
99  const int maxFactor = 50;
100 
101  int ProcTemp = machSize_;
102  int factors[maxFactor];
103  for (int jj = 0; jj < maxFactor; jj++) factors[jj] = 0;
104  for (int jj = 2; jj < maxFactor; jj++) {
105  bool flag = true;
106  while (flag) {
107  int temp = ProcTemp/jj;
108  if (temp*jj == ProcTemp) {
109  factors[jj]++;
110  ProcTemp = temp;
111 
112  } else {
113  flag = false;
114  }
115  }
116  }
117  xProcs_ = ProcTemp;
118  for (int jj = maxFactor-1; jj > 0; jj--) {
119  while (factors[jj] != 0) {
120  if ((xProcs_ <= yProcs_) && (xProcs_ <= zProcs_)) xProcs_ = xProcs_*jj;
121  else if ((yProcs_ <= xProcs_) && (yProcs_ <= zProcs_)) yProcs_ = yProcs_*jj;
122  else zProcs_ = zProcs_*jj;
123  factors[jj]--;
124  }
125  }
126  }
127 
128  } else if(xProcs_==-1) {
129  // default x only decomposition
130  xProcs_ = machSize_;
131  yProcs_ = 1;
132  zProcs_ = 1;
133  }
134  TEUCHOS_TEST_FOR_EXCEPTION(int(machSize_)!=xProcs_*yProcs_*zProcs_,std::logic_error,
135  "Cannot build CubeHexMeshFactory, the product of \"X Procs\", \"Y Procs\", and \"Z Procs\""
136  " must equal the number of processors.");
138 
139  // build meta information: blocks and side set setups
140  buildMetaData(parallelMach,*mesh);
141 
142  mesh->addPeriodicBCs(periodicBCVec_);
143 
144  return mesh;
145 }
146 
147 void CubeHexMeshFactory::completeMeshConstruction(STK_Interface & mesh,stk::ParallelMachine parallelMach) const
148 {
149  PANZER_FUNC_TIME_MONITOR("panzer::CubeHexMeshFactory::completeMeshConstruction()");
150 
151  if(not mesh.isInitialized())
152  mesh.initialize(parallelMach);
153 
154  // add node and element information
155  buildElements(parallelMach,mesh);
156 
157  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
158  out.setOutputToRootOnly(0);
159  out.setShowProcRank(true);
160 
161  // finish up the edges and faces
162  if(buildSubcells_) {
163  mesh.buildSubcells();
164 
165  out << "CubeHexMesh: Building sub cells" << std::endl;
166  }
167  else {
168  addSides(mesh);
169 
170  out << "CubeHexMesh: NOT building sub cells" << std::endl;
171  }
172 
173  mesh.buildLocalElementIDs();
174  if(createEdgeBlocks_) {
175  mesh.buildLocalEdgeIDs();
176  }
177  if(createFaceBlocks_) {
178  mesh.buildLocalFaceIDs();
179  }
180 
181  mesh.beginModification();
182 
183  // now that edges are built, side and node sets can be added
184  addSideSets(mesh);
185  addNodeSets(mesh);
186  if(createEdgeBlocks_) {
187  addEdgeBlocks(mesh);
188  }
189  if(createFaceBlocks_) {
190  addFaceBlocks(mesh);
191  }
192 
193  mesh.endModification();
194 
195  // calls Stk_MeshFactory::rebalance
196  this->rebalance(mesh);
197 }
198 
200 void CubeHexMeshFactory::setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & paramList)
201 {
202  paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
203 
204  setMyParamList(paramList);
205 
206  x0_ = paramList->get<double>("X0");
207  y0_ = paramList->get<double>("Y0");
208  z0_ = paramList->get<double>("Z0");
209 
210  xf_ = paramList->get<double>("Xf");
211  yf_ = paramList->get<double>("Yf");
212  zf_ = paramList->get<double>("Zf");
213 
214  xBlocks_ = paramList->get<int>("X Blocks");
215  yBlocks_ = paramList->get<int>("Y Blocks");
216  zBlocks_ = paramList->get<int>("Z Blocks");
217 
218  xProcs_ = paramList->get<int>("X Procs");
219  yProcs_ = paramList->get<int>("Y Procs");
220  zProcs_ = paramList->get<int>("Z Procs");
221 
222  nXElems_ = paramList->get<int>("X Elements");
223  nYElems_ = paramList->get<int>("Y Elements");
224  nZElems_ = paramList->get<int>("Z Elements");
225 
226  buildInterfaceSidesets_ = paramList->get<bool>("Build Interface Sidesets");
227 
228  buildSubcells_ = paramList->get<bool>("Build Subcells");
229 
230  createEdgeBlocks_ = paramList->get<bool>("Create Edge Blocks");
231  createFaceBlocks_ = paramList->get<bool>("Create Face Blocks");
232  if (not buildSubcells_ && createEdgeBlocks_) {
233  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
234  out.setOutputToRootOnly(0);
235  out.setShowProcRank(true);
236 
237  out << "CubeHexMesh: NOT creating edge blocks because building sub cells disabled" << std::endl;
238  createEdgeBlocks_ = false;
239  }
240  if (not buildSubcells_ && createFaceBlocks_) {
241  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
242  out.setOutputToRootOnly(0);
243  out.setShowProcRank(true);
244 
245  out << "CubeHexMesh: NOT creating face blocks because building sub cells disabled" << std::endl;
246  createFaceBlocks_ = false;
247  }
248 
249  // read in periodic boundary conditions
250  parsePeriodicBCList(Teuchos::rcpFromRef(paramList->sublist("Periodic BCs")),periodicBCVec_);
251 }
252 
254 Teuchos::RCP<const Teuchos::ParameterList> CubeHexMeshFactory::getValidParameters() const
255 {
256  static RCP<Teuchos::ParameterList> defaultParams;
257 
258  // fill with default values
259  if(defaultParams == Teuchos::null) {
260  defaultParams = rcp(new Teuchos::ParameterList);
261 
262  defaultParams->set<double>("X0",0.0);
263  defaultParams->set<double>("Y0",0.0);
264  defaultParams->set<double>("Z0",0.0);
265 
266  defaultParams->set<double>("Xf",1.0);
267  defaultParams->set<double>("Yf",1.0);
268  defaultParams->set<double>("Zf",1.0);
269 
270  defaultParams->set<int>("X Blocks",1);
271  defaultParams->set<int>("Y Blocks",1);
272  defaultParams->set<int>("Z Blocks",1);
273 
274  defaultParams->set<int>("X Procs",-1);
275  defaultParams->set<int>("Y Procs",1);
276  defaultParams->set<int>("Z Procs",1);
277 
278  defaultParams->set<int>("X Elements",5);
279  defaultParams->set<int>("Y Elements",5);
280  defaultParams->set<int>("Z Elements",5);
281 
282  defaultParams->set<bool>("Build Interface Sidesets",false);
283 
284  defaultParams->set<bool>("Build Subcells",true);
285 
286  // default to false for backward compatibility
287  defaultParams->set<bool>("Create Edge Blocks",false,"Create edge blocks in the mesh");
288  defaultParams->set<bool>("Create Face Blocks",false,"Create face blocks in the mesh");
289 
290  Teuchos::ParameterList & bcs = defaultParams->sublist("Periodic BCs");
291  bcs.set<int>("Count",0); // no default periodic boundary conditions
292  }
293 
294  return defaultParams;
295 }
296 
298 {
299  // get valid parameters
300  RCP<Teuchos::ParameterList> validParams = rcp(new Teuchos::ParameterList(*getValidParameters()));
301 
302  // set that parameter list
303  setParameterList(validParams);
304 }
305 
306 void CubeHexMeshFactory::buildMetaData(stk::ParallelMachine /* parallelMach */, STK_Interface & mesh) const
307 {
308  typedef shards::Hexahedron<8> HexTopo;
309  const CellTopologyData * ctd = shards::getCellTopologyData<HexTopo>();
310  const CellTopologyData * side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(2,0);
311 
312  // build meta data
313  //mesh.setDimension(2);
314  for(int bx=0;bx<xBlocks_;bx++) {
315  for(int by=0;by<yBlocks_;by++) {
316  for(int bz=0;bz<zBlocks_;bz++) {
317 
318  std::stringstream ebPostfix;
319  ebPostfix << "-" << bx << "_" << by << "_" << bz;
320 
321  // add element blocks
322  mesh.addElementBlock("eblock"+ebPostfix.str(),ctd);
323  }
324  }
325  }
326 
327  // add sidesets
328  mesh.addSideset("left",side_ctd);
329  mesh.addSideset("right",side_ctd);
330  mesh.addSideset("top",side_ctd);
331  mesh.addSideset("bottom",side_ctd);
332  mesh.addSideset("front",side_ctd);
333  mesh.addSideset("back",side_ctd);
334 
336  for(int bx=1;bx<xBlocks_;bx++) {
337  std::stringstream ss;
338  ss << "vertical_" << bx-1;
339  mesh.addSideset(ss.str(),side_ctd);
340  }
341  for(int by=1;by<yBlocks_;by++) {
342  std::stringstream ss;
343  ss << "horizontal_" << by-1;
344  mesh.addSideset(ss.str(),side_ctd);
345  }
346  for(int bz=1;bz<zBlocks_;bz++) {
347  std::stringstream ss;
348  ss << "transverse_" << bz-1;
349  mesh.addSideset(ss.str(),side_ctd);
350  }
351  }
352 
353  // add nodesets
354  mesh.addNodeset("origin");
355 
356  if(createEdgeBlocks_) {
357  const CellTopologyData * edge_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(1,0);
359  }
360  if(createFaceBlocks_) {
361  const CellTopologyData * face_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(2,0);
363  }
364 }
365 
366 void CubeHexMeshFactory::buildElements(stk::ParallelMachine parallelMach,STK_Interface & mesh) const
367 {
368  mesh.beginModification();
369  // build each block
370  for(int xBlock=0;xBlock<xBlocks_;xBlock++) {
371  for(int yBlock=0;yBlock<yBlocks_;yBlock++) {
372  for(int zBlock=0;zBlock<zBlocks_;zBlock++) {
373  buildBlock(parallelMach,xBlock,yBlock,zBlock,mesh);
374  }
375  }
376  }
377  mesh.endModification();
378 }
379 
380 void CubeHexMeshFactory::buildBlock(stk::ParallelMachine /* parallelMach */,int xBlock,int yBlock,int zBlock,STK_Interface & mesh) const
381 {
382  // grab this processors rank and machine size
383  std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> sizeAndStartX = determineXElemSizeAndStart(xBlock,xProcs_,machRank_);
384  std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> sizeAndStartY = determineYElemSizeAndStart(yBlock,yProcs_,machRank_);
385  std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> sizeAndStartZ = determineZElemSizeAndStart(zBlock,zProcs_,machRank_);
386 
387  panzer::GlobalOrdinal myXElems_start = sizeAndStartX.first;
388  panzer::GlobalOrdinal myXElems_end = myXElems_start+sizeAndStartX.second;
389  panzer::GlobalOrdinal myYElems_start = sizeAndStartY.first;
390  panzer::GlobalOrdinal myYElems_end = myYElems_start+sizeAndStartY.second;
391  panzer::GlobalOrdinal myZElems_start = sizeAndStartZ.first;
392  panzer::GlobalOrdinal myZElems_end = myZElems_start+sizeAndStartZ.second;
393 
394  panzer::GlobalOrdinal totalXElems = nXElems_*xBlocks_;
395  panzer::GlobalOrdinal totalYElems = nYElems_*yBlocks_;
396  panzer::GlobalOrdinal totalZElems = nZElems_*zBlocks_;
397 
398  double deltaX = (xf_-x0_)/double(totalXElems);
399  double deltaY = (yf_-y0_)/double(totalYElems);
400  double deltaZ = (zf_-z0_)/double(totalZElems);
401 
402  std::vector<double> coord(3,0.0);
403 
404  // build the nodes
405  for(panzer::GlobalOrdinal nx=myXElems_start;nx<myXElems_end+1;++nx) {
406  coord[0] = this->getMeshCoord(nx, deltaX, x0_);
407  for(panzer::GlobalOrdinal ny=myYElems_start;ny<myYElems_end+1;++ny) {
408  coord[1] = this->getMeshCoord(ny, deltaY, y0_);
409  for(panzer::GlobalOrdinal nz=myZElems_start;nz<myZElems_end+1;++nz) {
410  coord[2] = this->getMeshCoord(nz, deltaZ, z0_);
411 
412  mesh.addNode(nz*(totalYElems+1)*(totalXElems+1)+ny*(totalXElems+1)+nx+1,coord);
413  }
414  }
415  }
416 
417  std::stringstream blockName;
418  blockName << "eblock-" << xBlock << "_" << yBlock << "_" << zBlock;
419  stk::mesh::Part * block = mesh.getElementBlockPart(blockName.str());
420 
421  // build the elements
422  for(panzer::GlobalOrdinal nx=myXElems_start;nx<myXElems_end;++nx) {
423  for(panzer::GlobalOrdinal ny=myYElems_start;ny<myYElems_end;++ny) {
424  for(panzer::GlobalOrdinal nz=myZElems_start;nz<myZElems_end;++nz) {
425  stk::mesh::EntityId gid = totalXElems*totalYElems*nz+totalXElems*ny+nx+1;
426  std::vector<stk::mesh::EntityId> nodes(8);
427  nodes[0] = nx+1+ny*(totalXElems+1) +nz*(totalYElems+1)*(totalXElems+1);
428  nodes[1] = nodes[0]+1;
429  nodes[2] = nodes[1]+(totalXElems+1);
430  nodes[3] = nodes[2]-1;
431  nodes[4] = nodes[0]+(totalYElems+1)*(totalXElems+1);
432  nodes[5] = nodes[1]+(totalYElems+1)*(totalXElems+1);
433  nodes[6] = nodes[2]+(totalYElems+1)*(totalXElems+1);
434  nodes[7] = nodes[3]+(totalYElems+1)*(totalXElems+1);
435 
436  RCP<ElementDescriptor> ed = rcp(new ElementDescriptor(gid,nodes));
437  mesh.addElement(ed,block);
438  }
439  }
440  }
441 }
442 
443 std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> CubeHexMeshFactory::determineXElemSizeAndStart(int xBlock,unsigned int size,unsigned int /* rank */) const
444 {
445  std::size_t xProcLoc = procTuple_[0];
446  panzer::GlobalOrdinal minElements = nXElems_/size;
447  panzer::GlobalOrdinal extra = nXElems_ - minElements*size;
448 
449  TEUCHOS_ASSERT(minElements>0);
450 
451  // first "extra" elements get an extra column of elements
452  // this determines the starting X index and number of elements
453  panzer::GlobalOrdinal nume=0, start=0;
454  if(panzer::GlobalOrdinal(xProcLoc)<extra) {
455  nume = minElements+1;
456  start = xProcLoc*(minElements+1);
457  }
458  else {
459  nume = minElements;
460  start = extra*(minElements+1)+(xProcLoc-extra)*minElements;
461  }
462 
463  return std::make_pair(start+nXElems_*xBlock,nume);
464 }
465 
466 std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> CubeHexMeshFactory::determineYElemSizeAndStart(int yBlock,unsigned int size,unsigned int /* rank */) const
467 {
468  // int start = yBlock*nYElems_;
469  // return std::make_pair(start,nYElems_);
470 
471  std::size_t yProcLoc = procTuple_[1];
472  panzer::GlobalOrdinal minElements = nYElems_/size;
473  panzer::GlobalOrdinal extra = nYElems_ - minElements*size;
474 
475  TEUCHOS_ASSERT(minElements>0);
476 
477  // first "extra" elements get an extra column of elements
478  // this determines the starting X index and number of elements
479  panzer::GlobalOrdinal nume=0, start=0;
480  if(panzer::GlobalOrdinal(yProcLoc)<extra) {
481  nume = minElements+1;
482  start = yProcLoc*(minElements+1);
483  }
484  else {
485  nume = minElements;
486  start = extra*(minElements+1)+(yProcLoc-extra)*minElements;
487  }
488 
489  return std::make_pair(start+nYElems_*yBlock,nume);
490 }
491 
492 std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> CubeHexMeshFactory::determineZElemSizeAndStart(int zBlock,unsigned int size,unsigned int /* rank */) const
493 {
494  // int start = zBlock*nZElems_;
495  // return std::make_pair(start,nZElems_);
496  std::size_t zProcLoc = procTuple_[2];
497  panzer::GlobalOrdinal minElements = nZElems_/size;
498  panzer::GlobalOrdinal extra = nZElems_ - minElements*size;
499 
500  TEUCHOS_ASSERT(minElements>0);
501 
502  // first "extra" elements get an extra column of elements
503  // this determines the starting X index and number of elements
504  panzer::GlobalOrdinal nume=0, start=0;
505  if(zProcLoc<Teuchos::as<std::size_t>(extra)) {
506  nume = minElements+1;
507  start = zProcLoc*(minElements+1);
508  }
509  else {
510  nume = minElements;
511  start = extra*(minElements+1)+(zProcLoc-extra)*minElements;
512  }
513 
514  return std::make_pair(start+nZElems_*zBlock,nume);
515 }
516 
517 // this adds side entities only (does not inject them into side sets)
519 {
520  mesh.beginModification();
521 
522  std::size_t totalXElems = nXElems_*xBlocks_;
523  std::size_t totalYElems = nYElems_*yBlocks_;
524  std::size_t totalZElems = nZElems_*zBlocks_;
525 
526  std::vector<stk::mesh::Entity> localElmts;
527  mesh.getMyElements(localElmts);
528 
529  stk::mesh::EntityId offset[6];
530  offset[0] = 0;
531  offset[1] = offset[0] + totalXElems*totalZElems;
532  offset[2] = offset[1] + totalYElems*totalZElems;
533  offset[3] = offset[2] + totalXElems*totalZElems;
534  offset[4] = offset[3] + totalYElems*totalZElems;
535  offset[5] = offset[4] + totalXElems*totalYElems;
536 
537  // gid = totalXElems*totalYElems*nz+totalXElems*ny+nx+1
538 
539  // loop over elements adding sides to sidesets
540  std::vector<stk::mesh::Entity>::const_iterator itr;
541  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
542  stk::mesh::Entity element = (*itr);
543  stk::mesh::EntityId gid = mesh.elementGlobalId(element);
544 
545  std::size_t nx,ny,nz;
546  nz = (gid-1) / (totalXElems*totalYElems);
547  gid = (gid-1)-nz*(totalXElems*totalYElems);
548  ny = gid / totalXElems;
549  nx = gid-ny*totalXElems;
550 
551  std::vector<stk::mesh::Part*> parts;
552 
553  if(nz==0) {
554  // on the back
555  mesh.getBulkData()->declare_element_side(element, 4, parts);
556  }
557  if(nz+1==totalZElems) {
558  // on the front
559  mesh.getBulkData()->declare_element_side(element, 5, parts);
560  }
561 
562  if(ny==0) {
563  // on the bottom
564  mesh.getBulkData()->declare_element_side(element, 0, parts);
565  }
566  if(ny+1==totalYElems) {
567  // on the top
568  mesh.getBulkData()->declare_element_side(element, 2, parts);
569  }
570 
571  if(nx==0) {
572  // on the left
573  mesh.getBulkData()->declare_element_side(element, 3, parts);
574  }
575  if(nx+1==totalXElems) {
576  // on the right
577  mesh.getBulkData()->declare_element_side(element, 1, parts);
578  }
579  }
580 
581  mesh.endModification();
582 }
583 
584 // Pre-Condition: call beginModification() before entry
585 // Post-Condition: call endModification() after exit
587 {
588  const stk::mesh::EntityRank side_rank = mesh.getSideRank();
589 
590  std::size_t totalXElems = nXElems_*xBlocks_;
591  std::size_t totalYElems = nYElems_*yBlocks_;
592  std::size_t totalZElems = nZElems_*zBlocks_;
593 
594  // get all part vectors
595  stk::mesh::Part * left = mesh.getSideset("left");
596  stk::mesh::Part * right = mesh.getSideset("right");
597  stk::mesh::Part * top = mesh.getSideset("top");
598  stk::mesh::Part * bottom = mesh.getSideset("bottom");
599  stk::mesh::Part * front = mesh.getSideset("front");
600  stk::mesh::Part * back = mesh.getSideset("back");
601 
602  std::vector<stk::mesh::Part*> vertical;
603  std::vector<stk::mesh::Part*> horizontal;
604  std::vector<stk::mesh::Part*> transverse;
605 
607  for(int bx=1;bx<xBlocks_;bx++) {
608  std::stringstream ss;
609  ss << "vertical_" << bx-1;
610  vertical.push_back(mesh.getSideset(ss.str()));
611  }
612  for(int by=1;by<yBlocks_;by++) {
613  std::stringstream ss;
614  ss << "horizontal_" << by-1;
615  horizontal.push_back(mesh.getSideset(ss.str()));
616  }
617  for(int bz=1;bz<zBlocks_;bz++) {
618  std::stringstream ss;
619  ss << "transverse_" << bz-1;
620  transverse.push_back(mesh.getSideset(ss.str()));
621  }
622  }
623 
624  std::vector<stk::mesh::Entity> localElmts;
625  mesh.getMyElements(localElmts);
626 
627  // gid = totalXElems*totalYElems*nz+totalXElems*ny+nx+1
628 
629  // loop over elements adding sides to sidesets
630  std::vector<stk::mesh::Entity>::const_iterator itr;
631  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
632  stk::mesh::Entity element = (*itr);
633  stk::mesh::EntityId gid = mesh.elementGlobalId(element);
634 
635  std::size_t nx,ny,nz;
636  nz = (gid-1) / (totalXElems*totalYElems);
637  gid = (gid-1)-nz*(totalXElems*totalYElems);
638  ny = gid / totalXElems;
639  nx = gid-ny*totalXElems;
640 
641  if(nz % nZElems_==0) {
642  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 4);
643 
644  // on the back
645  if(mesh.entityOwnerRank(side)==machRank_) {
646  if(nz==0) {
647  mesh.addEntityToSideset(side,back);
648  } else {
650  int index = nz/nZElems_-1;
651  mesh.addEntityToSideset(side,transverse[index]);
652  }
653  }
654  }
655  }
656  if((nz+1) % nZElems_==0) {
657  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 5);
658 
659  // on the front
660  if(mesh.entityOwnerRank(side)==machRank_) {
661  if(nz+1==totalZElems) {
662  mesh.addEntityToSideset(side,front);
663  } else {
665  int index = (nz+1)/nZElems_-1;
666  mesh.addEntityToSideset(side,transverse[index]);
667  }
668  }
669  }
670  }
671 
672  if(ny % nYElems_==0) {
673  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 0);
674 
675  // on the bottom
676  if(mesh.entityOwnerRank(side)==machRank_) {
677  if(ny==0) {
678  mesh.addEntityToSideset(side,bottom);
679  } else {
681  int index = ny/nYElems_-1;
682  mesh.addEntityToSideset(side,horizontal[index]);
683  }
684  }
685  }
686  }
687  if((ny+1) % nYElems_==0) {
688  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 2);
689 
690  // on the top
691  if(mesh.entityOwnerRank(side)==machRank_) {
692  if(ny+1==totalYElems) {
693  mesh.addEntityToSideset(side,top);
694  } else {
696  int index = (ny+1)/nYElems_-1;
697  mesh.addEntityToSideset(side,horizontal[index]);
698  }
699  }
700  }
701  }
702 
703  if(nx % nXElems_==0) {
704  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 3);
705 
706  // on the left
707  if(mesh.entityOwnerRank(side)==machRank_) {
708  if(nx==0) {
709  mesh.addEntityToSideset(side,left);
710  } else {
712  int index = nx/nXElems_-1;
713  mesh.addEntityToSideset(side,vertical[index]);
714  }
715  }
716  }
717  }
718  if((nx+1) % nXElems_==0) {
719  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 1);
720 
721  // on the right
722  if(mesh.entityOwnerRank(side)==machRank_) {
723  if(nx+1==totalXElems) {
724  mesh.addEntityToSideset(side,right);
725  } else {
727  int index = (nx+1)/nXElems_-1;
728  mesh.addEntityToSideset(side,vertical[index]);
729  }
730  }
731  }
732  }
733  }
734 }
735 
736 // Pre-Condition: call beginModification() before entry
737 // Post-Condition: call endModification() after exit
739 {
740  // get all part vectors
741  stk::mesh::Part * origin = mesh.getNodeset("origin");
742 
743  Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
744  if(machRank_==0)
745  {
746  // add zero node to origin node set
747  stk::mesh::Entity node = bulkData->get_entity(mesh.getNodeRank(),1);
748  mesh.addEntityToNodeset(node,origin);
749  }
750 }
751 
752 // Pre-Condition: call beginModification() before entry
753 // Post-Condition: call endModification() after exit
755 {
756  stk::mesh::Part * edge_block = mesh.getEdgeBlock(panzer_stk::STK_Interface::edgeBlockString);
757 
758  Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
759  Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
760 
761  std::vector<stk::mesh::Entity> edges;
762  bulkData->get_entities(mesh.getEdgeRank(),metaData->locally_owned_part(),edges);
763  mesh.addEntitiesToEdgeBlock(edges, edge_block);
764 }
765 
766 // Pre-Condition: call beginModification() before entry
767 // Post-Condition: call endModification() after exit
769 {
770  stk::mesh::Part * face_block = mesh.getFaceBlock(panzer_stk::STK_Interface::faceBlockString);
771 
772  Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
773  Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
774 
775  std::vector<stk::mesh::Entity> faces;
776  bulkData->get_entities(mesh.getFaceRank(),metaData->locally_owned_part(),faces);
777  mesh.addEntitiesToFaceBlock(faces, face_block);
778 }
779 
781 Teuchos::Tuple<std::size_t,3> CubeHexMeshFactory::procRankToProcTuple(std::size_t procRank) const
782 {
783  std::size_t i=0,j=0,k=0;
784 
785  k = procRank/(xProcs_*yProcs_); procRank = procRank % (xProcs_*yProcs_);
786  j = procRank/xProcs_; procRank = procRank % xProcs_;
787  i = procRank;
788 
789  return Teuchos::tuple(i,j,k);
790 }
791 
792 } // end panzer_stk
std::pair< panzer::GlobalOrdinal, panzer::GlobalOrdinal > determineXElemSizeAndStart(int xBlock, unsigned int size, unsigned int rank) const
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 addSides(STK_Interface &mesh) const
void addEntityToNodeset(stk::mesh::Entity entity, stk::mesh::Part *nodeset)
stk::mesh::EntityRank getFaceRank() const
virtual Teuchos::RCP< STK_Interface > buildUncommitedMesh(stk::ParallelMachine parallelMach) const
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
void addSideSets(STK_Interface &mesh) const
virtual void completeMeshConstruction(STK_Interface &mesh, stk::ParallelMachine parallelMach) const
stk::mesh::EntityRank getSideRank() const
stk::mesh::Entity findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const
stk::mesh::Part * getEdgeBlock(const std::string &name) const
get the block part
std::pair< panzer::GlobalOrdinal, panzer::GlobalOrdinal > determineZElemSizeAndStart(int zBlock, unsigned int size, unsigned int rank) const
stk::mesh::EntityRank getEdgeRank() const
stk::mesh::Part * getSideset(const std::string &name) const
void addFaceBlocks(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 buildMetaData(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true, const bool buildRefinementSupport=false)
Teuchos::RCP< STK_Interface > buildMesh(stk::ParallelMachine parallelMach) const
Build the mesh object.
double getMeshCoord(const int nx, const double deltaX, const double x0) const
stk::mesh::Part * getFaceBlock(const std::string &name) const
get the block part
bool isInitialized() const
Has initialize been called on this mesh object?
static const std::string edgeBlockString
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
From ParameterListAcceptor.
Teuchos::Tuple< std::size_t, 3 > procRankToProcTuple(std::size_t procRank) const
what is the 3D tuple describe this processor distribution
void addFaceBlock(const std::string &name, const CellTopologyData *ctData)
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
Teuchos::Tuple< std::size_t, 3 > procTuple_
void buildElements(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
void addNodeSets(STK_Interface &mesh) const
void addSideset(const std::string &name, const CellTopologyData *ctData)
static const std::string faceBlockString
void addEdgeBlock(const std::string &name, const CellTopologyData *ctData)
std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > periodicBCVec_
void addEntitiesToFaceBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *faceblock)
unsigned entityOwnerRank(stk::mesh::Entity entity) const
void rebalance(STK_Interface &mesh) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
From ParameterListAcceptor.
std::pair< panzer::GlobalOrdinal, panzer::GlobalOrdinal > determineYElemSizeAndStart(int yBlock, unsigned int size, unsigned int rank) const
void addEdgeBlocks(STK_Interface &mesh) const
void addElementBlock(const std::string &name, const CellTopologyData *ctData)
stk::mesh::Part * getNodeset(const std::string &name) const
void buildBlock(stk::ParallelMachine machRank, int xBlock, int yBlock, int zBlock, STK_Interface &mesh) const
static void parsePeriodicBCList(const Teuchos::RCP< Teuchos::ParameterList > &pl, std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > &periodicBC)