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 <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 
175  // now that edges are built, side and node sets can be added
176  addSideSets(mesh);
177  addNodeSets(mesh);
178 
179  // calls Stk_MeshFactory::rebalance
180  this->rebalance(mesh);
181 }
182 
184 void CubeHexMeshFactory::setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & paramList)
185 {
186  paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
187 
188  setMyParamList(paramList);
189 
190  x0_ = paramList->get<double>("X0");
191  y0_ = paramList->get<double>("Y0");
192  z0_ = paramList->get<double>("Z0");
193 
194  xf_ = paramList->get<double>("Xf");
195  yf_ = paramList->get<double>("Yf");
196  zf_ = paramList->get<double>("Zf");
197 
198  xBlocks_ = paramList->get<int>("X Blocks");
199  yBlocks_ = paramList->get<int>("Y Blocks");
200  zBlocks_ = paramList->get<int>("Z Blocks");
201 
202  xProcs_ = paramList->get<int>("X Procs");
203  yProcs_ = paramList->get<int>("Y Procs");
204  zProcs_ = paramList->get<int>("Z Procs");
205 
206  nXElems_ = paramList->get<int>("X Elements");
207  nYElems_ = paramList->get<int>("Y Elements");
208  nZElems_ = paramList->get<int>("Z Elements");
209 
210  buildInterfaceSidesets_ = paramList->get<bool>("Build Interface Sidesets");
211 
212  buildSubcells_ = paramList->get<bool>("Build Subcells");
213 
214  // read in periodic boundary conditions
215  parsePeriodicBCList(Teuchos::rcpFromRef(paramList->sublist("Periodic BCs")),periodicBCVec_);
216 }
217 
219 Teuchos::RCP<const Teuchos::ParameterList> CubeHexMeshFactory::getValidParameters() const
220 {
221  static RCP<Teuchos::ParameterList> defaultParams;
222 
223  // fill with default values
224  if(defaultParams == Teuchos::null) {
225  defaultParams = rcp(new Teuchos::ParameterList);
226 
227  defaultParams->set<double>("X0",0.0);
228  defaultParams->set<double>("Y0",0.0);
229  defaultParams->set<double>("Z0",0.0);
230 
231  defaultParams->set<double>("Xf",1.0);
232  defaultParams->set<double>("Yf",1.0);
233  defaultParams->set<double>("Zf",1.0);
234 
235  defaultParams->set<int>("X Blocks",1);
236  defaultParams->set<int>("Y Blocks",1);
237  defaultParams->set<int>("Z Blocks",1);
238 
239  defaultParams->set<int>("X Procs",-1);
240  defaultParams->set<int>("Y Procs",1);
241  defaultParams->set<int>("Z Procs",1);
242 
243  defaultParams->set<int>("X Elements",5);
244  defaultParams->set<int>("Y Elements",5);
245  defaultParams->set<int>("Z Elements",5);
246 
247  defaultParams->set<bool>("Build Interface Sidesets",false);
248 
249  defaultParams->set<bool>("Build Subcells",true);
250 
251  Teuchos::ParameterList & bcs = defaultParams->sublist("Periodic BCs");
252  bcs.set<int>("Count",0); // no default periodic boundary conditions
253  }
254 
255  return defaultParams;
256 }
257 
259 {
260  // get valid parameters
261  RCP<Teuchos::ParameterList> validParams = rcp(new Teuchos::ParameterList(*getValidParameters()));
262 
263  // set that parameter list
264  setParameterList(validParams);
265 }
266 
267 void CubeHexMeshFactory::buildMetaData(stk::ParallelMachine /* parallelMach */, STK_Interface & mesh) const
268 {
269  typedef shards::Hexahedron<8> HexTopo;
270  const CellTopologyData * ctd = shards::getCellTopologyData<HexTopo>();
271  const CellTopologyData * side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(2,0);
272 
273  // build meta data
274  //mesh.setDimension(2);
275  for(int bx=0;bx<xBlocks_;bx++) {
276  for(int by=0;by<yBlocks_;by++) {
277  for(int bz=0;bz<zBlocks_;bz++) {
278 
279  std::stringstream ebPostfix;
280  ebPostfix << "-" << bx << "_" << by << "_" << bz;
281 
282  // add element blocks
283  mesh.addElementBlock("eblock"+ebPostfix.str(),ctd);
284  }
285  }
286  }
287 
288  // add sidesets
289  mesh.addSideset("left",side_ctd);
290  mesh.addSideset("right",side_ctd);
291  mesh.addSideset("top",side_ctd);
292  mesh.addSideset("bottom",side_ctd);
293  mesh.addSideset("front",side_ctd);
294  mesh.addSideset("back",side_ctd);
295 
297  for(int bx=1;bx<xBlocks_;bx++) {
298  std::stringstream ss;
299  ss << "vertical_" << bx-1;
300  mesh.addSideset(ss.str(),side_ctd);
301  }
302  for(int by=1;by<yBlocks_;by++) {
303  std::stringstream ss;
304  ss << "horizontal_" << by-1;
305  mesh.addSideset(ss.str(),side_ctd);
306  }
307  for(int bz=1;bz<zBlocks_;bz++) {
308  std::stringstream ss;
309  ss << "transverse_" << bz-1;
310  mesh.addSideset(ss.str(),side_ctd);
311  }
312  }
313 
314  // add nodesets
315  mesh.addNodeset("origin");
316 }
317 
318 void CubeHexMeshFactory::buildElements(stk::ParallelMachine parallelMach,STK_Interface & mesh) const
319 {
320  mesh.beginModification();
321  // build each block
322  for(int xBlock=0;xBlock<xBlocks_;xBlock++) {
323  for(int yBlock=0;yBlock<yBlocks_;yBlock++) {
324  for(int zBlock=0;zBlock<zBlocks_;zBlock++) {
325  buildBlock(parallelMach,xBlock,yBlock,zBlock,mesh);
326  }
327  }
328  }
329  mesh.endModification();
330 }
331 
332 void CubeHexMeshFactory::buildBlock(stk::ParallelMachine /* parallelMach */,int xBlock,int yBlock,int zBlock,STK_Interface & mesh) const
333 {
334  // grab this processors rank and machine size
335  std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> sizeAndStartX = determineXElemSizeAndStart(xBlock,xProcs_,machRank_);
336  std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> sizeAndStartY = determineYElemSizeAndStart(yBlock,yProcs_,machRank_);
337  std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> sizeAndStartZ = determineZElemSizeAndStart(zBlock,zProcs_,machRank_);
338 
339  panzer::GlobalOrdinal myXElems_start = sizeAndStartX.first;
340  panzer::GlobalOrdinal myXElems_end = myXElems_start+sizeAndStartX.second;
341  panzer::GlobalOrdinal myYElems_start = sizeAndStartY.first;
342  panzer::GlobalOrdinal myYElems_end = myYElems_start+sizeAndStartY.second;
343  panzer::GlobalOrdinal myZElems_start = sizeAndStartZ.first;
344  panzer::GlobalOrdinal myZElems_end = myZElems_start+sizeAndStartZ.second;
345 
346  panzer::GlobalOrdinal totalXElems = nXElems_*xBlocks_;
347  panzer::GlobalOrdinal totalYElems = nYElems_*yBlocks_;
348  panzer::GlobalOrdinal totalZElems = nZElems_*zBlocks_;
349 
350  double deltaX = (xf_-x0_)/double(totalXElems);
351  double deltaY = (yf_-y0_)/double(totalYElems);
352  double deltaZ = (zf_-z0_)/double(totalZElems);
353 
354  std::vector<double> coord(3,0.0);
355 
356  // build the nodes
357  for(panzer::GlobalOrdinal nx=myXElems_start;nx<myXElems_end+1;++nx) {
358  coord[0] = this->getMeshCoord(nx, deltaX, x0_);
359  for(panzer::GlobalOrdinal ny=myYElems_start;ny<myYElems_end+1;++ny) {
360  coord[1] = this->getMeshCoord(ny, deltaY, y0_);
361  for(panzer::GlobalOrdinal nz=myZElems_start;nz<myZElems_end+1;++nz) {
362  coord[2] = this->getMeshCoord(nz, deltaZ, z0_);
363 
364  mesh.addNode(nz*(totalYElems+1)*(totalXElems+1)+ny*(totalXElems+1)+nx+1,coord);
365  }
366  }
367  }
368 
369  std::stringstream blockName;
370  blockName << "eblock-" << xBlock << "_" << yBlock << "_" << zBlock;
371  stk::mesh::Part * block = mesh.getElementBlockPart(blockName.str());
372 
373  // build the elements
374  for(panzer::GlobalOrdinal nx=myXElems_start;nx<myXElems_end;++nx) {
375  for(panzer::GlobalOrdinal ny=myYElems_start;ny<myYElems_end;++ny) {
376  for(panzer::GlobalOrdinal nz=myZElems_start;nz<myZElems_end;++nz) {
377  stk::mesh::EntityId gid = totalXElems*totalYElems*nz+totalXElems*ny+nx+1;
378  std::vector<stk::mesh::EntityId> nodes(8);
379  nodes[0] = nx+1+ny*(totalXElems+1) +nz*(totalYElems+1)*(totalXElems+1);
380  nodes[1] = nodes[0]+1;
381  nodes[2] = nodes[1]+(totalXElems+1);
382  nodes[3] = nodes[2]-1;
383  nodes[4] = nodes[0]+(totalYElems+1)*(totalXElems+1);
384  nodes[5] = nodes[1]+(totalYElems+1)*(totalXElems+1);
385  nodes[6] = nodes[2]+(totalYElems+1)*(totalXElems+1);
386  nodes[7] = nodes[3]+(totalYElems+1)*(totalXElems+1);
387 
388  RCP<ElementDescriptor> ed = rcp(new ElementDescriptor(gid,nodes));
389  mesh.addElement(ed,block);
390  }
391  }
392  }
393 }
394 
395 std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> CubeHexMeshFactory::determineXElemSizeAndStart(int xBlock,unsigned int size,unsigned int /* rank */) const
396 {
397  std::size_t xProcLoc = procTuple_[0];
398  panzer::GlobalOrdinal minElements = nXElems_/size;
399  panzer::GlobalOrdinal extra = nXElems_ - minElements*size;
400 
401  TEUCHOS_ASSERT(minElements>0);
402 
403  // first "extra" elements get an extra column of elements
404  // this determines the starting X index and number of elements
405  panzer::GlobalOrdinal nume=0, start=0;
406  if(panzer::GlobalOrdinal(xProcLoc)<extra) {
407  nume = minElements+1;
408  start = xProcLoc*(minElements+1);
409  }
410  else {
411  nume = minElements;
412  start = extra*(minElements+1)+(xProcLoc-extra)*minElements;
413  }
414 
415  return std::make_pair(start+nXElems_*xBlock,nume);
416 }
417 
418 std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> CubeHexMeshFactory::determineYElemSizeAndStart(int yBlock,unsigned int size,unsigned int /* rank */) const
419 {
420  // int start = yBlock*nYElems_;
421  // return std::make_pair(start,nYElems_);
422 
423  std::size_t yProcLoc = procTuple_[1];
424  panzer::GlobalOrdinal minElements = nYElems_/size;
425  panzer::GlobalOrdinal extra = nYElems_ - minElements*size;
426 
427  TEUCHOS_ASSERT(minElements>0);
428 
429  // first "extra" elements get an extra column of elements
430  // this determines the starting X index and number of elements
431  panzer::GlobalOrdinal nume=0, start=0;
432  if(panzer::GlobalOrdinal(yProcLoc)<extra) {
433  nume = minElements+1;
434  start = yProcLoc*(minElements+1);
435  }
436  else {
437  nume = minElements;
438  start = extra*(minElements+1)+(yProcLoc-extra)*minElements;
439  }
440 
441  return std::make_pair(start+nYElems_*yBlock,nume);
442 }
443 
444 std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> CubeHexMeshFactory::determineZElemSizeAndStart(int zBlock,unsigned int size,unsigned int /* rank */) const
445 {
446  // int start = zBlock*nZElems_;
447  // return std::make_pair(start,nZElems_);
448  std::size_t zProcLoc = procTuple_[2];
449  panzer::GlobalOrdinal minElements = nZElems_/size;
450  panzer::GlobalOrdinal extra = nZElems_ - minElements*size;
451 
452  TEUCHOS_ASSERT(minElements>0);
453 
454  // first "extra" elements get an extra column of elements
455  // this determines the starting X index and number of elements
456  panzer::GlobalOrdinal nume=0, start=0;
457  if(zProcLoc<Teuchos::as<std::size_t>(extra)) {
458  nume = minElements+1;
459  start = zProcLoc*(minElements+1);
460  }
461  else {
462  nume = minElements;
463  start = extra*(minElements+1)+(zProcLoc-extra)*minElements;
464  }
465 
466  return std::make_pair(start+nZElems_*zBlock,nume);
467 }
468 
469 // this adds side entities only (does not inject them into side sets)
471 {
472  mesh.beginModification();
473 
474  std::size_t totalXElems = nXElems_*xBlocks_;
475  std::size_t totalYElems = nYElems_*yBlocks_;
476  std::size_t totalZElems = nZElems_*zBlocks_;
477 
478  std::vector<stk::mesh::Entity> localElmts;
479  mesh.getMyElements(localElmts);
480 
481  stk::mesh::EntityId offset[6];
482  offset[0] = 0;
483  offset[1] = offset[0] + totalXElems*totalZElems;
484  offset[2] = offset[1] + totalYElems*totalZElems;
485  offset[3] = offset[2] + totalXElems*totalZElems;
486  offset[4] = offset[3] + totalYElems*totalZElems;
487  offset[5] = offset[4] + totalXElems*totalYElems;
488 
489  // gid = totalXElems*totalYElems*nz+totalXElems*ny+nx+1
490 
491  // loop over elements adding sides to sidesets
492  std::vector<stk::mesh::Entity>::const_iterator itr;
493  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
494  stk::mesh::Entity element = (*itr);
495  stk::mesh::EntityId gid = mesh.elementGlobalId(element);
496 
497  std::size_t nx,ny,nz;
498  nz = (gid-1) / (totalXElems*totalYElems);
499  gid = (gid-1)-nz*(totalXElems*totalYElems);
500  ny = gid / totalXElems;
501  nx = gid-ny*totalXElems;
502 
503  std::vector<stk::mesh::Part*> parts;
504 
505  if(nz==0) {
506  // on the back
507  mesh.getBulkData()->declare_element_side(element, 4, parts);
508  }
509  if(nz+1==totalZElems) {
510  // on the front
511  mesh.getBulkData()->declare_element_side(element, 5, parts);
512  }
513 
514  if(ny==0) {
515  // on the bottom
516  mesh.getBulkData()->declare_element_side(element, 0, parts);
517  }
518  if(ny+1==totalYElems) {
519  // on the top
520  mesh.getBulkData()->declare_element_side(element, 2, parts);
521  }
522 
523  if(nx==0) {
524  // on the left
525  mesh.getBulkData()->declare_element_side(element, 3, parts);
526  }
527  if(nx+1==totalXElems) {
528  // on the right
529  mesh.getBulkData()->declare_element_side(element, 1, parts);
530  }
531  }
532 
533  mesh.endModification();
534 }
535 
537 {
538  mesh.beginModification();
539  const stk::mesh::EntityRank side_rank = mesh.getSideRank();
540 
541  std::size_t totalXElems = nXElems_*xBlocks_;
542  std::size_t totalYElems = nYElems_*yBlocks_;
543  std::size_t totalZElems = nZElems_*zBlocks_;
544 
545  // get all part vectors
546  stk::mesh::Part * left = mesh.getSideset("left");
547  stk::mesh::Part * right = mesh.getSideset("right");
548  stk::mesh::Part * top = mesh.getSideset("top");
549  stk::mesh::Part * bottom = mesh.getSideset("bottom");
550  stk::mesh::Part * front = mesh.getSideset("front");
551  stk::mesh::Part * back = mesh.getSideset("back");
552 
553  std::vector<stk::mesh::Part*> vertical;
554  std::vector<stk::mesh::Part*> horizontal;
555  std::vector<stk::mesh::Part*> transverse;
556 
558  for(int bx=1;bx<xBlocks_;bx++) {
559  std::stringstream ss;
560  ss << "vertical_" << bx-1;
561  vertical.push_back(mesh.getSideset(ss.str()));
562  }
563  for(int by=1;by<yBlocks_;by++) {
564  std::stringstream ss;
565  ss << "horizontal_" << by-1;
566  horizontal.push_back(mesh.getSideset(ss.str()));
567  }
568  for(int bz=1;bz<zBlocks_;bz++) {
569  std::stringstream ss;
570  ss << "transverse_" << bz-1;
571  transverse.push_back(mesh.getSideset(ss.str()));
572  }
573  }
574 
575  std::vector<stk::mesh::Entity> localElmts;
576  mesh.getMyElements(localElmts);
577 
578  // gid = totalXElems*totalYElems*nz+totalXElems*ny+nx+1
579 
580  // loop over elements adding sides to sidesets
581  std::vector<stk::mesh::Entity>::const_iterator itr;
582  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
583  stk::mesh::Entity element = (*itr);
584  stk::mesh::EntityId gid = mesh.elementGlobalId(element);
585 
586  std::size_t nx,ny,nz;
587  nz = (gid-1) / (totalXElems*totalYElems);
588  gid = (gid-1)-nz*(totalXElems*totalYElems);
589  ny = gid / totalXElems;
590  nx = gid-ny*totalXElems;
591 
592  if(nz % nZElems_==0) {
593  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 4);
594 
595  // on the back
596  if(mesh.entityOwnerRank(side)==machRank_) {
597  if(nz==0) {
598  mesh.addEntityToSideset(side,back);
599  } else {
601  int index = nz/nZElems_-1;
602  mesh.addEntityToSideset(side,transverse[index]);
603  }
604  }
605  }
606  }
607  if((nz+1) % nZElems_==0) {
608  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 5);
609 
610  // on the front
611  if(mesh.entityOwnerRank(side)==machRank_) {
612  if(nz+1==totalZElems) {
613  mesh.addEntityToSideset(side,front);
614  } else {
616  int index = (nz+1)/nZElems_-1;
617  mesh.addEntityToSideset(side,transverse[index]);
618  }
619  }
620  }
621  }
622 
623  if(ny % nYElems_==0) {
624  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 0);
625 
626  // on the bottom
627  if(mesh.entityOwnerRank(side)==machRank_) {
628  if(ny==0) {
629  mesh.addEntityToSideset(side,bottom);
630  } else {
632  int index = ny/nYElems_-1;
633  mesh.addEntityToSideset(side,horizontal[index]);
634  }
635  }
636  }
637  }
638  if((ny+1) % nYElems_==0) {
639  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 2);
640 
641  // on the top
642  if(mesh.entityOwnerRank(side)==machRank_) {
643  if(ny+1==totalYElems) {
644  mesh.addEntityToSideset(side,top);
645  } else {
647  int index = (ny+1)/nYElems_-1;
648  mesh.addEntityToSideset(side,horizontal[index]);
649  }
650  }
651  }
652  }
653 
654  if(nx % nXElems_==0) {
655  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 3);
656 
657  // on the left
658  if(mesh.entityOwnerRank(side)==machRank_) {
659  if(nx==0) {
660  mesh.addEntityToSideset(side,left);
661  } else {
663  int index = nx/nXElems_-1;
664  mesh.addEntityToSideset(side,vertical[index]);
665  }
666  }
667  }
668  }
669  if((nx+1) % nXElems_==0) {
670  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 1);
671 
672  // on the right
673  if(mesh.entityOwnerRank(side)==machRank_) {
674  if(nx+1==totalXElems) {
675  mesh.addEntityToSideset(side,right);
676  } else {
678  int index = (nx+1)/nXElems_-1;
679  mesh.addEntityToSideset(side,vertical[index]);
680  }
681  }
682  }
683  }
684  }
685 
686  mesh.endModification();
687 }
688 
690 {
691  mesh.beginModification();
692 
693  // get all part vectors
694  stk::mesh::Part * origin = mesh.getNodeset("origin");
695 
696  Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
697  if(machRank_==0)
698  {
699  // add zero node to origin node set
700  stk::mesh::Entity node = bulkData->get_entity(mesh.getNodeRank(),1);
701  mesh.addEntityToNodeset(node,origin);
702  }
703 
704  mesh.endModification();
705 }
706 
708 Teuchos::Tuple<std::size_t,3> CubeHexMeshFactory::procRankToProcTuple(std::size_t procRank) const
709 {
710  std::size_t i=0,j=0,k=0;
711 
712  k = procRank/(xProcs_*yProcs_); procRank = procRank % (xProcs_*yProcs_);
713  j = procRank/xProcs_; procRank = procRank % xProcs_;
714  i = procRank;
715 
716  return Teuchos::tuple(i,j,k);
717 }
718 
719 } // end panzer_stk
std::pair< panzer::GlobalOrdinal, panzer::GlobalOrdinal > determineXElemSizeAndStart(int xBlock, unsigned int size, unsigned int rank) const
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block count
void addNodeset(const std::string &name)
void addSides(STK_Interface &mesh) const
void addEntityToNodeset(stk::mesh::Entity entity, stk::mesh::Part *nodeset)
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
std::pair< panzer::GlobalOrdinal, panzer::GlobalOrdinal > determineZElemSizeAndStart(int zBlock, unsigned int size, unsigned int rank) 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
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
bool isInitialized() const
Has initialize been called on this mesh object?
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 addNode(stk::mesh::EntityId gid, const std::vector< double > &coord)
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)
std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > periodicBCVec_
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 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)