44 #include <Teuchos_TimeMonitor.hpp> 45 #include <PanzerAdaptersSTK_config.hpp> 46 #include <stk_mesh/base/FEMHelpers.hpp> 66 PANZER_FUNC_TIME_MONITOR(
"panzer::CubeHexMeshFactory::buildMesh()");
72 mesh->initialize(parallelMach);
82 PANZER_FUNC_TIME_MONITOR(
"panzer::CubeHexMeshFactory::buildUncomittedMesh()");
86 machRank_ = stk::parallel_machine_rank(parallelMach);
87 machSize_ = stk::parallel_machine_size(parallelMach);
99 const int maxFactor = 50;
102 int factors[maxFactor];
103 for (
int jj = 0; jj < maxFactor; jj++) factors[jj] = 0;
104 for (
int jj = 2; jj < maxFactor; jj++) {
107 int temp = ProcTemp/jj;
108 if (temp*jj == ProcTemp) {
118 for (
int jj = maxFactor-1; jj > 0; jj--) {
119 while (factors[jj] != 0) {
135 "Cannot build CubeHexMeshFactory, the product of \"X Procs\", \"Y Procs\", and \"Z Procs\"" 136 " must equal the number of processors.");
149 PANZER_FUNC_TIME_MONITOR(
"panzer::CubeHexMeshFactory::completeMeshConstruction()");
157 Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
158 out.setOutputToRootOnly(0);
159 out.setShowProcRank(
true);
165 out <<
"CubeHexMesh: Building sub cells" << std::endl;
170 out <<
"CubeHexMesh: NOT building sub cells" << std::endl;
204 setMyParamList(paramList);
206 x0_ = paramList->get<
double>(
"X0");
207 y0_ = paramList->get<
double>(
"Y0");
208 z0_ = paramList->get<
double>(
"Z0");
210 xf_ = paramList->get<
double>(
"Xf");
211 yf_ = paramList->get<
double>(
"Yf");
212 zf_ = paramList->get<
double>(
"Zf");
214 xBlocks_ = paramList->get<
int>(
"X Blocks");
215 yBlocks_ = paramList->get<
int>(
"Y Blocks");
216 zBlocks_ = paramList->get<
int>(
"Z Blocks");
218 xProcs_ = paramList->get<
int>(
"X Procs");
219 yProcs_ = paramList->get<
int>(
"Y Procs");
220 zProcs_ = paramList->get<
int>(
"Z Procs");
222 nXElems_ = paramList->get<
int>(
"X Elements");
223 nYElems_ = paramList->get<
int>(
"Y Elements");
224 nZElems_ = paramList->get<
int>(
"Z Elements");
233 Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
234 out.setOutputToRootOnly(0);
235 out.setShowProcRank(
true);
237 out <<
"CubeHexMesh: NOT creating edge blocks because building sub cells disabled" << std::endl;
241 Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
242 out.setOutputToRootOnly(0);
243 out.setShowProcRank(
true);
245 out <<
"CubeHexMesh: NOT creating face blocks because building sub cells disabled" << std::endl;
256 static RCP<Teuchos::ParameterList> defaultParams;
259 if(defaultParams == Teuchos::null) {
260 defaultParams = rcp(
new Teuchos::ParameterList);
262 defaultParams->set<
double>(
"X0",0.0);
263 defaultParams->set<
double>(
"Y0",0.0);
264 defaultParams->set<
double>(
"Z0",0.0);
266 defaultParams->set<
double>(
"Xf",1.0);
267 defaultParams->set<
double>(
"Yf",1.0);
268 defaultParams->set<
double>(
"Zf",1.0);
270 defaultParams->set<
int>(
"X Blocks",1);
271 defaultParams->set<
int>(
"Y Blocks",1);
272 defaultParams->set<
int>(
"Z Blocks",1);
274 defaultParams->set<
int>(
"X Procs",-1);
275 defaultParams->set<
int>(
"Y Procs",1);
276 defaultParams->set<
int>(
"Z Procs",1);
278 defaultParams->set<
int>(
"X Elements",5);
279 defaultParams->set<
int>(
"Y Elements",5);
280 defaultParams->set<
int>(
"Z Elements",5);
282 defaultParams->set<
bool>(
"Build Interface Sidesets",
false);
284 defaultParams->set<
bool>(
"Build Subcells",
true);
287 defaultParams->set<
bool>(
"Create Edge Blocks",
false,
"Create edge blocks in the mesh");
288 defaultParams->set<
bool>(
"Create Face Blocks",
false,
"Create face blocks in the mesh");
290 Teuchos::ParameterList & bcs = defaultParams->sublist(
"Periodic BCs");
291 bcs.set<
int>(
"Count",0);
294 return defaultParams;
300 RCP<Teuchos::ParameterList> validParams = rcp(
new Teuchos::ParameterList(*
getValidParameters()));
308 typedef shards::Hexahedron<8> HexTopo;
309 const CellTopologyData * ctd = shards::getCellTopologyData<HexTopo>();
310 const CellTopologyData * side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(2,0);
318 std::stringstream ebPostfix;
319 ebPostfix <<
"-" << bx <<
"_" << by <<
"_" << bz;
337 std::stringstream ss;
338 ss <<
"vertical_" << bx-1;
342 std::stringstream ss;
343 ss <<
"horizontal_" << by-1;
347 std::stringstream ss;
348 ss <<
"transverse_" << bz-1;
357 const CellTopologyData * edge_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(1,0);
361 const CellTopologyData * face_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(2,0);
370 for(
int xBlock=0;xBlock<
xBlocks_;xBlock++) {
371 for(
int yBlock=0;yBlock<
yBlocks_;yBlock++) {
372 for(
int zBlock=0;zBlock<
zBlocks_;zBlock++) {
373 buildBlock(parallelMach,xBlock,yBlock,zBlock,mesh);
387 panzer::GlobalOrdinal myXElems_start = sizeAndStartX.first;
388 panzer::GlobalOrdinal myXElems_end = myXElems_start+sizeAndStartX.second;
389 panzer::GlobalOrdinal myYElems_start = sizeAndStartY.first;
390 panzer::GlobalOrdinal myYElems_end = myYElems_start+sizeAndStartY.second;
391 panzer::GlobalOrdinal myZElems_start = sizeAndStartZ.first;
392 panzer::GlobalOrdinal myZElems_end = myZElems_start+sizeAndStartZ.second;
398 double deltaX = (
xf_-
x0_)/
double(totalXElems);
399 double deltaY = (
yf_-
y0_)/
double(totalYElems);
400 double deltaZ = (
zf_-
z0_)/
double(totalZElems);
402 std::vector<double> coord(3,0.0);
405 for(panzer::GlobalOrdinal nx=myXElems_start;nx<myXElems_end+1;++nx) {
407 for(panzer::GlobalOrdinal ny=myYElems_start;ny<myYElems_end+1;++ny) {
409 for(panzer::GlobalOrdinal nz=myZElems_start;nz<myZElems_end+1;++nz) {
412 mesh.
addNode(nz*(totalYElems+1)*(totalXElems+1)+ny*(totalXElems+1)+nx+1,coord);
417 std::stringstream blockName;
418 blockName <<
"eblock-" << xBlock <<
"_" << yBlock <<
"_" << zBlock;
422 for(panzer::GlobalOrdinal nx=myXElems_start;nx<myXElems_end;++nx) {
423 for(panzer::GlobalOrdinal ny=myYElems_start;ny<myYElems_end;++ny) {
424 for(panzer::GlobalOrdinal nz=myZElems_start;nz<myZElems_end;++nz) {
425 stk::mesh::EntityId gid = totalXElems*totalYElems*nz+totalXElems*ny+nx+1;
426 std::vector<stk::mesh::EntityId> nodes(8);
427 nodes[0] = nx+1+ny*(totalXElems+1) +nz*(totalYElems+1)*(totalXElems+1);
428 nodes[1] = nodes[0]+1;
429 nodes[2] = nodes[1]+(totalXElems+1);
430 nodes[3] = nodes[2]-1;
431 nodes[4] = nodes[0]+(totalYElems+1)*(totalXElems+1);
432 nodes[5] = nodes[1]+(totalYElems+1)*(totalXElems+1);
433 nodes[6] = nodes[2]+(totalYElems+1)*(totalXElems+1);
434 nodes[7] = nodes[3]+(totalYElems+1)*(totalXElems+1);
446 panzer::GlobalOrdinal minElements =
nXElems_/size;
447 panzer::GlobalOrdinal extra =
nXElems_ - minElements*size;
449 TEUCHOS_ASSERT(minElements>0);
453 panzer::GlobalOrdinal nume=0, start=0;
454 if(panzer::GlobalOrdinal(xProcLoc)<extra) {
455 nume = minElements+1;
456 start = xProcLoc*(minElements+1);
460 start = extra*(minElements+1)+(xProcLoc-extra)*minElements;
463 return std::make_pair(start+
nXElems_*xBlock,nume);
472 panzer::GlobalOrdinal minElements =
nYElems_/size;
473 panzer::GlobalOrdinal extra =
nYElems_ - minElements*size;
475 TEUCHOS_ASSERT(minElements>0);
479 panzer::GlobalOrdinal nume=0, start=0;
480 if(panzer::GlobalOrdinal(yProcLoc)<extra) {
481 nume = minElements+1;
482 start = yProcLoc*(minElements+1);
486 start = extra*(minElements+1)+(yProcLoc-extra)*minElements;
489 return std::make_pair(start+
nYElems_*yBlock,nume);
497 panzer::GlobalOrdinal minElements =
nZElems_/size;
498 panzer::GlobalOrdinal extra =
nZElems_ - minElements*size;
500 TEUCHOS_ASSERT(minElements>0);
504 panzer::GlobalOrdinal nume=0, start=0;
505 if(zProcLoc<Teuchos::as<std::size_t>(extra)) {
506 nume = minElements+1;
507 start = zProcLoc*(minElements+1);
511 start = extra*(minElements+1)+(zProcLoc-extra)*minElements;
514 return std::make_pair(start+
nZElems_*zBlock,nume);
526 std::vector<stk::mesh::Entity> localElmts;
529 stk::mesh::EntityId offset[6];
531 offset[1] = offset[0] + totalXElems*totalZElems;
532 offset[2] = offset[1] + totalYElems*totalZElems;
533 offset[3] = offset[2] + totalXElems*totalZElems;
534 offset[4] = offset[3] + totalYElems*totalZElems;
535 offset[5] = offset[4] + totalXElems*totalYElems;
540 std::vector<stk::mesh::Entity>::const_iterator itr;
541 for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
542 stk::mesh::Entity element = (*itr);
545 std::size_t nx,ny,nz;
546 nz = (gid-1) / (totalXElems*totalYElems);
547 gid = (gid-1)-nz*(totalXElems*totalYElems);
548 ny = gid / totalXElems;
549 nx = gid-ny*totalXElems;
551 std::vector<stk::mesh::Part*> parts;
555 mesh.
getBulkData()->declare_element_side(element, 4, parts);
557 if(nz+1==totalZElems) {
559 mesh.
getBulkData()->declare_element_side(element, 5, parts);
564 mesh.
getBulkData()->declare_element_side(element, 0, parts);
566 if(ny+1==totalYElems) {
568 mesh.
getBulkData()->declare_element_side(element, 2, parts);
573 mesh.
getBulkData()->declare_element_side(element, 3, parts);
575 if(nx+1==totalXElems) {
577 mesh.
getBulkData()->declare_element_side(element, 1, parts);
588 const stk::mesh::EntityRank side_rank = mesh.
getSideRank();
595 stk::mesh::Part * left = mesh.
getSideset(
"left");
596 stk::mesh::Part * right = mesh.
getSideset(
"right");
597 stk::mesh::Part * top = mesh.
getSideset(
"top");
598 stk::mesh::Part * bottom = mesh.
getSideset(
"bottom");
599 stk::mesh::Part * front = mesh.
getSideset(
"front");
600 stk::mesh::Part * back = mesh.
getSideset(
"back");
602 std::vector<stk::mesh::Part*> vertical;
603 std::vector<stk::mesh::Part*> horizontal;
604 std::vector<stk::mesh::Part*> transverse;
608 std::stringstream ss;
609 ss <<
"vertical_" << bx-1;
610 vertical.push_back(mesh.
getSideset(ss.str()));
613 std::stringstream ss;
614 ss <<
"horizontal_" << by-1;
615 horizontal.push_back(mesh.
getSideset(ss.str()));
618 std::stringstream ss;
619 ss <<
"transverse_" << bz-1;
620 transverse.push_back(mesh.
getSideset(ss.str()));
624 std::vector<stk::mesh::Entity> localElmts;
630 std::vector<stk::mesh::Entity>::const_iterator itr;
631 for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
632 stk::mesh::Entity element = (*itr);
635 std::size_t nx,ny,nz;
636 nz = (gid-1) / (totalXElems*totalYElems);
637 gid = (gid-1)-nz*(totalXElems*totalYElems);
638 ny = gid / totalXElems;
639 nx = gid-ny*totalXElems;
661 if(nz+1==totalZElems) {
692 if(ny+1==totalYElems) {
723 if(nx+1==totalXElems) {
741 stk::mesh::Part * origin = mesh.
getNodeset(
"origin");
743 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.
getBulkData();
747 stk::mesh::Entity node = bulkData->get_entity(mesh.
getNodeRank(),1);
758 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.
getBulkData();
759 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.
getMetaData();
761 std::vector<stk::mesh::Entity> edges;
762 bulkData->get_entities(mesh.
getEdgeRank(),metaData->locally_owned_part(),edges);
772 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.
getBulkData();
773 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.
getMetaData();
775 std::vector<stk::mesh::Entity> faces;
776 bulkData->get_entities(mesh.
getFaceRank(),metaData->locally_owned_part(),faces);
783 std::size_t i=0,j=0,k=0;
789 return Teuchos::tuple(i,j,k);
std::pair< panzer::GlobalOrdinal, panzer::GlobalOrdinal > determineXElemSizeAndStart(int xBlock, unsigned int size, unsigned int rank) const
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block part
void initializeWithDefaults()
virtual ~CubeHexMeshFactory()
Destructor.
void addNodeset(const std::string &name)
void addSides(STK_Interface &mesh) const
void addEntityToNodeset(stk::mesh::Entity entity, stk::mesh::Part *nodeset)
stk::mesh::EntityRank getFaceRank() const
virtual Teuchos::RCP< STK_Interface > buildUncommitedMesh(stk::ParallelMachine parallelMach) const
panzer::GlobalOrdinal nXElems_
void getMyElements(std::vector< stk::mesh::Entity > &elements) const
void addEntityToSideset(stk::mesh::Entity entity, stk::mesh::Part *sideset)
stk::mesh::EntityRank getNodeRank() const
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
void addSideSets(STK_Interface &mesh) const
virtual void completeMeshConstruction(STK_Interface &mesh, stk::ParallelMachine parallelMach) const
stk::mesh::EntityRank getSideRank() const
stk::mesh::Entity findConnectivityById(stk::mesh::Entity src, stk::mesh::EntityRank tgt_rank, unsigned rel_id) const
stk::mesh::Part * getEdgeBlock(const std::string &name) const
get the block part
std::pair< panzer::GlobalOrdinal, panzer::GlobalOrdinal > determineZElemSizeAndStart(int zBlock, unsigned int size, unsigned int rank) const
stk::mesh::EntityRank getEdgeRank() const
stk::mesh::Part * getSideset(const std::string &name) const
void addFaceBlocks(STK_Interface &mesh) const
void addElement(const Teuchos::RCP< ElementDescriptor > &ed, stk::mesh::Part *block)
stk::mesh::EntityId elementGlobalId(std::size_t lid) const
void buildMetaData(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true, const bool buildRefinementSupport=false)
Teuchos::RCP< STK_Interface > buildMesh(stk::ParallelMachine parallelMach) const
Build the mesh object.
double getMeshCoord(const int nx, const double deltaX, const double x0) const
stk::mesh::Part * getFaceBlock(const std::string &name) const
get the block part
bool isInitialized() const
Has initialize been called on this mesh object?
static const std::string edgeBlockString
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > ¶mList)
From ParameterListAcceptor.
Teuchos::Tuple< std::size_t, 3 > procRankToProcTuple(std::size_t procRank) const
what is the 3D tuple describe this processor distribution
void addFaceBlock(const std::string &name, const CellTopologyData *ctData)
void addNode(stk::mesh::EntityId gid, const std::vector< double > &coord)
void addEntitiesToEdgeBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *edgeblock)
void buildSubcells()
force the mesh to build subcells: edges and faces
Teuchos::Tuple< std::size_t, 3 > procTuple_
CubeHexMeshFactory()
Constructor.
void buildElements(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
void addNodeSets(STK_Interface &mesh) const
void addSideset(const std::string &name, const CellTopologyData *ctData)
static const std::string faceBlockString
void addEdgeBlock(const std::string &name, const CellTopologyData *ctData)
std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > periodicBCVec_
void addEntitiesToFaceBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *faceblock)
bool buildInterfaceSidesets_
void buildLocalElementIDs()
unsigned entityOwnerRank(stk::mesh::Entity entity) const
panzer::GlobalOrdinal nYElems_
void rebalance(STK_Interface &mesh) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
From ParameterListAcceptor.
std::pair< panzer::GlobalOrdinal, panzer::GlobalOrdinal > determineYElemSizeAndStart(int yBlock, unsigned int size, unsigned int rank) const
void addEdgeBlocks(STK_Interface &mesh) const
panzer::GlobalOrdinal nZElems_
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)