Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_STK_MultiBlockMeshFactory.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 
45 using Teuchos::RCP;
46 using Teuchos::rcp;
47 
48 namespace panzer_stk {
49 
51 {
53 }
54 
57 {
58 }
59 
61 Teuchos::RCP<STK_Interface> MultiBlockMeshFactory::buildMesh(stk::ParallelMachine parallelMach) const
62 {
63  // build all meta data
64  RCP<STK_Interface> mesh = buildUncommitedMesh(parallelMach);
65 
66  // commit meta data
67  mesh->initialize(parallelMach);
68 
69  // build bulk data
70  completeMeshConstruction(*mesh,parallelMach);
71 
72  return mesh;
73 }
74 
76 {
77  RCP<STK_Interface> mesh = rcp(new STK_Interface(2));
78 
79  machRank_ = stk::parallel_machine_rank(parallelMach);
80  machSize_ = stk::parallel_machine_size(parallelMach);
81 
82  // build meta information: blocks and side set setups
83  buildMetaData(parallelMach,*mesh);
84 
85  return mesh;
86 }
87 
88 void MultiBlockMeshFactory::completeMeshConstruction(STK_Interface & mesh,stk::ParallelMachine parallelMach) const
89 {
90  if(not mesh.isInitialized())
91  mesh.initialize(parallelMach);
92 
93  // add node and element information
94  buildElements(parallelMach,mesh);
95 
96  // finish up the edges
97  mesh.buildSubcells();
98  mesh.buildLocalElementIDs();
99 
100  // now that edges are built, sidets can be added
101  addSideSets(mesh);
102 
103  // calls Stk_MeshFactory::rebalance
104  this->rebalance(mesh);
105 }
106 
109 {
110 }
111 
114 {
115  static RCP<Teuchos::ParameterList> defaultParams;
116 
117  // fill with default values
118  if(defaultParams == Teuchos::null) {
119  defaultParams = rcp(new Teuchos::ParameterList);
120 
121  defaultParams->set<int>("X Elements",2);
122  defaultParams->set<int>("Y Elements",2);
123  }
124 
125  return defaultParams;
126 }
127 
129 {
130  nXElems_ = 2;
131  nYElems_ = 2;
132 }
133 
134 void MultiBlockMeshFactory::buildMetaData(stk::ParallelMachine /* parallelMach */, STK_Interface& mesh) const
135 {
136  typedef shards::Quadrilateral<4> QuadTopo;
137  const CellTopologyData * ctd = shards::getCellTopologyData<QuadTopo>();
138  const CellTopologyData * side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(1,0);
139 
140  // build meta data
141  //mesh.setDimension(2);
142  for(int bx=0;bx<2;bx++) {
143  // add this element block
144  std::stringstream ebPostfix;
145  ebPostfix << "-" << "0_" << bx;
146 
147  // add element blocks
148  mesh.addElementBlock("eblock"+ebPostfix.str(),ctd);
149  }
150 
151  // add sidesets
152  mesh.addSideset("left",side_ctd);
153  mesh.addSideset("right",side_ctd);
154  mesh.addSideset("top",side_ctd);
155  mesh.addSideset("bottom",side_ctd);
156 }
157 
158 void MultiBlockMeshFactory::buildElements(stk::ParallelMachine parallelMach,STK_Interface & mesh) const
159 {
160  mesh.beginModification();
161 
162  if(machRank_==0) {
163  buildBlock(parallelMach,0,0,mesh);
164  }
165  else if(machRank_==1) {
166  buildBlock(parallelMach,0,1,mesh);
167  }
168  else TEUCHOS_ASSERT(false);
169 
170  mesh.endModification();
171 }
172 
173 void MultiBlockMeshFactory::buildBlock(stk::ParallelMachine /* parallelMach */, int xBlock, int yBlock, STK_Interface& mesh) const
174 {
175  int myXElems_start = (machRank_==0 ? 0 : 2);
176  int myXElems_end = (machRank_==0 ? 1 : 3);
177  int myYElems_start = 0;
178  int myYElems_end = 1;
179  int totalXElems = nXElems_*2;
180  int totalYElems = nYElems_*1;
181 
182  double x0_ = 0.0;
183  double xf_ = 1.0;
184  double y0_ = 0.0;
185  double yf_ = 1.0;
186  double deltaX = (xf_-x0_)/double(totalXElems);
187  double deltaY = (yf_-y0_)/double(totalYElems);
188 
189  std::vector<double> coord(2,0.0);
190 
191  // build the nodes
192  for(int nx=myXElems_start;nx<myXElems_end+1;++nx) {
193  coord[0] = double(nx)*deltaX+x0_;
194  for(int ny=myYElems_start;ny<myYElems_end+1;++ny) {
195  coord[1] = double(ny)*deltaY+y0_;
196 
197  mesh.addNode(ny*(totalXElems+1)+nx+1,coord);
198  }
199  }
200 
201  std::stringstream blockName;
202  blockName << "eblock-" << xBlock << "_" << yBlock;
203  stk::mesh::Part * block = mesh.getElementBlockPart(blockName.str());
204 
205  // build the elements
206  for(int nx=myXElems_start;nx<myXElems_end;++nx) {
207  for(int ny=myYElems_start;ny<myYElems_end;++ny) {
208  stk::mesh::EntityId gid = totalXElems*ny+nx+1;
209  std::vector<stk::mesh::EntityId> nodes(4);
210  nodes[0] = nx+1+ny*(totalXElems+1);
211  nodes[1] = nodes[0]+1;
212  nodes[2] = nodes[1]+(totalXElems+1);
213  nodes[3] = nodes[2]-1;
214 
215  RCP<ElementDescriptor> ed = rcp(new ElementDescriptor(gid,nodes));
216  mesh.addElement(ed,block);
217  }
218  }
219 }
220 
221 std::pair<int,int> MultiBlockMeshFactory::determineXElemSizeAndStart(int xBlock,unsigned int size,unsigned int rank) const
222 {
223  unsigned int minElements = nXElems_/size;
224  unsigned int extra = nXElems_ - minElements*size;
225 
226  TEUCHOS_ASSERT(minElements>0);
227 
228  // first "extra" elements get an extra column of elements
229  // this determines the starting X index and number of elements
230  int nume=0, start=0;
231  if(rank<extra) {
232  nume = minElements+1;
233  start = rank*(minElements+1);
234  }
235  else {
236  nume = minElements;
237  start = extra*(minElements+1)+(rank-extra)*minElements;
238  }
239 
240  return std::make_pair(start+nXElems_*xBlock,nume);
241 }
242 
243 std::pair<int,int> MultiBlockMeshFactory::determineYElemSizeAndStart(int yBlock, unsigned int /* size */, unsigned int /* rank */) const
244 {
245  int start = yBlock*nYElems_;
246 
247  return std::make_pair(start,nYElems_);
248 }
249 
251 {
252  mesh.beginModification();
253 
254  std::size_t totalXElems = nXElems_*2;
255  std::size_t totalYElems = nYElems_*1;
256 
257  // get all part vectors
258  stk::mesh::Part * left = mesh.getSideset("left");
259  stk::mesh::Part * right = mesh.getSideset("right");
260  stk::mesh::Part * top = mesh.getSideset("top");
261  stk::mesh::Part * bottom = mesh.getSideset("bottom");
262 
263  std::vector<stk::mesh::Entity> localElmts;
264  mesh.getMyElements(localElmts);
265 
266  // loop over elements adding edges to sidesets
267  std::vector<stk::mesh::Entity>::const_iterator itr;
268  for(itr=localElmts.begin();itr!=localElmts.end();++itr) {
269  stk::mesh::Entity element = (*itr);
270  stk::mesh::EntityId gid = mesh.elementGlobalId(element);
271 
272  std::size_t nx,ny;
273  ny = (gid-1) / totalXElems;
274  nx = gid-ny*totalXElems-1;
275 
276  // vertical boundaries
278 
279  if(nx+1==totalXElems) {
280  stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 1);
281 
282  // on the right
283  if(mesh.entityOwnerRank(edge)==machRank_)
284  mesh.addEntityToSideset(edge,right);
285  }
286 
287  if(nx==0) {
288  stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 3);
289 
290  // on the left
291  if(mesh.entityOwnerRank(edge)==machRank_)
292  mesh.addEntityToSideset(edge,left);
293  }
294 
295  // horizontal boundaries
297 
298  if(ny==0) {
299  stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 0);
300 
301  // on the bottom
302  if(mesh.entityOwnerRank(edge)==machRank_)
303  mesh.addEntityToSideset(edge,bottom);
304  }
305 
306  if(ny+1==totalYElems) {
307  stk::mesh::Entity edge = mesh.findConnectivityById(element, stk::topology::EDGE_RANK, 2);
308 
309  // on the top
310  if(mesh.entityOwnerRank(edge)==machRank_)
311  mesh.addEntityToSideset(edge,top);
312  }
313  }
314 
315  mesh.endModification();
316 }
317 
318 } // end panzer_stk
virtual Teuchos::RCP< STK_Interface > buildUncommitedMesh(stk::ParallelMachine parallelMach) const
unsigned entityOwnerRank(stk::mesh::Entity entity) const
stk::mesh::Part * getElementBlockPart(const std::string &name) const
get the block count
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
void addEntityToSideset(stk::mesh::Entity entity, stk::mesh::Part *sideset)
stk::mesh::Part * getSideset(const std::string &name) const
void buildBlock(stk::ParallelMachine parallelMach, int xBlock, int yBlock, STK_Interface &mesh) 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)
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true, const bool buildRefinementSupport=false)
void rebalance(STK_Interface &mesh) const
stk::mesh::EntityId elementGlobalId(std::size_t lid) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
From ParameterListAcceptor.
Teuchos::RCP< STK_Interface > buildMesh(stk::ParallelMachine parallelMach) const
Build the mesh object.
void addNode(stk::mesh::EntityId gid, const std::vector< double > &coord)
void buildSubcells()
force the mesh to build subcells: edges and faces
bool isInitialized() const
Has initialize been called on this mesh object?
void buildElements(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
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::pair< int, int > determineYElemSizeAndStart(int yBlock, unsigned int size, unsigned int rank) const
std::pair< int, int > determineXElemSizeAndStart(int xBlock, unsigned int size, unsigned int rank) const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
From ParameterListAcceptor.
virtual void completeMeshConstruction(STK_Interface &mesh, stk::ParallelMachine parallelMach) const
void getMyElements(std::vector< stk::mesh::Entity > &elements) const
#define TEUCHOS_ASSERT(assertion_test)
void buildMetaData(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
void addElementBlock(const std::string &name, const CellTopologyData *ctData)