Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_STK_CustomMeshFactory.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 
47 using Teuchos::RCP;
48 using Teuchos::rcp;
49 
50 namespace panzer_stk {
51 
53  {
55  }
56 
59  {
60  }
61 
64  CustomMeshFactory::buildMesh(stk::ParallelMachine parallelMach) const
65  {
66  PANZER_FUNC_TIME_MONITOR("panzer::CustomMeshFactory::buildMesh()");
67 
68  // build all meta data
69  RCP<STK_Interface> mesh = buildUncommitedMesh(parallelMach);
70 
71  // commit meta data; allocate field variables
72  mesh->initialize(parallelMach);
73 
74  // build bulk data
75  completeMeshConstruction(*mesh,parallelMach);
76 
77  return mesh;
78  }
79 
81  CustomMeshFactory::buildUncommitedMesh(stk::ParallelMachine parallelMach) const
82  {
83  PANZER_FUNC_TIME_MONITOR("panzer::CustomMeshFactory::buildUncomittedMesh()");
84 
85  RCP<STK_Interface> mesh = rcp(new STK_Interface(3));
86 
87  machRank_ = stk::parallel_machine_rank(parallelMach);
88  machSize_ = stk::parallel_machine_size(parallelMach);
89 
90  // add blocks and side sets (global geometry setup)
91  buildMetaData(*mesh);
92 
93  mesh->addPeriodicBCs(periodicBCVec_);
94 
95  return mesh;
96  }
97 
98  void
100  stk::ParallelMachine parallelMach) const
101  {
102  PANZER_FUNC_TIME_MONITOR("panzer::CustomMeshFactory::completeMeshConstruction()");
103 
104  if (not mesh.isInitialized())
105  mesh.initialize(parallelMach);
106 
107  // add node and element information
108  buildElements(mesh);
109 
110  // build edges and faces; fyi: addSides(mesh) builds only edges
111  mesh.buildSubcells();
112  mesh.buildLocalElementIDs();
113 
114 
115 
116  // now that edges are built, side and node sets can be added
117  addSideSets(mesh);
118 
119  // set solution fields
120  fillSolutionFieldData(mesh);
121 
122  // calls Stk_MeshFactory::rebalance
123  //this->rebalance(mesh);
124  }
125 
127  void
129  {
131 
132  setMyParamList(paramList);
133 
134  Dimension_ = paramList->get<int>("Dimension");
135 
136  NumBlocks_ = paramList->get<int>("NumBlocks");
137 
138  NumNodesPerProc_ = paramList->get<int>("NumNodesPerProc");
139  Nodes_ = paramList->get<int*>("Nodes");
140 
141  Coords_ = paramList->get<double*>("Coords");
142 
143  NumElementsPerProc_ = paramList->get<int>("NumElementsPerProc");
144 
145  BlockIDs_ = paramList->get<int*>("BlockIDs");
146  Element2Nodes_ = paramList->get<int*>("Element2Nodes");
147 
148  OffsetToGlobalElementIDs_ = paramList->get<int>("OffsetToGlobalElementIDs");
149 
150  // solution fields from tramonto
151  ChargeDensity_ = paramList->get<double*>("ChargeDensity");
152  ElectricPotential_ = paramList->get<double*>("ElectricPotential");
153 
154  // read in periodic boundary conditions
155  parsePeriodicBCList(Teuchos::rcpFromRef(paramList->sublist("Periodic BCs")),periodicBCVec_);
156  }
157 
161  {
162  static RCP<Teuchos::ParameterList> defaultParams;
163 
164  // fill with default values
165  if(defaultParams == Teuchos::null) {
166  defaultParams = rcp(new Teuchos::ParameterList);
167 
168  defaultParams->set<int>("Dimension",3);
169 
170  defaultParams->set<int>("NumBlocks",0);
171 
172  defaultParams->set<int>("NumNodesPerProc",0);
173  defaultParams->set<int*>("Nodes",NULL);
174 
175  defaultParams->set<double*>("Coords",NULL);
176 
177  defaultParams->set<int>("NumElementsPerProc",0);
178 
179  defaultParams->set<int*>("BlockIDs",NULL);
180  defaultParams->set<int*>("Element2Nodes",NULL);
181 
182  defaultParams->set<int>("OffsetToGlobalElementIDs", 0);
183 
184  defaultParams->set<double*>("ChargeDensity",NULL);
185  defaultParams->set<double*>("ElectricPotential",NULL);
186 
187  Teuchos::ParameterList &bcs = defaultParams->sublist("Periodic BCs");
188  bcs.set<int>("Count",0); // no default periodic boundary conditions
189  }
190 
191  return defaultParams;
192  }
193 
194  void
196  {
197  // get valid parameters
199 
200  // set that parameter list
201  setParameterList(validParams);
202  }
203 
204  void
206  {
207  typedef shards::Hexahedron<8> HexTopo;
208  const CellTopologyData *ctd = shards::getCellTopologyData<HexTopo>();
209  const CellTopologyData *side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(2,0);
210 
211  for (int blk=0;blk<NumBlocks_;++blk) {
212  std::stringstream block_id;
213  block_id << "eblock-" << blk;
214 
215  // add element blocks
216  mesh.addElementBlock(block_id.str(),ctd);
217 
218  mesh.addSolutionField("CHARGE_DENSITY",block_id.str());
219  mesh.addSolutionField("ELECTRIC_POTENTIAL",block_id.str());
220  }
221 
222  mesh.addSideset("left", side_ctd);
223  mesh.addSideset("right", side_ctd);
224  mesh.addSideset("top", side_ctd);
225  mesh.addSideset("bottom",side_ctd);
226  mesh.addSideset("front", side_ctd);
227  mesh.addSideset("back", side_ctd);
228 
229  mesh.addSideset("wall", side_ctd);
230  }
231 
232  void
234  {
235  mesh.beginModification();
236 
237  const int dim = mesh.getDimension();
238 
239  // build the nodes
240  std::vector<double> coords(dim,0.0);
241  for (int i=0;i<NumNodesPerProc_;++i) {
242  for (int k=0;k<dim;++k)
243  coords[k] = Coords_[i*dim+k];
244  mesh.addNode(Nodes_[i], coords);
245  }
246 
247  // build the elements
248  std::vector<std::string> block_ids;
249  mesh.getElementBlockNames(block_ids);
250 
251  for (int i=0;i<NumElementsPerProc_;++i) {
252 
253  // get block by its name
254  std::stringstream block_id;
255  block_id << "eblock-" << BlockIDs_[i];
256 
257  stk::mesh::Part *block = mesh.getElementBlockPart(block_id.str());
258 
259  // construct element and its nodal connectivity
260  stk::mesh::EntityId elt = i + OffsetToGlobalElementIDs_;
261  std::vector<stk::mesh::EntityId> elt2nodes(8);
262 
263  for (int k=0;k<8;++k)
264  elt2nodes[k] = Element2Nodes_[i*8+k];
265 
266  RCP<ElementDescriptor> ed = rcp(new ElementDescriptor(elt,elt2nodes));
267  mesh.addElement(ed,block);
268  }
269 
270  mesh.endModification();
271  }
272 
273  void
275  {
276  mesh.beginModification();
277  stk::mesh::BulkData& bulkData = *mesh.getBulkData();
278  const stk::mesh::EntityRank sideRank = mesh.getSideRank();
279 
280  // get all part vectors
281  stk::mesh::Part *box[6];
282  box[0] = mesh.getSideset("front");
283  box[1] = mesh.getSideset("right");
284  box[2] = mesh.getSideset("back");
285  box[3] = mesh.getSideset("left");
286  box[4] = mesh.getSideset("bottom");
287  box[5] = mesh.getSideset("top");
288 
289  stk::mesh::Part *wall = mesh.getSideset("wall");
290 
291  std::vector<stk::mesh::Entity> elements;
292  mesh.getMyElements(elements);
293 
294  // loop over elements adding sides to sidesets
295  for (std::vector<stk::mesh::Entity>::const_iterator
296  itr=elements.begin();itr!=elements.end();++itr) {
297  stk::mesh::Entity element = (*itr);
298  const size_t numSides = bulkData.num_connectivity(element, sideRank);
299  stk::mesh::Entity const* relations = bulkData.begin(element, sideRank);
300 
301  // loop over side id checking element neighbors
302  for (std::size_t i=0;i<numSides;++i) {
303  stk::mesh::Entity side = relations[i];
304  const size_t numNeighbors = bulkData.num_elements(side);
305  stk::mesh::Entity const* neighbors = bulkData.begin_elements(side);
306 
307  if (numNeighbors == 1) {
308  if (mesh.entityOwnerRank(side) == machRank_)
309  mesh.addEntityToSideset(side, box[i]);
310  }
311  else if (numNeighbors == 2) {
312  std::string neig_block_id_0 = mesh.containingBlockId(neighbors[0]);
313  std::string neig_block_id_1 = mesh.containingBlockId(neighbors[1]);
314  if (neig_block_id_0 != neig_block_id_1 && mesh.entityOwnerRank(side) == machRank_)
315  mesh.addEntityToSideset(side, wall);
316  }
317  else {
318  // runtime exception
319  }
320  }
321  }
322 
323  mesh.endModification();
324  }
325 
326  void
328  {
329  for (int blk=0;blk<NumBlocks_;++blk) {
330 
331  std::stringstream block_id;
332  block_id << "eblock-" << blk;
333 
334  // elements in this processor for this block
335  std::vector<stk::mesh::Entity> elements;
336  mesh.getMyElements(block_id.str(), elements);
337 
338  // size of elements in the current block
339  std::size_t n_elements = elements.size();
340 
341  // build local element index
342  std::vector<std::size_t> local_ids;
343  for (std::vector<stk::mesh::Entity>::const_iterator
344  itr=elements.begin();itr!=elements.end();++itr)
345  local_ids.push_back(mesh.elementLocalId(*itr));
346 
347  // re-index solution fields in the same order of local_ids
348  std::vector<double> charge_density_by_local_ids, electric_potential_by_local_ids;
349  for (std::vector<stk::mesh::Entity>::const_iterator
350  itr=elements.begin();itr!=elements.end();++itr) {
351  int q = mesh.elementGlobalId(*itr) - OffsetToGlobalElementIDs_;
352  for (int k=0;k<8;++k) {
353  int loc = q*8 + k;
354  charge_density_by_local_ids.push_back(ChargeDensity_[loc]);
355  electric_potential_by_local_ids.push_back(ElectricPotential_[loc]);
356  }
357  }
358 
359  // wrap the buffer with a proper container
360  FieldContainer charge_density(n_elements, 8, &charge_density_by_local_ids[0]),
361  electric_potential(n_elements, 8, &electric_potential_by_local_ids[0]);
362 
363  // write out to stk mesh
364  mesh.setSolutionFieldData("CHARGE_DENSITY",
365  block_id.str(),
366  local_ids,
367  charge_density, 1.0);
368 
369  mesh.setSolutionFieldData("ELECTRIC_POTENTIAL",
370  block_id.str(),
371  local_ids,
372  electric_potential, 1.0);
373  }
374  }
375 } // end panzer_stk
unsigned entityOwnerRank(stk::mesh::Entity entity) const
void setSolutionFieldData(const std::string &fieldName, const std::string &blockId, const std::vector< std::size_t > &localElementIds, const ArrayT &solutionValues, double scaleValue=1.0)
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
From ParameterListAcceptor.
void getElementBlockNames(std::vector< std::string > &names) const
T & get(const std::string &name, T def_value)
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)
void addSolutionField(const std::string &fieldName, const std::string &blockId)
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true)
stk::mesh::Part * getSideset(const std::string &name) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
From ParameterListAcceptor.
Teuchos::RCP< STK_Interface > buildMesh(stk::ParallelMachine parallelMach) const
Build the mesh object.
std::size_t elementLocalId(stk::mesh::Entity elmt) const
void addSideSets(STK_Interface &mesh) const
unsigned getDimension() const
get the dimension
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
virtual void completeMeshConstruction(STK_Interface &mesh, stk::ParallelMachine parallelMach) const
Teuchos::RCP< stk::mesh::BulkData > getBulkData() 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)
std::string containingBlockId(stk::mesh::Entity elmt) const
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 fillSolutionFieldData(STK_Interface &mesh) const
void buildElements(STK_Interface &mesh) 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
virtual Teuchos::RCP< STK_Interface > buildUncommitedMesh(stk::ParallelMachine parallelMach) const
void buildMetaData(STK_Interface &mesh) 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)