44#include <Teuchos_TimeMonitor.hpp>
45#include <PanzerAdaptersSTK_config.hpp>
65 PANZER_FUNC_TIME_MONITOR(
"panzer::CubeTetMeshFactory::buildMesh()");
71 mesh->initialize(parallelMach);
81 PANZER_FUNC_TIME_MONITOR(
"panzer::CubeTetMeshFactory::buildUncomittedMesh()");
85 machRank_ = stk::parallel_machine_rank(parallelMach);
86 machSize_ = stk::parallel_machine_size(parallelMach);
98 const int maxFactor = 50;
101 int factors[maxFactor];
102 for (
int jj = 0; jj < maxFactor; jj++) factors[jj] = 0;
103 for (
int jj = 2; jj < maxFactor; jj++) {
106 int temp = ProcTemp/jj;
107 if (temp*jj == ProcTemp) {
117 for (
int jj = maxFactor-1; jj > 0; jj--) {
118 while (factors[jj] != 0) {
134 "Cannot build CubeTetMeshFactory, the product of \"X Procs\", \"Y Procs\", and \"Z Procs\""
135 " must equal the number of processors.");
149 PANZER_FUNC_TIME_MONITOR(
"panzer::CubeTetMeshFactory::completeMeshConstruction()");
189 setMyParamList(paramList);
191 x0_ = paramList->get<
double>(
"X0");
192 y0_ = paramList->get<
double>(
"Y0");
193 z0_ = paramList->get<
double>(
"Z0");
195 xf_ = paramList->get<
double>(
"Xf");
196 yf_ = paramList->get<
double>(
"Yf");
197 zf_ = paramList->get<
double>(
"Zf");
199 xBlocks_ = paramList->get<
int>(
"X Blocks");
200 yBlocks_ = paramList->get<
int>(
"Y Blocks");
201 zBlocks_ = paramList->get<
int>(
"Z Blocks");
203 xProcs_ = paramList->get<
int>(
"X Procs");
204 yProcs_ = paramList->get<
int>(
"Y Procs");
205 zProcs_ = paramList->get<
int>(
"Z Procs");
207 nXElems_ = paramList->get<
int>(
"X Elements");
208 nYElems_ = paramList->get<
int>(
"Y Elements");
209 nZElems_ = paramList->get<
int>(
"Z Elements");
221 static RCP<Teuchos::ParameterList> defaultParams;
224 if(defaultParams == Teuchos::null) {
225 defaultParams = rcp(
new Teuchos::ParameterList);
227 defaultParams->set<
double>(
"X0",0.0);
228 defaultParams->set<
double>(
"Y0",0.0);
229 defaultParams->set<
double>(
"Z0",0.0);
231 defaultParams->set<
double>(
"Xf",1.0);
232 defaultParams->set<
double>(
"Yf",1.0);
233 defaultParams->set<
double>(
"Zf",1.0);
235 defaultParams->set<
int>(
"X Blocks",1);
236 defaultParams->set<
int>(
"Y Blocks",1);
237 defaultParams->set<
int>(
"Z Blocks",1);
239 defaultParams->set<
int>(
"X Procs",-1);
240 defaultParams->set<
int>(
"Y Procs",1);
241 defaultParams->set<
int>(
"Z Procs",1);
243 defaultParams->set<
int>(
"X Elements",5);
244 defaultParams->set<
int>(
"Y Elements",5);
245 defaultParams->set<
int>(
"Z Elements",5);
248 defaultParams->set<
bool>(
"Create Edge Blocks",
false,
"Create edge blocks in the mesh");
249 defaultParams->set<
bool>(
"Create Face Blocks",
false,
"Create face blocks in the mesh");
251 Teuchos::ParameterList & bcs = defaultParams->sublist(
"Periodic BCs");
252 bcs.set<
int>(
"Count",0);
255 return defaultParams;
261 RCP<Teuchos::ParameterList> validParams = rcp(
new Teuchos::ParameterList(*
getValidParameters()));
277 typedef shards::Tetrahedron<4> TetTopo;
278 const CellTopologyData * ctd = shards::getCellTopologyData<TetTopo>();
279 const CellTopologyData * side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(2,0);
280 const CellTopologyData * edge_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(1,0);
281 const CellTopologyData * face_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(2,0);
289 std::stringstream ebPostfix;
290 ebPostfix <<
"-" << bx <<
"_" << by <<
"_" << bz;
323 for(
int xBlock=0;xBlock<
xBlocks_;xBlock++) {
324 for(
int yBlock=0;yBlock<
yBlocks_;yBlock++) {
325 for(
int zBlock=0;zBlock<
zBlocks_;zBlock++) {
326 buildBlock(parallelMach,xBlock,yBlock,zBlock,mesh);
340 int myXElems_start = sizeAndStartX.first;
341 int myXElems_end = myXElems_start+sizeAndStartX.second;
342 int myYElems_start = sizeAndStartY.first;
343 int myYElems_end = myYElems_start+sizeAndStartY.second;
344 int myZElems_start = sizeAndStartZ.first;
345 int myZElems_end = myZElems_start+sizeAndStartZ.second;
351 double deltaX = (
xf_-
x0_)/
double(totalXElems);
352 double deltaY = (
yf_-
y0_)/
double(totalYElems);
353 double deltaZ = (
zf_-
z0_)/
double(totalZElems);
355 std::vector<double> coord(3,0.0);
358 for(
int nx=myXElems_start;nx<myXElems_end+1;++nx) {
360 for(
int ny=myYElems_start;ny<myYElems_end+1;++ny) {
362 for(
int nz=myZElems_start;nz<myZElems_end+1;++nz) {
365 mesh.
addNode(nz*(totalYElems+1)*(totalXElems+1)+ny*(totalXElems+1)+nx+1,coord);
370 std::stringstream blockName;
371 blockName <<
"eblock-" << xBlock <<
"_" << yBlock <<
"_" << zBlock;
375 for(
int nx=myXElems_start;nx<myXElems_end;++nx) {
376 for(
int ny=myYElems_start;ny<myYElems_end;++ny) {
377 for(
int nz=myZElems_start;nz<myZElems_end;++nz) {
379 std::vector<stk::mesh::EntityId> nodes(8);
380 nodes[0] = nx+1+ny*(totalXElems+1) +nz*(totalYElems+1)*(totalXElems+1);
381 nodes[1] = nodes[0]+1;
382 nodes[2] = nodes[1]+(totalXElems+1);
383 nodes[3] = nodes[2]-1;
384 nodes[4] = nodes[0]+(totalYElems+1)*(totalXElems+1);
385 nodes[5] = nodes[1]+(totalYElems+1)*(totalXElems+1);
386 nodes[6] = nodes[2]+(totalYElems+1)*(totalXElems+1);
387 nodes[7] = nodes[3]+(totalYElems+1)*(totalXElems+1);
389 buildTetsOnHex(Teuchos::tuple(totalXElems,totalYElems,totalZElems),
390 Teuchos::tuple(nx,ny,nz),
398 const Teuchos::Tuple<int,3> & element,
399 stk::mesh::Part * block,
400 const std::vector<stk::mesh::EntityId> & h_nodes,
403 Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
404 out.setShowProcRank(
true);
405 out.setOutputToRootOnly(-1);
407 int totalXElems = meshDesc[0];
int totalYElems = meshDesc[1];
int totalZElems = meshDesc[2];
408 int nx = element[0];
int ny = element[1];
int nz = element[2];
410 stk::mesh::EntityId hex_id = totalXElems*totalYElems*nz+totalXElems*ny+nx+1;
411 stk::mesh::EntityId gid_0 = 12*(hex_id-1)+1;
412 std::vector<stk::mesh::EntityId> nodes(4);
415 stk::mesh::EntityId centroid = 0;
417 stk::mesh::EntityId largestNode = (totalXElems+1)*(totalYElems+1)*(totalZElems+1);
418 centroid = hex_id+largestNode;
421 std::vector<double> coord(3,0.0);
422 for(std::size_t i=0;i<h_nodes.size();i++) {
424 coord[0] += node_coord[0];
425 coord[1] += node_coord[1];
426 coord[2] += node_coord[2];
436 int idSet[][3] = { { 0, 1, 2},
449 for(
int i=0;i<12;i++) {
450 nodes[0] = h_nodes[idSet[i][0]];
451 nodes[1] = h_nodes[idSet[i][1]];
452 nodes[2] = h_nodes[idSet[i][2]];
463 unsigned int minElements =
nXElems_/size;
464 unsigned int extra =
nXElems_ - minElements*size;
466 TEUCHOS_ASSERT(minElements>0);
472 nume = minElements+1;
473 start = xProcLoc*(minElements+1);
477 start = extra*(minElements+1)+(xProcLoc-extra)*minElements;
480 return std::make_pair(start+
nXElems_*xBlock,nume);
489 unsigned int minElements =
nYElems_/size;
490 unsigned int extra =
nYElems_ - minElements*size;
492 TEUCHOS_ASSERT(minElements>0);
498 nume = minElements+1;
499 start = yProcLoc*(minElements+1);
503 start = extra*(minElements+1)+(yProcLoc-extra)*minElements;
506 return std::make_pair(start+
nYElems_*yBlock,nume);
514 unsigned int minElements =
nZElems_/size;
515 unsigned int extra =
nZElems_ - minElements*size;
517 TEUCHOS_ASSERT(minElements>0);
523 nume = minElements+1;
524 start = zProcLoc*(minElements+1);
528 start = extra*(minElements+1)+(zProcLoc-extra)*minElements;
531 return std::make_pair(start+
nZElems_*zBlock,nume);
537 const stk::mesh::EntityRank side_rank = mesh.
getSideRank();
544 stk::mesh::Part * left = mesh.
getSideset(
"left");
545 stk::mesh::Part * right = mesh.
getSideset(
"right");
546 stk::mesh::Part * top = mesh.
getSideset(
"top");
547 stk::mesh::Part * bottom = mesh.
getSideset(
"bottom");
548 stk::mesh::Part * front = mesh.
getSideset(
"front");
549 stk::mesh::Part * back = mesh.
getSideset(
"back");
551 std::vector<stk::mesh::Entity> localElmts;
557 std::vector<stk::mesh::Entity>::const_iterator itr;
558 for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
559 stk::mesh::Entity element = (*itr);
563 stk::mesh::EntityId h_gid = (gid-1)/12+1;
564 stk::mesh::EntityId t_offset = gid - (12*(h_gid-1)+1);
566 std::size_t nx,ny,nz;
567 nz = (h_gid-1) / (totalXElems*totalYElems);
568 h_gid = (h_gid-1)-nz*(totalXElems*totalYElems);
569 ny = h_gid / totalXElems;
570 nx = h_gid-ny*totalXElems;
572 if(nz==0 && (t_offset==0 || t_offset==1)) {
579 if(nz+1==totalZElems && (t_offset==10 || t_offset==11)) {
587 if(ny==0 && (t_offset==2 || t_offset==3)) {
594 if(ny+1==totalYElems && (t_offset==8 || t_offset==9)) {
602 if(nx==0 && (t_offset==4 || t_offset==5)) {
609 if(nx+1==totalXElems && (t_offset==6 || t_offset==7)) {
626 stk::mesh::Part * origin = mesh.
getNodeset(
"origin");
628 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.
getBulkData();
632 stk::mesh::Entity node = bulkData->get_entity(mesh.
getNodeRank(),1);
643 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.
getBulkData();
644 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.
getMetaData();
648 stk::mesh::Selector owned_block = metaData->locally_owned_part();
650 std::vector<stk::mesh::Entity> edges;
651 bulkData->get_entities(mesh.
getEdgeRank(), owned_block, edges);
659 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.
getBulkData();
660 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.
getMetaData();
664 stk::mesh::Selector owned_block = metaData->locally_owned_part();
666 std::vector<stk::mesh::Entity> faces;
667 bulkData->get_entities(mesh.
getFaceRank(), owned_block, faces);
674 std::size_t i=0,j=0,k=0;
680 return Teuchos::tuple(i,j,k);
virtual ~CubeTetMeshFactory()
Destructor.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
From ParameterListAcceptor.
void addFaceBlocks(STK_Interface &mesh) const
void addEdgeBlocks(STK_Interface &mesh) const
virtual Teuchos::RCP< STK_Interface > buildUncommitedMesh(stk::ParallelMachine parallelMach) 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.
Teuchos::Tuple< std::size_t, 3 > procTuple_
void buildMetaData(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
void addSideSets(STK_Interface &mesh) const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > ¶mList)
From ParameterListAcceptor.
std::string faceBlockName_
void buildElements(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
std::string edgeBlockName_
std::pair< int, int > determineZElemSizeAndStart(int zBlock, unsigned int size, unsigned int rank) const
std::pair< int, int > determineXElemSizeAndStart(int xBlock, unsigned int size, unsigned int rank) const
virtual void completeMeshConstruction(STK_Interface &mesh, stk::ParallelMachine parallelMach) const
void addNodeSets(STK_Interface &mesh) const
void initializeWithDefaults()
void buildTetsOnHex(const Teuchos::Tuple< int, 3 > &meshDesc, const Teuchos::Tuple< int, 3 > &element, stk::mesh::Part *block, const std::vector< stk::mesh::EntityId > &h_nodes, STK_Interface &mesh) const
CubeTetMeshFactory()
Constructor.
std::pair< int, int > determineYElemSizeAndStart(int yBlock, unsigned int size, unsigned int rank) const
Teuchos::Tuple< std::size_t, 3 > procRankToProcTuple(std::size_t procRank) const
what is the 3D tuple describe this processor distribution
void buildLocalElementIDs()
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
const double * getNodeCoordinates(stk::mesh::EntityId nodeId) const
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_