Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_STK_Quad8ToQuad4MeshFactory.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
43#include "PanzerAdaptersSTK_config.hpp"
47#include "Teuchos_TimeMonitor.hpp"
48#include "Teuchos_StandardParameterEntryValidators.hpp" // for plist validation
49#include <stk_mesh/base/FieldBase.hpp>
50
51// #define ENABLE_UNIFORM
52
53using Teuchos::RCP;
54using Teuchos::rcp;
55
56namespace panzer_stk {
57
58Quad8ToQuad4MeshFactory::Quad8ToQuad4MeshFactory(const std::string& quad8MeshFileName,
59 stk::ParallelMachine mpi_comm,
60 const bool print_debug)
61 : createEdgeBlocks_(false),
62 print_debug_(print_debug)
63{
64 panzer_stk::STK_ExodusReaderFactory factory;
65 RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList);
66 pl->set("File Name",quad8MeshFileName);
67 factory.setParameterList(pl);
68 quad8Mesh_ = factory.buildMesh(mpi_comm);
69
71}
72
73Quad8ToQuad4MeshFactory::Quad8ToQuad4MeshFactory(const Teuchos::RCP<panzer_stk::STK_Interface>& quad8Mesh,
74 const bool print_debug)
75 : quad8Mesh_(quad8Mesh),
76 createEdgeBlocks_(false),
77 print_debug_(print_debug)
78{
80}
81
83Teuchos::RCP<STK_Interface> Quad8ToQuad4MeshFactory::buildMesh(stk::ParallelMachine parallelMach) const
84{
85 PANZER_FUNC_TIME_MONITOR("panzer::Quad8ToQuad4MeshFactory::buildMesh()");
86
87 // Make sure the Quad8 and Quad4 Comms match
88 {
89 int result = MPI_UNEQUAL;
90 MPI_Comm_compare(parallelMach, quad8Mesh_->getBulkData()->parallel(), &result);
91 TEUCHOS_ASSERT(result != MPI_UNEQUAL);
92 }
93
94 // build all meta data
95 RCP<STK_Interface> mesh = this->buildUncommitedMesh(parallelMach);
96
97 // commit meta data
98 mesh->initialize(parallelMach);
99
100 // build bulk data
101 this->completeMeshConstruction(*mesh,parallelMach);
102
103 return mesh;
104}
105
106Teuchos::RCP<STK_Interface> Quad8ToQuad4MeshFactory::buildUncommitedMesh(stk::ParallelMachine parallelMach) const
107{
108 PANZER_FUNC_TIME_MONITOR("panzer::Quad8ToQuad4MeshFactory::buildUncomittedMesh()");
109
110 RCP<STK_Interface> mesh = rcp(new STK_Interface(2));
111
112 machRank_ = stk::parallel_machine_rank(parallelMach);
113 machSize_ = stk::parallel_machine_size(parallelMach);
114
115 // build meta information: blocks and side set setups
116 this->buildMetaData(parallelMach,*mesh);
117
118 mesh->addPeriodicBCs(periodicBCVec_);
119 mesh->setBoundingBoxSearchFlag(useBBoxSearch_);
120
121 return mesh;
122}
123
124void Quad8ToQuad4MeshFactory::completeMeshConstruction(STK_Interface & mesh,stk::ParallelMachine parallelMach) const
125{
126 PANZER_FUNC_TIME_MONITOR("panzer::Quad8ToQuad4MeshFactory::completeMeshConstruction()");
127
128 if(not mesh.isInitialized())
129 mesh.initialize(parallelMach);
130
131 // add node and element information
132 this->buildElements(parallelMach,mesh);
133
134 // finish up the edges
135#ifndef ENABLE_UNIFORM
136 mesh.buildSubcells();
137#endif
140 mesh.buildLocalEdgeIDs();
141 }
142
143 // now that edges are built, sidsets can be added
144#ifndef ENABLE_UNIFORM
145 this->addSideSets(mesh);
146#endif
147
148 // add nodesets
149 this->addNodeSets(mesh);
150
152 this->addEdgeBlocks(mesh);
153 }
154
155 // Copy field data
156 this->copyCellFieldData(mesh);
157
158 // calls Stk_MeshFactory::rebalance
159 // this->rebalance(mesh);
160}
161
163void Quad8ToQuad4MeshFactory::setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & paramList)
164{
165 paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
166
167 setMyParamList(paramList);
168
169 // offsetGIDs_ = (paramList->get<std::string>("Offset mesh GIDs above 32-bit int limit") == "ON") ? true : false;
170
171 createEdgeBlocks_ = paramList->get<bool>("Create Edge Blocks");
172
173 // read in periodic boundary conditions
174 parsePeriodicBCList(Teuchos::rcpFromRef(paramList->sublist("Periodic BCs")),periodicBCVec_,useBBoxSearch_);
175}
176
178Teuchos::RCP<const Teuchos::ParameterList> Quad8ToQuad4MeshFactory::getValidParameters() const
179{
180 static RCP<Teuchos::ParameterList> defaultParams;
181
182 // fill with default values
183 if(defaultParams == Teuchos::null) {
184 defaultParams = rcp(new Teuchos::ParameterList);
185
186 Teuchos::setStringToIntegralParameter<int>(
187 "Offset mesh GIDs above 32-bit int limit",
188 "OFF",
189 "If 64-bit GIDs are supported, the mesh element and node global indices will start at a value greater than 32-bit limit.",
190 Teuchos::tuple<std::string>("OFF", "ON"),
191 defaultParams.get());
192
193 // default to false for backward compatibility
194 defaultParams->set<bool>("Create Edge Blocks",false,"Create edge blocks in the mesh");
195
196 Teuchos::ParameterList & bcs = defaultParams->sublist("Periodic BCs");
197 bcs.set<int>("Count",0); // no default periodic boundary conditions
198 }
199
200 return defaultParams;
201}
202
203void Quad8ToQuad4MeshFactory::buildMetaData(stk::ParallelMachine /* parallelMach */, STK_Interface & mesh) const
204{
205 if (print_debug_) {
206 std::cout << "\n\n**** DEBUG: begin printing source Quad8 exodus file metadata ****\n";
207 quad8Mesh_->getMetaData()->dump_all_meta_info(std::cout);
208 std::cout << "\n\n**** DEBUG: end printing source Quad8 exodus file metadata ****\n";
209 }
210
211 // Required topologies
212 using QuadTopo = shards::Quadrilateral<4>;
213 const CellTopologyData * ctd = shards::getCellTopologyData<QuadTopo>();
214 const CellTopologyData * side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(1,0);
215
216 // Add in element blocks
217 std::vector<std::string> element_block_names;
218 quad8Mesh_->getElementBlockNames(element_block_names);
219 {
220 for (const auto& n : element_block_names)
221 mesh.addElementBlock(n,ctd);
222 }
223
224 // Add in sidesets
225 {
226 std::vector<std::string> sideset_names;
227 quad8Mesh_->getSidesetNames(sideset_names);
228 for (const auto& n : sideset_names)
229 mesh.addSideset(n,side_ctd);
230 }
231
232 // Add in nodesets
233 {
234 std::vector<std::string> nodeset_names;
235 quad8Mesh_->getNodesetNames(nodeset_names);
236 for (const auto& n : nodeset_names)
237 mesh.addNodeset(n);
238 }
239
241 const CellTopologyData * edge_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(1,0);
242 std::vector<std::string> element_block_names;
243 quad8Mesh_->getElementBlockNames(element_block_names);
244 for (const auto& block_name : element_block_names)
245 mesh.addEdgeBlock(block_name,edgeBlockName_,edge_ctd);
246 }
247
248 // Add in nodal fields
249 {
250 const auto& fields = quad8Mesh_->getMetaData()->get_fields(mesh.getNodeRank());
251 for (const auto& f : fields) {
252 if (print_debug_)
253 std::cout << "Field=" << f->name() << ", rank=" << f->entity_rank() << std::endl;
254
255 // Cull out the coordinate fields. That is a default field in
256 // stk that is automatically created by stk.
257 if (f->name() != "coordinates") {
258 for (const auto& n : element_block_names)
259 mesh.addSolutionField(f->name(),n);
260 }
261 }
262 }
263
264 // Add in element fields
265 {
266 const auto& fields = quad8Mesh_->getMetaData()->get_fields(mesh.getElementRank());
267 for (const auto& f : fields) {
268 if (print_debug_)
269 std::cout << "Add Cell Field=" << f->name() << ", rank=" << f->entity_rank() << std::endl;
270
271 for (const auto& n : element_block_names)
272 mesh.addCellField(f->name(),n);
273 }
274 }
275
276 // NOTE: skipping edge and face fields for now. Can't error out since sidesets count as edge fields.
277 // TEUCHOS_TEST_FOR_EXCEPTION(quad8Mesh_->getMetaData()->get_fields(mesh.getEdgeRank()).size() != 0,std::runtime_error,
278 // "ERROR: the Quad8 mesh contains edge fields\""
279 // << quad8Mesh_->getMetaData()->get_fields(mesh.getEdgeRank())[0]->name()
280 // << "\". Edge fields are not supported in Quad8ToQuad4!");
281
282 if (print_debug_) {
283 std::cout << "\n\n**** DEBUG: begin printing source Quad4 exodus file metadata ****\n";
284 mesh.getMetaData()->dump_all_meta_info(std::cout);
285 std::cout << "\n\n**** DEBUG: end printing source Quad4 exodus file metadata ****\n";
286 }
287}
288
289void Quad8ToQuad4MeshFactory::buildElements(stk::ParallelMachine parallelMach,STK_Interface & mesh) const
290{
291 mesh.beginModification();
292
293 auto metadata = mesh.getMetaData();
294 auto bulkdata = mesh.getBulkData();
295
296 // Loop over element blocks
297 std::vector<std::string> block_names;
298 quad8Mesh_->getElementBlockNames(block_names);
299 for (const auto& block_name : block_names) {
300
301 // Get the elements on this process
302 std::vector<stk::mesh::Entity> elements;
303 quad8Mesh_->getMyElements(block_name,elements);
304
305 if (print_debug_) {
306 std::cout << "*************************************************" << std::endl;
307 std::cout << "block_name=" << block_name << ", num_my_elements=" << elements.size() << std::endl;
308 std::cout << "*************************************************" << std::endl;
309 }
310
311 for (const auto& element : elements) {
312
313 const auto element_gid = quad8Mesh_->getBulkData()->identifier(element);
314
315 if (print_debug_) {
316 std::cout << "rank=" << machRank_
317 << ", block=" << block_name
318 << ", element entity_id=" << element
319 << ", gid=" << element_gid << std::endl;
320 }
321
322 // Register nodes with the mesh
323 std::vector<stk::mesh::EntityId> nodes(4);
324 for (int i=0; i < 4; ++i) {
325 // NOTE: this assumes that the Quad4 topology is nested in
326 // the Quad8 as an extended topology - i.e. the first four
327 // nodes of the quad8 topology are the corresponding quad4
328 // nodes. Shards topologies enforce this throught the concept
329 // of extended topologies.
330 stk::mesh::Entity node_entity = quad8Mesh_->findConnectivityById(element,mesh.getNodeRank(),i);
331 TEUCHOS_ASSERT(node_entity.is_local_offset_valid());
332 const auto node_gid = quad8Mesh_->getBulkData()->identifier(node_entity);
333 const double* node_coords = quad8Mesh_->getNodeCoordinates(node_entity);
334 std::vector<double> coords_vec(2);
335 coords_vec[0] = node_coords[0];
336 coords_vec[1] = node_coords[1];
337 mesh.addNode(node_gid,coords_vec);
338 nodes[i]=node_gid;
339 if (print_debug_) {
340 std::cout << "elem gid=" << quad8Mesh_->getBulkData()->identifier(element)
341 << ", node_gid=" << node_gid << ", (" << coords_vec[0] << "," << coords_vec[1] << ")\n";
342 }
343 }
344
345 // Register element with the element block
346 auto element_descriptor = panzer_stk::buildElementDescriptor(element_gid,nodes);
347 auto element_block_part = mesh.getMetaData()->get_part(block_name);
348 mesh.addElement(element_descriptor,element_block_part);
349 }
350 }
351 mesh.endModification();
352}
353
355{
356 mesh.beginModification();
357
358 // Loop over sidesets
359 std::vector<std::string> sideset_names;
360 quad8Mesh_->getSidesetNames(sideset_names);
361 for (const auto& sideset_name : sideset_names) {
362
363 stk::mesh::Part* sideset_part = mesh.getSideset(sideset_name);
364
365 std::vector<stk::mesh::Entity> q8_sides;
366 quad8Mesh_->getMySides(sideset_name,q8_sides);
367
368 // Loop over edges
369 for (const auto q8_edge : q8_sides) {
370 // The edge numbering scheme uses the element/node gids, so it
371 // should be consistent between the quad8 and quad4 meshes
372 // since we used the same gids. We use this fact to populate
373 // the quad4 sidesets.
374 stk::mesh::EntityId edge_gid = quad8Mesh_->getBulkData()->identifier(q8_edge);
375 stk::mesh::Entity q4_edge = mesh.getBulkData()->get_entity(mesh.getEdgeRank(),edge_gid);
376 mesh.addEntityToSideset(q4_edge,sideset_part);
377 }
378 }
379
380 mesh.endModification();
381}
382
385
387{
388 mesh.beginModification();
389
390 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
391 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
392
393 stk::mesh::Part * edge_block = mesh.getEdgeBlock(edgeBlockName_);
394
395 stk::mesh::Selector owned_block = metaData->locally_owned_part();
396
397 std::vector<stk::mesh::Entity> edges;
398 bulkData->get_entities(mesh.getEdgeRank(), owned_block, edges);
399 mesh.addEntitiesToEdgeBlock(edges, edge_block);
400
401 mesh.endModification();
402}
403
405{
406 // Vector of pointers to field data
407 const auto& fields = q4_mesh.getMetaData()->get_fields(q4_mesh.getElementRank());
408
409 // loop over fields and add the data to the new mesh.
410 for (const auto& field : fields) {
411
412 if (print_debug_)
413 std::cout << "Adding field values for \"" << *field << "\" to the Quad4 mesh!\n";
414
415 auto field_name = field->name();
416
417 // Divide into scalar and vector fields, ignoring any other types
418 // for now.
419 if (field->type_is<double>() &&
420 field_name != "coordinates" &&
421 field_name != "PROC_ID" &&
422 field_name != "LOAD_BAL") {
423
424 // Loop over element blocks
425 std::vector<std::string> block_names;
426 quad8Mesh_->getElementBlockNames(block_names);
427 for (const auto& block : block_names) {
428
429 auto* q4_field = q4_mesh.getCellField(field_name,block);
430 TEUCHOS_ASSERT(q4_field != nullptr);
431 // The q8 mesh doesn't have the field names set up, so a query
432 // fails. Go to stk directly in this case.
433 auto* q8_field = quad8Mesh_->getMetaData()->get_field(quad8Mesh_->getElementRank(),field_name);
434#ifdef PANZER_DEBUG
435 TEUCHOS_ASSERT(q8_field != nullptr);
436#endif
437
438 // Get the elements for this block.
439 std::vector<stk::mesh::Entity> q4_elements;
440 q4_mesh.getMyElements(block,q4_elements);
441 std::vector<stk::mesh::Entity> q8_elements;
442 quad8Mesh_->getMyElements(block,q8_elements);
443 TEUCHOS_ASSERT(q4_elements.size() == q8_elements.size());
444
445 for (size_t i=0; i < q4_elements.size(); ++i) {
446#ifdef PANZER_DEBUG
447 TEUCHOS_ASSERT(q4_mesh.getBulkData()->identifier(q4_elements[i]) ==
448 quad8Mesh_->getBulkData()->identifier(q8_elements[i]));
449#endif
450
451 double* const q4_val = static_cast<double*>(stk::mesh::field_data(*q4_field,q4_elements[i]));
452 const double* const q8_val = static_cast<double*>(stk::mesh::field_data(*q8_field,q8_elements[i]));
453 *q4_val = *q8_val;
454
455 if (print_debug_) {
456 std::cout << "field=" << field_name << ", block=" << block
457 << ", q4e(" << q4_mesh.getBulkData()->identifier(q4_elements[i]) << ") = " << *q4_val
458 << ", q8e(" << quad8Mesh_->getBulkData()->identifier(q8_elements[i]) << ") = " << *q8_val
459 << std::endl;
460 }
461
462 }
463 }
464 }
465 }
466}
467
468} // end panzer_stk
PHX::MDField< ScalarT, panzer::Cell, panzer::IP > result
A field that will be used to build up the result of the integral we're performing.
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral.
virtual void completeMeshConstruction(STK_Interface &mesh, stk::ParallelMachine parallelMach) const
Quad8ToQuad4MeshFactory(const std::string &quad8MeshFileName, stk::ParallelMachine mpi_comm=MPI_COMM_WORLD, const bool print_debug=false)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Derived from ParameterListAcceptor.
Teuchos::RCP< STK_Interface > buildMesh(stk::ParallelMachine parallelMach) const
Build the mesh object.
void buildElements(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
void buildMetaData(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
Teuchos::RCP< panzer_stk::STK_Interface > quad8Mesh_
virtual Teuchos::RCP< STK_Interface > buildUncommitedMesh(stk::ParallelMachine parallelMach) const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
Derived from ParameterListAcceptor.
void initialize(stk::ParallelMachine parallelMach, bool setupIO=true, const bool buildRefinementSupport=false)
static const std::string edgeBlockString
stk::mesh::Part * getEdgeBlock(const std::string &name) const
get the block part
bool isInitialized() const
Has initialize been called on this mesh object?
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
stk::mesh::EntityRank getElementRank() const
void addElement(const Teuchos::RCP< ElementDescriptor > &ed, stk::mesh::Part *block)
void addNode(stk::mesh::EntityId gid, const std::vector< double > &coord)
void addCellField(const std::string &fieldName, const std::string &blockId)
void addNodeset(const std::string &name)
void addSideset(const std::string &name, const CellTopologyData *ctData)
Teuchos::RCP< stk::mesh::MetaData > getMetaData() const
stk::mesh::Field< double > * getCellField(const std::string &fieldName, const std::string &blockId) const
void addEntitiesToEdgeBlock(std::vector< stk::mesh::Entity > entities, stk::mesh::Part *edgeblock)
void addSolutionField(const std::string &fieldName, const std::string &blockId)
void addEntityToSideset(stk::mesh::Entity entity, stk::mesh::Part *sideset)
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::EntityRank getEdgeRank() 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_
Teuchos::RCP< ElementDescriptor > buildElementDescriptor(stk::mesh::EntityId elmtId, std::vector< stk::mesh::EntityId > &nodes)