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  mesh->setBoundingBoxSearchFlag(useBBoxSearch_);
95 
96  return mesh;
97  }
98 
99  void
101  stk::ParallelMachine parallelMach) const
102  {
103  PANZER_FUNC_TIME_MONITOR("panzer::CustomMeshFactory::completeMeshConstruction()");
104 
105  if (not mesh.isInitialized())
106  mesh.initialize(parallelMach);
107 
108  // add node and element information
109  buildElements(mesh);
110 
111  // build edges and faces; fyi: addSides(mesh) builds only edges
112  mesh.buildSubcells();
113  mesh.buildLocalElementIDs();
114 
115 
116 
117  // now that edges are built, side and node sets can be added
118  addSideSets(mesh);
119 
120  // set solution fields
121  fillSolutionFieldData(mesh);
122 
123  // calls Stk_MeshFactory::rebalance
124  //this->rebalance(mesh);
125  }
126 
128  void
130  {
132 
133  setMyParamList(paramList);
134 
135  Dimension_ = paramList->get<int>("Dimension");
136 
137  NumBlocks_ = paramList->get<int>("NumBlocks");
138 
139  NumNodesPerProc_ = paramList->get<int>("NumNodesPerProc");
140  Nodes_ = paramList->get<int*>("Nodes");
141 
142  Coords_ = paramList->get<double*>("Coords");
143 
144  NumElementsPerProc_ = paramList->get<int>("NumElementsPerProc");
145 
146  BlockIDs_ = paramList->get<int*>("BlockIDs");
147  Element2Nodes_ = paramList->get<int*>("Element2Nodes");
148 
149  OffsetToGlobalElementIDs_ = paramList->get<int>("OffsetToGlobalElementIDs");
150 
151  // solution fields from tramonto
152  ChargeDensity_ = paramList->get<double*>("ChargeDensity");
153  ElectricPotential_ = paramList->get<double*>("ElectricPotential");
154 
155  // read in periodic boundary conditions
156  parsePeriodicBCList(Teuchos::rcpFromRef(paramList->sublist("Periodic BCs")),periodicBCVec_,useBBoxSearch_);
157  }
158 
162  {
163  static RCP<Teuchos::ParameterList> defaultParams;
164 
165  // fill with default values
166  if(defaultParams == Teuchos::null) {
167  defaultParams = rcp(new Teuchos::ParameterList);
168 
169  defaultParams->set<int>("Dimension",3);
170 
171  defaultParams->set<int>("NumBlocks",0);
172 
173  defaultParams->set<int>("NumNodesPerProc",0);
174  defaultParams->set<int*>("Nodes",NULL);
175 
176  defaultParams->set<double*>("Coords",NULL);
177 
178  defaultParams->set<int>("NumElementsPerProc",0);
179 
180  defaultParams->set<int*>("BlockIDs",NULL);
181  defaultParams->set<int*>("Element2Nodes",NULL);
182 
183  defaultParams->set<int>("OffsetToGlobalElementIDs", 0);
184 
185  defaultParams->set<double*>("ChargeDensity",NULL);
186  defaultParams->set<double*>("ElectricPotential",NULL);
187 
188  Teuchos::ParameterList &bcs = defaultParams->sublist("Periodic BCs");
189  bcs.set<int>("Count",0); // no default periodic boundary conditions
190  }
191 
192  return defaultParams;
193  }
194 
195  void
197  {
198  // get valid parameters
200 
201  // set that parameter list
202  setParameterList(validParams);
203  }
204 
205  void
207  {
208  typedef shards::Hexahedron<8> HexTopo;
209  const CellTopologyData *ctd = shards::getCellTopologyData<HexTopo>();
210  const CellTopologyData *side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(2,0);
211 
212  for (int blk=0;blk<NumBlocks_;++blk) {
213  std::stringstream block_id;
214  block_id << "eblock-" << blk;
215 
216  // add element blocks
217  mesh.addElementBlock(block_id.str(),ctd);
218 
219  mesh.addSolutionField("CHARGE_DENSITY",block_id.str());
220  mesh.addSolutionField("ELECTRIC_POTENTIAL",block_id.str());
221  }
222 
223  mesh.addSideset("left", side_ctd);
224  mesh.addSideset("right", side_ctd);
225  mesh.addSideset("top", side_ctd);
226  mesh.addSideset("bottom",side_ctd);
227  mesh.addSideset("front", side_ctd);
228  mesh.addSideset("back", side_ctd);
229 
230  mesh.addSideset("wall", side_ctd);
231  }
232 
233  void
235  {
236  mesh.beginModification();
237 
238  const int dim = mesh.getDimension();
239 
240  // build the nodes
241  std::vector<double> coords(dim,0.0);
242  for (int i=0;i<NumNodesPerProc_;++i) {
243  for (int k=0;k<dim;++k)
244  coords[k] = Coords_[i*dim+k];
245  mesh.addNode(Nodes_[i], coords);
246  }
247 
248  // build the elements
249  std::vector<std::string> block_ids;
250  mesh.getElementBlockNames(block_ids);
251 
252  for (int i=0;i<NumElementsPerProc_;++i) {
253 
254  // get block by its name
255  std::stringstream block_id;
256  block_id << "eblock-" << BlockIDs_[i];
257 
258  stk::mesh::Part *block = mesh.getElementBlockPart(block_id.str());
259 
260  // construct element and its nodal connectivity
261  stk::mesh::EntityId elt = i + OffsetToGlobalElementIDs_;
262  std::vector<stk::mesh::EntityId> elt2nodes(8);
263 
264  for (int k=0;k<8;++k)
265  elt2nodes[k] = Element2Nodes_[i*8+k];
266 
267  RCP<ElementDescriptor> ed = rcp(new ElementDescriptor(elt,elt2nodes));
268  mesh.addElement(ed,block);
269  }
270 
271  mesh.endModification();
272  }
273 
274  void
276  {
277  mesh.beginModification();
278  stk::mesh::BulkData& bulkData = *mesh.getBulkData();
279  const stk::mesh::EntityRank sideRank = mesh.getSideRank();
280 
281  // get all part vectors
282  stk::mesh::Part *box[6];
283  box[0] = mesh.getSideset("front");
284  box[1] = mesh.getSideset("right");
285  box[2] = mesh.getSideset("back");
286  box[3] = mesh.getSideset("left");
287  box[4] = mesh.getSideset("bottom");
288  box[5] = mesh.getSideset("top");
289 
290  stk::mesh::Part *wall = mesh.getSideset("wall");
291 
292  std::vector<stk::mesh::Entity> elements;
293  mesh.getMyElements(elements);
294 
295  // loop over elements adding sides to sidesets
296  for (std::vector<stk::mesh::Entity>::const_iterator
297  itr=elements.begin();itr!=elements.end();++itr) {
298  stk::mesh::Entity element = (*itr);
299  const size_t numSides = bulkData.num_connectivity(element, sideRank);
300  stk::mesh::Entity const* relations = bulkData.begin(element, sideRank);
301 
302  // loop over side id checking element neighbors
303  for (std::size_t i=0;i<numSides;++i) {
304  stk::mesh::Entity side = relations[i];
305  const size_t numNeighbors = bulkData.num_elements(side);
306  stk::mesh::Entity const* neighbors = bulkData.begin_elements(side);
307 
308  if (numNeighbors == 1) {
309  if (mesh.entityOwnerRank(side) == machRank_)
310  mesh.addEntityToSideset(side, box[i]);
311  }
312  else if (numNeighbors == 2) {
313  std::string neig_block_id_0 = mesh.containingBlockId(neighbors[0]);
314  std::string neig_block_id_1 = mesh.containingBlockId(neighbors[1]);
315  if (neig_block_id_0 != neig_block_id_1 && mesh.entityOwnerRank(side) == machRank_)
316  mesh.addEntityToSideset(side, wall);
317  }
318  else {
319  // runtime exception
320  }
321  }
322  }
323 
324  mesh.endModification();
325  }
326 
327  void
329  {
330  constexpr std::size_t dim_1 = 8;
331 
332  for (int blk=0;blk<NumBlocks_;++blk) {
333 
334  std::stringstream block_id;
335  block_id << "eblock-" << blk;
336 
337  // elements in this processor for this block
338  std::vector<stk::mesh::Entity> elements;
339  mesh.getMyElements(block_id.str(), elements);
340 
341  // size of elements in the current block
342  const auto n_elements = elements.size();
343 
344  // build local element index
345  std::vector<std::size_t> local_ids;
346 
347  Kokkos::View<double**, Kokkos::HostSpace> charge_density_by_local_ids("charge_density_by_local_ids", n_elements, dim_1);
348  Kokkos::View<double**, Kokkos::HostSpace> electric_potential_by_local_ids("electric_potential_by_local_ids", n_elements, dim_1);
349 
350  for (const auto& elem : elements){
351  local_ids.push_back(mesh.elementLocalId(elem));
352  const auto q = mesh.elementGlobalId(elem) - OffsetToGlobalElementIDs_;
353 
354  for (std::size_t k=0;k<dim_1;++k) {
355  const auto loc = q*dim_1 + k;
356  charge_density_by_local_ids(q,k) = ChargeDensity_[loc];
357  electric_potential_by_local_ids(q,k) = ElectricPotential_[loc];
358  }
359  }
360 
361  // write out to stk mesh
362  mesh.setSolutionFieldData("CHARGE_DENSITY",
363  block_id.str(),
364  local_ids,
365  charge_density_by_local_ids, 1.0);
366 
367  mesh.setSolutionFieldData("ELECTRIC_POTENTIAL",
368  block_id.str(),
369  local_ids,
370  electric_potential_by_local_ids, 1.0);
371  }
372  }
373 } // 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.
static void parsePeriodicBCList(const Teuchos::RCP< Teuchos::ParameterList > &pl, std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > &periodicBC, bool &useBBoxSearch)
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 part
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)
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
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true, const bool buildRefinementSupport=false)
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)