Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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  }
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  mesh->setBoundingBoxSearchFlag(useBBoxSearch_);
144 
145  return mesh;
146 }
147 
148 void CubeHexMeshFactory::completeMeshConstruction(STK_Interface & mesh,stk::ParallelMachine parallelMach) const
149 {
150  PANZER_FUNC_TIME_MONITOR("panzer::CubeHexMeshFactory::completeMeshConstruction()");
151 
152  if(not mesh.isInitialized())
153  mesh.initialize(parallelMach);
154 
155  // add node and element information
156  buildElements(parallelMach,mesh);
157 
158  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
159  out.setOutputToRootOnly(0);
160  out.setShowProcRank(true);
161 
162  // finish up the edges and faces
163  if(buildSubcells_) {
164  mesh.buildSubcells();
165 
166  out << "CubeHexMesh: Building sub cells" << std::endl;
167  }
168  else {
169  addSides(mesh);
170 
171  out << "CubeHexMesh: NOT building sub cells" << std::endl;
172  }
173 
174  mesh.buildLocalElementIDs();
175  if(createEdgeBlocks_) {
176  mesh.buildLocalEdgeIDs();
177  }
178  if(createFaceBlocks_) {
179  mesh.buildLocalFaceIDs();
180  }
181 
182  mesh.beginModification();
183 
184  // now that edges are built, side and node sets can be added
185  addSideSets(mesh);
186  addNodeSets(mesh);
187  if(createEdgeBlocks_) {
188  addEdgeBlocks(mesh);
189  }
190  if(createFaceBlocks_) {
191  addFaceBlocks(mesh);
192  }
193 
194  mesh.endModification();
195 
196  // calls Stk_MeshFactory::rebalance
197  this->rebalance(mesh);
198 }
199 
202 {
204 
205  setMyParamList(paramList);
206 
207  x0_ = paramList->get<double>("X0");
208  y0_ = paramList->get<double>("Y0");
209  z0_ = paramList->get<double>("Z0");
210 
211  xf_ = paramList->get<double>("Xf");
212  yf_ = paramList->get<double>("Yf");
213  zf_ = paramList->get<double>("Zf");
214 
215  xBlocks_ = paramList->get<int>("X Blocks");
216  yBlocks_ = paramList->get<int>("Y Blocks");
217  zBlocks_ = paramList->get<int>("Z Blocks");
218 
219  xProcs_ = paramList->get<int>("X Procs");
220  yProcs_ = paramList->get<int>("Y Procs");
221  zProcs_ = paramList->get<int>("Z Procs");
222 
223  nXElems_ = paramList->get<int>("X Elements");
224  nYElems_ = paramList->get<int>("Y Elements");
225  nZElems_ = paramList->get<int>("Z Elements");
226 
227  buildInterfaceSidesets_ = paramList->get<bool>("Build Interface Sidesets");
228 
229  buildSubcells_ = paramList->get<bool>("Build Subcells");
230 
231  createEdgeBlocks_ = paramList->get<bool>("Create Edge Blocks");
232  createFaceBlocks_ = paramList->get<bool>("Create Face Blocks");
233  if (not buildSubcells_ && createEdgeBlocks_) {
234  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
235  out.setOutputToRootOnly(0);
236  out.setShowProcRank(true);
237 
238  out << "CubeHexMesh: NOT creating edge blocks because building sub cells disabled" << std::endl;
239  createEdgeBlocks_ = false;
240  }
241  if (not buildSubcells_ && createFaceBlocks_) {
242  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
243  out.setOutputToRootOnly(0);
244  out.setShowProcRank(true);
245 
246  out << "CubeHexMesh: NOT creating face blocks because building sub cells disabled" << std::endl;
247  createFaceBlocks_ = false;
248  }
249 
250  // read in periodic boundary conditions
251  parsePeriodicBCList(Teuchos::rcpFromRef(paramList->sublist("Periodic BCs")),periodicBCVec_,useBBoxSearch_);
252 }
253 
256 {
257  static RCP<Teuchos::ParameterList> defaultParams;
258 
259  // fill with default values
260  if(defaultParams == Teuchos::null) {
261  defaultParams = rcp(new Teuchos::ParameterList);
262 
263  defaultParams->set<double>("X0",0.0);
264  defaultParams->set<double>("Y0",0.0);
265  defaultParams->set<double>("Z0",0.0);
266 
267  defaultParams->set<double>("Xf",1.0);
268  defaultParams->set<double>("Yf",1.0);
269  defaultParams->set<double>("Zf",1.0);
270 
271  defaultParams->set<int>("X Blocks",1);
272  defaultParams->set<int>("Y Blocks",1);
273  defaultParams->set<int>("Z Blocks",1);
274 
275  defaultParams->set<int>("X Procs",-1);
276  defaultParams->set<int>("Y Procs",1);
277  defaultParams->set<int>("Z Procs",1);
278 
279  defaultParams->set<int>("X Elements",5);
280  defaultParams->set<int>("Y Elements",5);
281  defaultParams->set<int>("Z Elements",5);
282 
283  defaultParams->set<bool>("Build Interface Sidesets",false);
284 
285  defaultParams->set<bool>("Build Subcells",true);
286 
287  // default to false for backward compatibility
288  defaultParams->set<bool>("Create Edge Blocks",false,"Create edge blocks in the mesh");
289  defaultParams->set<bool>("Create Face Blocks",false,"Create face blocks in the mesh");
290 
291  Teuchos::ParameterList & bcs = defaultParams->sublist("Periodic BCs");
292  bcs.set<int>("Count",0); // no default periodic boundary conditions
293  }
294 
295  return defaultParams;
296 }
297 
299 {
300  // get valid parameters
302 
303  // set that parameter list
304  setParameterList(validParams);
305 
306  /* This is a hex mesh factory so all elements in all element blocks
307  * will be hex8. This means that all the edges will be line2 and
308  * all the faces will be quad4. The edge and face block names are
309  * hard coded to reflect this.
310  */
313 }
314 
315 void CubeHexMeshFactory::buildMetaData(stk::ParallelMachine /* parallelMach */, STK_Interface & mesh) const
316 {
317  typedef shards::Hexahedron<8> HexTopo;
318  const CellTopologyData * ctd = shards::getCellTopologyData<HexTopo>();
319  const CellTopologyData * side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(2,0);
320  const CellTopologyData * edge_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(1,0);
321  const CellTopologyData * face_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(2,0);
322 
323  // build meta data
324  //mesh.setDimension(2);
325  for(int bx=0;bx<xBlocks_;bx++) {
326  for(int by=0;by<yBlocks_;by++) {
327  for(int bz=0;bz<zBlocks_;bz++) {
328 
329  std::stringstream ebPostfix;
330  ebPostfix << "-" << bx << "_" << by << "_" << bz;
331 
332  // add element blocks
333  mesh.addElementBlock("eblock"+ebPostfix.str(),ctd);
334 
335  if(createEdgeBlocks_) {
336  mesh.addEdgeBlock("eblock"+ebPostfix.str(),
338  edge_ctd);
339  }
340  if(createFaceBlocks_) {
341  mesh.addFaceBlock("eblock"+ebPostfix.str(),
343  face_ctd);
344  }
345  }
346  }
347  }
348 
349  // add sidesets
350  mesh.addSideset("left",side_ctd);
351  mesh.addSideset("right",side_ctd);
352  mesh.addSideset("top",side_ctd);
353  mesh.addSideset("bottom",side_ctd);
354  mesh.addSideset("front",side_ctd);
355  mesh.addSideset("back",side_ctd);
356 
358  for(int bx=1;bx<xBlocks_;bx++) {
359  std::stringstream ss;
360  ss << "vertical_" << bx-1;
361  mesh.addSideset(ss.str(),side_ctd);
362  }
363  for(int by=1;by<yBlocks_;by++) {
364  std::stringstream ss;
365  ss << "horizontal_" << by-1;
366  mesh.addSideset(ss.str(),side_ctd);
367  }
368  for(int bz=1;bz<zBlocks_;bz++) {
369  std::stringstream ss;
370  ss << "transverse_" << bz-1;
371  mesh.addSideset(ss.str(),side_ctd);
372  }
373  }
374 
375  // add nodesets
376  mesh.addNodeset("origin");
377 }
378 
379 void CubeHexMeshFactory::buildElements(stk::ParallelMachine parallelMach,STK_Interface & mesh) const
380 {
381  mesh.beginModification();
382  // build each block
383  for(int xBlock=0;xBlock<xBlocks_;xBlock++) {
384  for(int yBlock=0;yBlock<yBlocks_;yBlock++) {
385  for(int zBlock=0;zBlock<zBlocks_;zBlock++) {
386  buildBlock(parallelMach,xBlock,yBlock,zBlock,mesh);
387  }
388  }
389  }
390  mesh.endModification();
391 }
392 
393 void CubeHexMeshFactory::buildBlock(stk::ParallelMachine /* parallelMach */,int xBlock,int yBlock,int zBlock,STK_Interface & mesh) const
394 {
395  // grab this processors rank and machine size
396  std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> sizeAndStartX = determineXElemSizeAndStart(xBlock,xProcs_,machRank_);
397  std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> sizeAndStartY = determineYElemSizeAndStart(yBlock,yProcs_,machRank_);
398  std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> sizeAndStartZ = determineZElemSizeAndStart(zBlock,zProcs_,machRank_);
399 
400  panzer::GlobalOrdinal myXElems_start = sizeAndStartX.first;
401  panzer::GlobalOrdinal myXElems_end = myXElems_start+sizeAndStartX.second;
402  panzer::GlobalOrdinal myYElems_start = sizeAndStartY.first;
403  panzer::GlobalOrdinal myYElems_end = myYElems_start+sizeAndStartY.second;
404  panzer::GlobalOrdinal myZElems_start = sizeAndStartZ.first;
405  panzer::GlobalOrdinal myZElems_end = myZElems_start+sizeAndStartZ.second;
406 
407  panzer::GlobalOrdinal totalXElems = nXElems_*xBlocks_;
408  panzer::GlobalOrdinal totalYElems = nYElems_*yBlocks_;
409  panzer::GlobalOrdinal totalZElems = nZElems_*zBlocks_;
410 
411  double deltaX = (xf_-x0_)/double(totalXElems);
412  double deltaY = (yf_-y0_)/double(totalYElems);
413  double deltaZ = (zf_-z0_)/double(totalZElems);
414 
415  std::vector<double> coord(3,0.0);
416 
417  // build the nodes
418  for(panzer::GlobalOrdinal nx=myXElems_start;nx<myXElems_end+1;++nx) {
419  coord[0] = this->getMeshCoord(nx, deltaX, x0_);
420  for(panzer::GlobalOrdinal ny=myYElems_start;ny<myYElems_end+1;++ny) {
421  coord[1] = this->getMeshCoord(ny, deltaY, y0_);
422  for(panzer::GlobalOrdinal nz=myZElems_start;nz<myZElems_end+1;++nz) {
423  coord[2] = this->getMeshCoord(nz, deltaZ, z0_);
424 
425  mesh.addNode(nz*(totalYElems+1)*(totalXElems+1)+ny*(totalXElems+1)+nx+1,coord);
426  }
427  }
428  }
429 
430  std::stringstream blockName;
431  blockName << "eblock-" << xBlock << "_" << yBlock << "_" << zBlock;
432  stk::mesh::Part * block = mesh.getElementBlockPart(blockName.str());
433 
434  // build the elements
435  for(panzer::GlobalOrdinal nx=myXElems_start;nx<myXElems_end;++nx) {
436  for(panzer::GlobalOrdinal ny=myYElems_start;ny<myYElems_end;++ny) {
437  for(panzer::GlobalOrdinal nz=myZElems_start;nz<myZElems_end;++nz) {
438  stk::mesh::EntityId gid = totalXElems*totalYElems*nz+totalXElems*ny+nx+1;
439  std::vector<stk::mesh::EntityId> nodes(8);
440  nodes[0] = nx+1+ny*(totalXElems+1) +nz*(totalYElems+1)*(totalXElems+1);
441  nodes[1] = nodes[0]+1;
442  nodes[2] = nodes[1]+(totalXElems+1);
443  nodes[3] = nodes[2]-1;
444  nodes[4] = nodes[0]+(totalYElems+1)*(totalXElems+1);
445  nodes[5] = nodes[1]+(totalYElems+1)*(totalXElems+1);
446  nodes[6] = nodes[2]+(totalYElems+1)*(totalXElems+1);
447  nodes[7] = nodes[3]+(totalYElems+1)*(totalXElems+1);
448 
449  RCP<ElementDescriptor> ed = rcp(new ElementDescriptor(gid,nodes));
450  mesh.addElement(ed,block);
451  }
452  }
453  }
454 }
455 
456 std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> CubeHexMeshFactory::determineXElemSizeAndStart(int xBlock,unsigned int size,unsigned int /* rank */) const
457 {
458  std::size_t xProcLoc = procTuple_[0];
459  panzer::GlobalOrdinal minElements = nXElems_/size;
460  panzer::GlobalOrdinal extra = nXElems_ - minElements*size;
461 
462  TEUCHOS_ASSERT(minElements>0);
463 
464  // first "extra" elements get an extra column of elements
465  // this determines the starting X index and number of elements
466  panzer::GlobalOrdinal nume=0, start=0;
467  if(panzer::GlobalOrdinal(xProcLoc)<extra) {
468  nume = minElements+1;
469  start = xProcLoc*(minElements+1);
470  }
471  else {
472  nume = minElements;
473  start = extra*(minElements+1)+(xProcLoc-extra)*minElements;
474  }
475 
476  return std::make_pair(start+nXElems_*xBlock,nume);
477 }
478 
479 std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> CubeHexMeshFactory::determineYElemSizeAndStart(int yBlock,unsigned int size,unsigned int /* rank */) const
480 {
481  // int start = yBlock*nYElems_;
482  // return std::make_pair(start,nYElems_);
483 
484  std::size_t yProcLoc = procTuple_[1];
485  panzer::GlobalOrdinal minElements = nYElems_/size;
486  panzer::GlobalOrdinal extra = nYElems_ - minElements*size;
487 
488  TEUCHOS_ASSERT(minElements>0);
489 
490  // first "extra" elements get an extra column of elements
491  // this determines the starting X index and number of elements
492  panzer::GlobalOrdinal nume=0, start=0;
493  if(panzer::GlobalOrdinal(yProcLoc)<extra) {
494  nume = minElements+1;
495  start = yProcLoc*(minElements+1);
496  }
497  else {
498  nume = minElements;
499  start = extra*(minElements+1)+(yProcLoc-extra)*minElements;
500  }
501 
502  return std::make_pair(start+nYElems_*yBlock,nume);
503 }
504 
505 std::pair<panzer::GlobalOrdinal,panzer::GlobalOrdinal> CubeHexMeshFactory::determineZElemSizeAndStart(int zBlock,unsigned int size,unsigned int /* rank */) const
506 {
507  // int start = zBlock*nZElems_;
508  // return std::make_pair(start,nZElems_);
509  std::size_t zProcLoc = procTuple_[2];
510  panzer::GlobalOrdinal minElements = nZElems_/size;
511  panzer::GlobalOrdinal extra = nZElems_ - minElements*size;
512 
513  TEUCHOS_ASSERT(minElements>0);
514 
515  // first "extra" elements get an extra column of elements
516  // this determines the starting X index and number of elements
517  panzer::GlobalOrdinal nume=0, start=0;
518  if(zProcLoc<Teuchos::as<std::size_t>(extra)) {
519  nume = minElements+1;
520  start = zProcLoc*(minElements+1);
521  }
522  else {
523  nume = minElements;
524  start = extra*(minElements+1)+(zProcLoc-extra)*minElements;
525  }
526 
527  return std::make_pair(start+nZElems_*zBlock,nume);
528 }
529 
530 // this adds side entities only (does not inject them into side sets)
532 {
533  mesh.beginModification();
534 
535  std::size_t totalXElems = nXElems_*xBlocks_;
536  std::size_t totalYElems = nYElems_*yBlocks_;
537  std::size_t totalZElems = nZElems_*zBlocks_;
538 
539  std::vector<stk::mesh::Entity> localElmts;
540  mesh.getMyElements(localElmts);
541 
542  stk::mesh::EntityId offset[6];
543  offset[0] = 0;
544  offset[1] = offset[0] + totalXElems*totalZElems;
545  offset[2] = offset[1] + totalYElems*totalZElems;
546  offset[3] = offset[2] + totalXElems*totalZElems;
547  offset[4] = offset[3] + totalYElems*totalZElems;
548  offset[5] = offset[4] + totalXElems*totalYElems;
549 
550  // gid = totalXElems*totalYElems*nz+totalXElems*ny+nx+1
551 
552  // loop over elements adding sides to sidesets
553  std::vector<stk::mesh::Entity>::const_iterator itr;
554  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
555  stk::mesh::Entity element = (*itr);
556  stk::mesh::EntityId gid = mesh.elementGlobalId(element);
557 
558  std::size_t nx,ny,nz;
559  nz = (gid-1) / (totalXElems*totalYElems);
560  gid = (gid-1)-nz*(totalXElems*totalYElems);
561  ny = gid / totalXElems;
562  nx = gid-ny*totalXElems;
563 
564  std::vector<stk::mesh::Part*> parts;
565 
566  if(nz==0) {
567  // on the back
568  mesh.getBulkData()->declare_element_side(element, 4, parts);
569  }
570  if(nz+1==totalZElems) {
571  // on the front
572  mesh.getBulkData()->declare_element_side(element, 5, parts);
573  }
574 
575  if(ny==0) {
576  // on the bottom
577  mesh.getBulkData()->declare_element_side(element, 0, parts);
578  }
579  if(ny+1==totalYElems) {
580  // on the top
581  mesh.getBulkData()->declare_element_side(element, 2, parts);
582  }
583 
584  if(nx==0) {
585  // on the left
586  mesh.getBulkData()->declare_element_side(element, 3, parts);
587  }
588  if(nx+1==totalXElems) {
589  // on the right
590  mesh.getBulkData()->declare_element_side(element, 1, parts);
591  }
592  }
593 
594  mesh.endModification();
595 }
596 
597 // Pre-Condition: call beginModification() before entry
598 // Post-Condition: call endModification() after exit
600 {
601  const stk::mesh::EntityRank side_rank = mesh.getSideRank();
602 
603  std::size_t totalXElems = nXElems_*xBlocks_;
604  std::size_t totalYElems = nYElems_*yBlocks_;
605  std::size_t totalZElems = nZElems_*zBlocks_;
606 
607  // get all part vectors
608  stk::mesh::Part * left = mesh.getSideset("left");
609  stk::mesh::Part * right = mesh.getSideset("right");
610  stk::mesh::Part * top = mesh.getSideset("top");
611  stk::mesh::Part * bottom = mesh.getSideset("bottom");
612  stk::mesh::Part * front = mesh.getSideset("front");
613  stk::mesh::Part * back = mesh.getSideset("back");
614 
615  std::vector<stk::mesh::Part*> vertical;
616  std::vector<stk::mesh::Part*> horizontal;
617  std::vector<stk::mesh::Part*> transverse;
618 
620  for(int bx=1;bx<xBlocks_;bx++) {
621  std::stringstream ss;
622  ss << "vertical_" << bx-1;
623  vertical.push_back(mesh.getSideset(ss.str()));
624  }
625  for(int by=1;by<yBlocks_;by++) {
626  std::stringstream ss;
627  ss << "horizontal_" << by-1;
628  horizontal.push_back(mesh.getSideset(ss.str()));
629  }
630  for(int bz=1;bz<zBlocks_;bz++) {
631  std::stringstream ss;
632  ss << "transverse_" << bz-1;
633  transverse.push_back(mesh.getSideset(ss.str()));
634  }
635  }
636 
637  std::vector<stk::mesh::Entity> localElmts;
638  mesh.getMyElements(localElmts);
639 
640  // gid = totalXElems*totalYElems*nz+totalXElems*ny+nx+1
641 
642  // loop over elements adding sides to sidesets
643  std::vector<stk::mesh::Entity>::const_iterator itr;
644  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
645  stk::mesh::Entity element = (*itr);
646  stk::mesh::EntityId gid = mesh.elementGlobalId(element);
647 
648  std::size_t nx,ny,nz;
649  nz = (gid-1) / (totalXElems*totalYElems);
650  gid = (gid-1)-nz*(totalXElems*totalYElems);
651  ny = gid / totalXElems;
652  nx = gid-ny*totalXElems;
653 
654  if(nz % nZElems_==0) {
655  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 4);
656 
657  // on the back
658  if(mesh.entityOwnerRank(side)==machRank_) {
659  if(nz==0) {
660  mesh.addEntityToSideset(side,back);
661  } else {
663  int index = nz/nZElems_-1;
664  mesh.addEntityToSideset(side,transverse[index]);
665  }
666  }
667  }
668  }
669  if((nz+1) % nZElems_==0) {
670  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 5);
671 
672  // on the front
673  if(mesh.entityOwnerRank(side)==machRank_) {
674  if(nz+1==totalZElems) {
675  mesh.addEntityToSideset(side,front);
676  } else {
678  int index = (nz+1)/nZElems_-1;
679  mesh.addEntityToSideset(side,transverse[index]);
680  }
681  }
682  }
683  }
684 
685  if(ny % nYElems_==0) {
686  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 0);
687 
688  // on the bottom
689  if(mesh.entityOwnerRank(side)==machRank_) {
690  if(ny==0) {
691  mesh.addEntityToSideset(side,bottom);
692  } else {
694  int index = ny/nYElems_-1;
695  mesh.addEntityToSideset(side,horizontal[index]);
696  }
697  }
698  }
699  }
700  if((ny+1) % nYElems_==0) {
701  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 2);
702 
703  // on the top
704  if(mesh.entityOwnerRank(side)==machRank_) {
705  if(ny+1==totalYElems) {
706  mesh.addEntityToSideset(side,top);
707  } else {
709  int index = (ny+1)/nYElems_-1;
710  mesh.addEntityToSideset(side,horizontal[index]);
711  }
712  }
713  }
714  }
715 
716  if(nx % nXElems_==0) {
717  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 3);
718 
719  // on the left
720  if(mesh.entityOwnerRank(side)==machRank_) {
721  if(nx==0) {
722  mesh.addEntityToSideset(side,left);
723  } else {
725  int index = nx/nXElems_-1;
726  mesh.addEntityToSideset(side,vertical[index]);
727  }
728  }
729  }
730  }
731  if((nx+1) % nXElems_==0) {
732  stk::mesh::Entity side = mesh.findConnectivityById(element, side_rank, 1);
733 
734  // on the right
735  if(mesh.entityOwnerRank(side)==machRank_) {
736  if(nx+1==totalXElems) {
737  mesh.addEntityToSideset(side,right);
738  } else {
740  int index = (nx+1)/nXElems_-1;
741  mesh.addEntityToSideset(side,vertical[index]);
742  }
743  }
744  }
745  }
746  }
747 }
748 
749 // Pre-Condition: call beginModification() before entry
750 // Post-Condition: call endModification() after exit
752 {
753  // get all part vectors
754  stk::mesh::Part * origin = mesh.getNodeset("origin");
755 
757  if(machRank_==0)
758  {
759  // add zero node to origin node set
760  stk::mesh::Entity node = bulkData->get_entity(mesh.getNodeRank(),1);
761  mesh.addEntityToNodeset(node,origin);
762  }
763 }
764 
765 // Pre-Condition: call beginModification() before entry
766 // Post-Condition: call endModification() after exit
768 {
771 
772  stk::mesh::Part * edge_block = mesh.getEdgeBlock(edgeBlockName_);
773 
774  stk::mesh::Selector owned_block = metaData->locally_owned_part();
775 
776  std::vector<stk::mesh::Entity> edges;
777  bulkData->get_entities(mesh.getEdgeRank(), owned_block, edges);
778  mesh.addEntitiesToEdgeBlock(edges, edge_block);
779 }
780 
781 // Pre-Condition: call beginModification() before entry
782 // Post-Condition: call endModification() after exit
784 {
787 
788  stk::mesh::Part * face_block = mesh.getFaceBlock(faceBlockName_);
789 
790  stk::mesh::Selector owned_block = metaData->locally_owned_part();
791 
792  std::vector<stk::mesh::Entity> faces;
793  bulkData->get_entities(mesh.getFaceRank(), owned_block, faces);
794  mesh.addEntitiesToFaceBlock(faces, face_block);
795 }
796 
799 {
800  std::size_t i=0,j=0,k=0;
801 
802  k = procRank/(xProcs_*yProcs_); procRank = procRank % (xProcs_*yProcs_);
803  j = procRank/xProcs_; procRank = procRank % xProcs_;
804  i = procRank;
805 
806  return Teuchos::tuple(i,j,k);
807 }
808 
809 } // end panzer_stk
void addSides(STK_Interface &mesh) const
void addNodeset(const std::string &name)
unsigned entityOwnerRank(stk::mesh::Entity entity) const
void addEntityToNodeset(stk::mesh::Entity entity, stk::mesh::Part *nodeset)
void addEdgeBlocks(STK_Interface &mesh) const
static void parsePeriodicBCList(const Teuchos::RCP< Teuchos::ParameterList > &pl, std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > &periodicBC, bool &useBBoxSearch)
void addFaceBlocks(STK_Interface &mesh) const
void buildMetaData(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
basic_FancyOStream & setShowProcRank(const bool showProcRank)
T & get(const std::string &name, T def_value)
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block part
void buildElements(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void addEntityToSideset(stk::mesh::Entity entity, stk::mesh::Part *sideset)
void addSideSets(STK_Interface &mesh) const
stk::mesh::Part * getNodeset(const std::string &name) const
stk::mesh::Part * getFaceBlock(const std::string &name) const
get the block part
void addEdgeBlock(const std::string &elemBlockName, const std::string &edgeBlockName, const stk::topology &topology)
stk::mesh::Part * getSideset(const std::string &name) const
void addElement(const Teuchos::RCP< ElementDescriptor > &ed, stk::mesh::Part *block)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
stk::mesh::EntityRank getSideRank() const
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true, const bool buildRefinementSupport=false)
void rebalance(STK_Interface &mesh) const
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
void buildBlock(stk::ParallelMachine machRank, int xBlock, int yBlock, int zBlock, STK_Interface &mesh) const
void addNodeSets(STK_Interface &mesh) const
stk::mesh::EntityId elementGlobalId(std::size_t lid) const
void setMyParamList(const RCP< ParameterList > &paramList)
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
static const std::string edgeBlockString
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
From ParameterListAcceptor.
std::pair< panzer::GlobalOrdinal, panzer::GlobalOrdinal > determineZElemSizeAndStart(int zBlock, unsigned int size, unsigned int rank) const
stk::mesh::EntityRank getFaceRank() 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
Teuchos::Tuple< std::size_t, 3 > procTuple_
Teuchos::RCP< STK_Interface > buildMesh(stk::ParallelMachine parallelMach) const
Build the mesh object.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
From ParameterListAcceptor.
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
bool isInitialized() const
Has initialize been called on this mesh object?
stk::mesh::Entity findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const
void addSideset(const std::string &name, const CellTopologyData *ctData)
static const std::string faceBlockString
double getMeshCoord(const int nx, const double deltaX, const double x0) const
std::pair< panzer::GlobalOrdinal, panzer::GlobalOrdinal > determineXElemSizeAndStart(int xBlock, unsigned int size, unsigned int rank) const
stk::mesh::Part * getEdgeBlock(const std::string &name) const
get the block part
std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > periodicBCVec_
stk::mesh::EntityRank getEdgeRank() const
void addEntitiesToFaceBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *faceblock)
void addFaceBlock(const std::string &elemBlockName, const std::string &faceBlockName, const stk::topology &topology)
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
std::pair< panzer::GlobalOrdinal, panzer::GlobalOrdinal > determineYElemSizeAndStart(int yBlock, unsigned int size, unsigned int rank) const
void getMyElements(std::vector< stk::mesh::Entity > &elements) const
#define TEUCHOS_ASSERT(assertion_test)
stk::mesh::EntityRank getNodeRank() const
virtual Teuchos::RCP< STK_Interface > buildUncommitedMesh(stk::ParallelMachine parallelMach) const
void addElementBlock(const std::string &name, const CellTopologyData *ctData)
virtual void completeMeshConstruction(STK_Interface &mesh, stk::ParallelMachine parallelMach) const
Teuchos::Tuple< std::size_t, 3 > procRankToProcTuple(std::size_t procRank) const
what is the 3D tuple describe this processor distribution