Panzer Version of the Day
Loading...
Searching...
No Matches
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
48using Teuchos::RCP;
49using Teuchos::rcp;
50
51namespace panzer_stk {
52
57
62
64Teuchos::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
80Teuchos::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
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 mesh->setBoundingBoxSearchFlag(useBBoxSearch_);
144
145 return mesh;
146}
147
148void 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
176 mesh.buildLocalEdgeIDs();
177 }
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);
188 addEdgeBlocks(mesh);
189 }
191 addFaceBlocks(mesh);
192 }
193
194 mesh.endModification();
195
196 // calls Stk_MeshFactory::rebalance
197 this->rebalance(mesh);
198}
199
201void CubeHexMeshFactory::setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & paramList)
202{
203 paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
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");
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 }
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
255Teuchos::RCP<const Teuchos::ParameterList> CubeHexMeshFactory::getValidParameters() const
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
301 RCP<Teuchos::ParameterList> validParams = rcp(new Teuchos::ParameterList(*getValidParameters()));
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
315void 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
336 mesh.addEdgeBlock("eblock"+ebPostfix.str(),
338 edge_ctd);
339 }
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
379void 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
393void 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
456std::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
479std::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
505std::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
756 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
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{
769 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
770 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
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{
785 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
786 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
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
798Teuchos::Tuple<std::size_t,3> CubeHexMeshFactory::procRankToProcTuple(std::size_t procRank) const
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 buildElements(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
virtual Teuchos::RCP< STK_Interface > buildUncommitedMesh(stk::ParallelMachine parallelMach) const
std::pair< panzer::GlobalOrdinal, panzer::GlobalOrdinal > determineZElemSizeAndStart(int zBlock, unsigned int size, unsigned int rank) const
void buildMetaData(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
Teuchos::Tuple< std::size_t, 3 > procRankToProcTuple(std::size_t procRank) const
what is the 3D tuple describe this processor distribution
std::pair< panzer::GlobalOrdinal, panzer::GlobalOrdinal > determineYElemSizeAndStart(int yBlock, unsigned int size, unsigned int rank) const
std::pair< panzer::GlobalOrdinal, panzer::GlobalOrdinal > determineXElemSizeAndStart(int xBlock, unsigned int size, unsigned int rank) const
void addNodeSets(STK_Interface &mesh) const
Teuchos::Tuple< std::size_t, 3 > procTuple_
void addEdgeBlocks(STK_Interface &mesh) const
void buildBlock(stk::ParallelMachine machRank, int xBlock, int yBlock, int zBlock, STK_Interface &mesh) const
Teuchos::RCP< STK_Interface > buildMesh(stk::ParallelMachine parallelMach) const
Build the mesh object.
void addFaceBlocks(STK_Interface &mesh) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
From ParameterListAcceptor.
void addSideSets(STK_Interface &mesh) const
virtual void completeMeshConstruction(STK_Interface &mesh, stk::ParallelMachine parallelMach) const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
From ParameterListAcceptor.
stk::mesh::Entity findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true, const bool buildRefinementSupport=false)
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block part
static const std::string edgeBlockString
stk::mesh::EntityId elementGlobalId(std::size_t lid) const
stk::mesh::Part * getEdgeBlock(const std::string &name) const
get the block part
bool isInitialized() const
Has initialize been called on this mesh object?
void addEntitiesToFaceBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *faceblock)
stk::mesh::EntityRank getNodeRank() const
void addEdgeBlock(const std::string &elemBlockName, const std::string &edgeBlockName, const stk::topology &topology)
void buildSubcells()
force the mesh to build subcells: edges and faces
void addEntityToNodeset(stk::mesh::Entity entity, stk::mesh::Part *nodeset)
void addElement(const Teuchos::RCP< ElementDescriptor > &ed, stk::mesh::Part *block)
void addNode(stk::mesh::EntityId gid, const std::vector< double > &coord)
void addNodeset(const std::string &name)
void addSideset(const std::string &name, const CellTopologyData *ctData)
stk::mesh::EntityRank getFaceRank() const
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
stk::mesh::Part * getNodeset(const std::string &name) const
unsigned entityOwnerRank(stk::mesh::Entity entity) const
void addEntitiesToEdgeBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *edgeblock)
void addEntityToSideset(stk::mesh::Entity entity, stk::mesh::Part *sideset)
stk::mesh::EntityRank getSideRank() const
void getMyElements(std::vector< stk::mesh::Entity > &elements) const
void addElementBlock(const std::string &name, const CellTopologyData *ctData)
void addFaceBlock(const std::string &elemBlockName, const std::string &faceBlockName, const stk::topology &topology)
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
static const std::string faceBlockString
stk::mesh::Part * getFaceBlock(const std::string &name) const
get the block part
stk::mesh::EntityRank getEdgeRank() const
stk::mesh::Part * getSideset(const std::string &name) const
void rebalance(STK_Interface &mesh) const
static void parsePeriodicBCList(const Teuchos::RCP< Teuchos::ParameterList > &pl, std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > &periodicBC, bool &useBBoxSearch)
double getMeshCoord(const int nx, const double deltaX, const double x0) const
std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > periodicBCVec_