Panzer Version of the Day
Loading...
Searching...
No Matches
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
47using Teuchos::RCP;
48using Teuchos::rcp;
49
50namespace panzer_stk {
51
56
61
63 Teuchos::RCP<STK_Interface>
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
80 Teuchos::RCP<STK_Interface>
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();
114
115
116
117 // now that edges are built, side and node sets can be added
118 addSideSets(mesh);
119
120 // set solution fields
122
123 // calls Stk_MeshFactory::rebalance
124 //this->rebalance(mesh);
125 }
126
128 void
129 CustomMeshFactory::setParameterList(const Teuchos::RCP<Teuchos::ParameterList> &paramList)
130 {
131 paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
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
160 Teuchos::RCP<const Teuchos::ParameterList>
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
199 RCP<Teuchos::ParameterList> validParams = rcp(new Teuchos::ParameterList(*getValidParameters()));
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
void buildMetaData(STK_Interface &mesh) const
void addSideSets(STK_Interface &mesh) const
void fillSolutionFieldData(STK_Interface &mesh) const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
From ParameterListAcceptor.
virtual void completeMeshConstruction(STK_Interface &mesh, stk::ParallelMachine parallelMach) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
From ParameterListAcceptor.
virtual Teuchos::RCP< STK_Interface > buildUncommitedMesh(stk::ParallelMachine parallelMach) const
void buildElements(STK_Interface &mesh) const
Teuchos::RCP< STK_Interface > buildMesh(stk::ParallelMachine parallelMach) const
Build the mesh object.
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
stk::mesh::EntityId elementGlobalId(std::size_t lid) const
bool isInitialized() const
Has initialize been called on this mesh object?
std::string containingBlockId(stk::mesh::Entity elmt) const
void buildSubcells()
force the mesh to build subcells: edges and faces
void addElement(const Teuchos::RCP< ElementDescriptor > &ed, stk::mesh::Part *block)
void addNode(stk::mesh::EntityId gid, const std::vector< double > &coord)
void setSolutionFieldData(const std::string &fieldName, const std::string &blockId, const std::vector< std::size_t > &localElementIds, const ArrayT &solutionValues, double scaleValue=1.0)
std::size_t elementLocalId(stk::mesh::Entity elmt) const
void getElementBlockNames(std::vector< std::string > &names) const
void addSideset(const std::string &name, const CellTopologyData *ctData)
unsigned getDimension() const
get the dimension
unsigned entityOwnerRank(stk::mesh::Entity entity) const
void addSolutionField(const std::string &fieldName, const std::string &blockId)
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)
Teuchos::RCP< stk::mesh::BulkData > getBulkData() const
stk::mesh::Part * getSideset(const std::string &name) const
static void parsePeriodicBCList(const Teuchos::RCP< Teuchos::ParameterList > &pl, std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > &periodicBC, bool &useBBoxSearch)
std::vector< Teuchos::RCP< const PeriodicBC_MatcherBase > > periodicBCVec_