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 <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 
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 
185 {
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 
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
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::Ordinal64,panzer::Ordinal64> sizeAndStartX = determineXElemSizeAndStart(xBlock,xProcs_,machRank_);
336  std::pair<panzer::Ordinal64,panzer::Ordinal64> sizeAndStartY = determineYElemSizeAndStart(yBlock,yProcs_,machRank_);
337  std::pair<panzer::Ordinal64,panzer::Ordinal64> sizeAndStartZ = determineZElemSizeAndStart(zBlock,zProcs_,machRank_);
338 
339  panzer::Ordinal64 myXElems_start = sizeAndStartX.first;
340  panzer::Ordinal64 myXElems_end = myXElems_start+sizeAndStartX.second;
341  panzer::Ordinal64 myYElems_start = sizeAndStartY.first;
342  panzer::Ordinal64 myYElems_end = myYElems_start+sizeAndStartY.second;
343  panzer::Ordinal64 myZElems_start = sizeAndStartZ.first;
344  panzer::Ordinal64 myZElems_end = myZElems_start+sizeAndStartZ.second;
345 
346  panzer::Ordinal64 totalXElems = nXElems_*xBlocks_;
347  panzer::Ordinal64 totalYElems = nYElems_*yBlocks_;
348  panzer::Ordinal64 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::Ordinal64 nx=myXElems_start;nx<myXElems_end+1;++nx) {
358  coord[0] = double(nx)*deltaX+x0_;
359  for(panzer::Ordinal64 ny=myYElems_start;ny<myYElems_end+1;++ny) {
360  coord[1] = double(ny)*deltaY+y0_;
361  for(panzer::Ordinal64 nz=myZElems_start;nz<myZElems_end+1;++nz) {
362  coord[2] = double(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::Ordinal64 nx=myXElems_start;nx<myXElems_end;++nx) {
375  for(panzer::Ordinal64 ny=myYElems_start;ny<myYElems_end;++ny) {
376  for(panzer::Ordinal64 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::Ordinal64,panzer::Ordinal64> CubeHexMeshFactory::determineXElemSizeAndStart(int xBlock,unsigned int size,unsigned int /* rank */) const
396 {
397  std::size_t xProcLoc = procTuple_[0];
398  panzer::Ordinal64 minElements = nXElems_/size;
399  panzer::Ordinal64 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::Ordinal64 nume=0, start=0;
406  if(panzer::Ordinal64(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::Ordinal64,panzer::Ordinal64> 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::Ordinal64 minElements = nYElems_/size;
425  panzer::Ordinal64 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::Ordinal64 nume=0, start=0;
432  if(panzer::Ordinal64(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::Ordinal64,panzer::Ordinal64> 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::Ordinal64 minElements = nZElems_/size;
450  panzer::Ordinal64 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::Ordinal64 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 
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 
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
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 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 count
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
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true)
stk::mesh::Part * getSideset(const std::string &name) const
std::pair< panzer::Ordinal64, panzer::Ordinal64 > determineXElemSizeAndStart(int xBlock, unsigned int size, unsigned int rank) 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
std::pair< panzer::Ordinal64, panzer::Ordinal64 > determineYElemSizeAndStart(int yBlock, unsigned int size, unsigned int rank) const
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)
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
From ParameterListAcceptor.
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_
std::pair< panzer::Ordinal64, panzer::Ordinal64 > determineZElemSizeAndStart(int zBlock, unsigned int size, unsigned int rank) const
Teuchos::RCP< STK_Interface > buildMesh(stk::ParallelMachine parallelMach) const
Build the mesh object.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
From ParameterListAcceptor.
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)
std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > periodicBCVec_
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
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)
static void parsePeriodicBCList(const Teuchos::RCP< Teuchos::ParameterList > &pl, std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > &periodicBC)
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