Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_STK_QuadraticToLinearMeshFactory.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
51using Teuchos::RCP;
52using Teuchos::rcp;
53
54namespace panzer_stk {
55
57 stk::ParallelMachine mpi_comm,
58 const bool print_debug)
59 : createEdgeBlocks_(false),
60 print_debug_(print_debug)
61{
62 panzer_stk::STK_ExodusReaderFactory factory;
63 RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList);
64 pl->set("File Name",quadMeshFileName);
65 factory.setParameterList(pl);
66 quadMesh_ = factory.buildMesh(mpi_comm);
67
68 // TODO for currently supported input topologies, this should always be true
69 // but may need to be generalized in the future
71
72 // check the conversion is supported
73 // and get output topology
74 this->getOutputTopology();
75}
76
77QuadraticToLinearMeshFactory::QuadraticToLinearMeshFactory(const Teuchos::RCP<panzer_stk::STK_Interface>& quadMesh,
78 const bool print_debug)
79 : quadMesh_(quadMesh),
80 createEdgeBlocks_(false),
81 print_debug_(print_debug)
82{
83 // TODO for currently supported input topologies, this should always be true
84 // but may need to be generalized in the future
86
87 // check the conversion is supported
88 // and get output topology
89 this->getOutputTopology();
90}
91
93Teuchos::RCP<STK_Interface> QuadraticToLinearMeshFactory::buildMesh(stk::ParallelMachine parallelMach) const
94{
95 PANZER_FUNC_TIME_MONITOR("panzer::QuadraticToLinearMeshFactory::buildMesh()");
96
97 // Make sure the Comms match between the meshes
98 {
99 int result = MPI_UNEQUAL;
100 MPI_Comm_compare(parallelMach, quadMesh_->getBulkData()->parallel(), &result);
101 TEUCHOS_ASSERT(result != MPI_UNEQUAL);
102 }
103
104 // build all meta data
105 RCP<STK_Interface> mesh = this->buildUncommitedMesh(parallelMach);
106
107 // commit meta data
108 mesh->initialize(parallelMach);
109
110 // build bulk data
111 this->completeMeshConstruction(*mesh,parallelMach);
112
113 return mesh;
114}
115
117{
118 bool errFlag = false;
119
120 std::vector<std::string> eblock_names;
121 quadMesh_->getElementBlockNames(eblock_names);
122
123 // check that we have a supported topology
124 auto inputTopo = quadMesh_->getCellTopology(eblock_names[0]);
125 if (std::find(supportedInputTopos_.begin(),
126 supportedInputTopos_.end(),*inputTopo) == supportedInputTopos_.end()) errFlag = true;
127 TEUCHOS_TEST_FOR_EXCEPTION(errFlag,std::logic_error,
128 "ERROR :: Input topology " << inputTopo->getName() << " currently unsupported by QuadraticToLinearMeshFactory!");
129
130 // check that the topology is the same over blocks
131 // not sure this is 100% foolproof
132 for (auto & eblock : eblock_names) {
133 auto cellTopo = quadMesh_->getCellTopology(eblock);
134 if (*cellTopo != *inputTopo) errFlag = true;
135 }
136
137 TEUCHOS_TEST_FOR_EXCEPTION(errFlag, std::logic_error,
138 "ERROR :: The mesh has different topologies on different blocks!");
139
140 outputTopoData_ = outputTopoMap_[inputTopo->getName()];
141
142 nDim_ = outputTopoData_->dimension;
143 nNodes_ = outputTopoData_->node_count;
144
145 return;
146}
147
148Teuchos::RCP<STK_Interface> QuadraticToLinearMeshFactory::buildUncommitedMesh(stk::ParallelMachine parallelMach) const
149{
150 PANZER_FUNC_TIME_MONITOR("panzer::QuadraticToLinearMeshFactory::buildUncomittedMesh()");
151
152 RCP<STK_Interface> mesh = rcp(new STK_Interface(nDim_));
153
154 machRank_ = stk::parallel_machine_rank(parallelMach);
155 machSize_ = stk::parallel_machine_size(parallelMach);
156
157 // build meta information: blocks and side set setups
158 this->buildMetaData(parallelMach,*mesh);
159
160 mesh->addPeriodicBCs(periodicBCVec_);
161 mesh->setBoundingBoxSearchFlag(useBBoxSearch_);
162
163 return mesh;
164}
165
166void QuadraticToLinearMeshFactory::completeMeshConstruction(STK_Interface & mesh,stk::ParallelMachine parallelMach) const
167{
168 PANZER_FUNC_TIME_MONITOR("panzer::QuadraticToLinearMeshFactory::completeMeshConstruction()");
169
170 if(not mesh.isInitialized())
171 mesh.initialize(parallelMach);
172
173 // add node and element information
174 this->buildElements(parallelMach,mesh);
175
176 // finish up the edges
177#ifndef ENABLE_UNIFORM
178 mesh.buildSubcells();
179#endif
182 mesh.buildLocalEdgeIDs();
183 }
184
185 // now that edges are built, sidesets can be added
186#ifndef ENABLE_UNIFORM
187 this->addSideSets(mesh);
188#endif
189
190 // add nodesets
191 this->addNodeSets(mesh);
192
193 // TODO this functionality may be untested
195 this->addEdgeBlocks(mesh);
196 }
197
198 // Copy field data
199 this->copyCellFieldData(mesh);
200
201 // calls Stk_MeshFactory::rebalance
202 // this->rebalance(mesh);
203}
204
206void QuadraticToLinearMeshFactory::setParameterList(const Teuchos::RCP<Teuchos::ParameterList> & paramList)
207{
208 paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
209
210 setMyParamList(paramList);
211
212 // offsetGIDs_ = (paramList->get<std::string>("Offset mesh GIDs above 32-bit int limit") == "ON") ? true : false;
213
214 createEdgeBlocks_ = paramList->get<bool>("Create Edge Blocks");
215
216 // read in periodic boundary conditions
217 parsePeriodicBCList(Teuchos::rcpFromRef(paramList->sublist("Periodic BCs")),periodicBCVec_,useBBoxSearch_);
218}
219
221Teuchos::RCP<const Teuchos::ParameterList> QuadraticToLinearMeshFactory::getValidParameters() const
222{
223 static RCP<Teuchos::ParameterList> defaultParams;
224
225 // fill with default values
226 if(defaultParams == Teuchos::null) {
227 defaultParams = rcp(new Teuchos::ParameterList);
228
229 Teuchos::setStringToIntegralParameter<int>(
230 "Offset mesh GIDs above 32-bit int limit",
231 "OFF",
232 "If 64-bit GIDs are supported, the mesh element and node global indices will start at a value greater than 32-bit limit.",
233 Teuchos::tuple<std::string>("OFF", "ON"),
234 defaultParams.get());
235
236 // default to false for backward compatibility
237 defaultParams->set<bool>("Create Edge Blocks",false,"Create edge blocks in the mesh");
238
239 Teuchos::ParameterList & bcs = defaultParams->sublist("Periodic BCs");
240 bcs.set<int>("Count",0); // no default periodic boundary conditions
241 }
242
243 return defaultParams;
244}
245
246void QuadraticToLinearMeshFactory::buildMetaData(stk::ParallelMachine /* parallelMach */, STK_Interface & mesh) const
247{
248 if (print_debug_) {
249 std::cout << "\n\n**** DEBUG: begin printing source quad mesh exodus file metadata ****\n";
250 quadMesh_->getMetaData()->dump_all_meta_info(std::cout);
251 std::cout << "\n\n**** DEBUG: end printing source quad mesh exodus file metadata ****\n";
252 }
253
254 // Required topologies
255 auto ctd = outputTopoData_;
256 // will only work for 2D and 3D meshes
257 const CellTopologyData * side_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(nDim_-1,0);
258
259 // Add in element blocks
260 std::vector<std::string> element_block_names;
261 quadMesh_->getElementBlockNames(element_block_names);
262 {
263 for (const auto& n : element_block_names)
264 mesh.addElementBlock(n,ctd);
265 }
266
267 // Add in sidesets
268 {
269 std::vector<std::string> sideset_names;
270 quadMesh_->getSidesetNames(sideset_names);
271 for (const auto& n : sideset_names)
272 mesh.addSideset(n,side_ctd);
273 }
274
275 // Add in nodesets
276 {
277 std::vector<std::string> nodeset_names;
278 quadMesh_->getNodesetNames(nodeset_names);
279 for (const auto& n : nodeset_names)
280 mesh.addNodeset(n);
281 }
282
284 const CellTopologyData * edge_ctd = shards::CellTopology(ctd).getBaseCellTopologyData(1,0);
285 std::vector<std::string> element_block_names;
286 quadMesh_->getElementBlockNames(element_block_names);
287 for (const auto& block_name : element_block_names)
288 mesh.addEdgeBlock(block_name,edgeBlockName_,edge_ctd);
289 }
290
291 // Add in nodal fields
292 {
293 const auto& fields = quadMesh_->getMetaData()->get_fields(mesh.getNodeRank());
294 for (const auto& f : fields) {
295 if (print_debug_)
296 std::cout << "Field=" << f->name() << ", rank=" << f->entity_rank() << std::endl;
297
298 // Cull out the coordinate fields. That is a default field in
299 // stk that is automatically created by stk.
300 if (f->name() != "coordinates") {
301 for (const auto& n : element_block_names)
302 mesh.addSolutionField(f->name(),n);
303 }
304 }
305 }
306
307 // Add in element fields
308 {
309 const auto& fields = quadMesh_->getMetaData()->get_fields(mesh.getElementRank());
310 for (const auto& f : fields) {
311 if (print_debug_)
312 std::cout << "Add Cell Field=" << f->name() << ", rank=" << f->entity_rank() << std::endl;
313
314 for (const auto& n : element_block_names)
315 mesh.addCellField(f->name(),n);
316 }
317 }
318
319 // NOTE: skipping edge and face fields for now. Can't error out since sidesets count as edge fields.
320 // TEUCHOS_TEST_FOR_EXCEPTION(quadMesh_->getMetaData()->get_fields(mesh.getEdgeRank()).size() != 0,std::runtime_error,
321 // "ERROR: the Quad8 mesh contains edge fields\""
322 // << quadMesh_->getMetaData()->get_fields(mesh.getEdgeRank())[0]->name()
323 // << "\". Edge fields are not supported in Quad8ToQuad4!");
324
325 if (print_debug_) {
326 std::cout << "\n\n**** DEBUG: begin printing source linear mesh exodus file metadata ****\n";
327 mesh.getMetaData()->dump_all_meta_info(std::cout);
328 std::cout << "\n\n**** DEBUG: end printing source linear mesh exodus file metadata ****\n";
329 }
330}
331
332void QuadraticToLinearMeshFactory::buildElements(stk::ParallelMachine parallelMach,STK_Interface & mesh) const
333{
334 mesh.beginModification();
335
336 auto metadata = mesh.getMetaData();
337 auto bulkdata = mesh.getBulkData();
338
339 // Loop over element blocks
340 std::vector<std::string> block_names;
341 quadMesh_->getElementBlockNames(block_names);
342 for (const auto& block_name : block_names) {
343
344 // Get the elements on this process
345 std::vector<stk::mesh::Entity> elements;
346 quadMesh_->getMyElements(block_name,elements);
347
348 if (print_debug_) {
349 std::cout << "*************************************************" << std::endl;
350 std::cout << "block_name=" << block_name << ", num_my_elements=" << elements.size() << std::endl;
351 std::cout << "*************************************************" << std::endl;
352 }
353
354 for (const auto& element : elements) {
355
356 const auto element_gid = quadMesh_->getBulkData()->identifier(element);
357
358 if (print_debug_) {
359 std::cout << "rank=" << machRank_
360 << ", block=" << block_name
361 << ", element entity_id=" << element
362 << ", gid=" << element_gid << std::endl;
363 }
364
365 // Register nodes with the mesh
366 std::vector<stk::mesh::EntityId> nodes(nNodes_);
367 for (size_t i=0; i < nNodes_; ++i) {
368 // NOTE: this assumes that the linear topology is nested in
369 // the quadratic topology as an extended topology - i.e. the first n
370 // nodes of the quadratic topology are the corresponding linear
371 // nodes. Shards topologies enforce this through the concept
372 // of extended topologies.
373 stk::mesh::Entity node_entity = quadMesh_->findConnectivityById(element,mesh.getNodeRank(),i);
374 TEUCHOS_ASSERT(node_entity.is_local_offset_valid());
375 const auto node_gid = quadMesh_->getBulkData()->identifier(node_entity);
376 const double* node_coords = quadMesh_->getNodeCoordinates(node_entity);
377 std::vector<double> coords_vec(nDim_);
378 for (size_t j=0; j < nDim_; ++j) coords_vec[j] = node_coords[j];
379 mesh.addNode(node_gid,coords_vec);
380 nodes[i]=node_gid;
381 if (print_debug_) {
382 if (nDim_==2) {
383 std::cout << "elem gid=" << quadMesh_->getBulkData()->identifier(element)
384 << ", node_gid=" << node_gid << ", ("
385 << coords_vec[0] << "," << coords_vec[1] << ")\n";
386 } else {
387 std::cout << "elem gid=" << quadMesh_->getBulkData()->identifier(element)
388 << ", node_gid=" << node_gid << ", ("
389 << coords_vec[0] << "," << coords_vec[1] << "," << coords_vec[2] << ")\n";
390 }
391 }
392 }
393
394 // Register element with the element block
395 auto element_descriptor = panzer_stk::buildElementDescriptor(element_gid,nodes);
396 auto element_block_part = mesh.getMetaData()->get_part(block_name);
397 mesh.addElement(element_descriptor,element_block_part);
398 }
399 }
400 mesh.endModification();
401}
402
404{
405 mesh.beginModification();
406
407 // Loop over sidesets
408 std::vector<std::string> sideset_names;
409 quadMesh_->getSidesetNames(sideset_names);
410 for (const auto& sideset_name : sideset_names) {
411
412 stk::mesh::Part* sideset_part = mesh.getSideset(sideset_name);
413
414 std::vector<stk::mesh::Entity> q_sides;
415 quadMesh_->getMySides(sideset_name,q_sides);
416
417 // Loop over edges
418 for (const auto q_ent : q_sides) {
419 // The edge numbering scheme uses the element/node gids, so it
420 // should be consistent between the quadratic and linear meshes
421 // since we used the same gids. We use this fact to populate
422 // the quadratic sidesets.
423 stk::mesh::EntityId ent_gid = quadMesh_->getBulkData()->identifier(q_ent);
424 stk::mesh::Entity lin_ent = mesh.getBulkData()->get_entity(mesh.getSideRank(),ent_gid);
425 mesh.addEntityToSideset(lin_ent,sideset_part);
426 }
427 }
428
429 mesh.endModification();
430}
431
433{}
434
436{
437 mesh.beginModification();
438
439 Teuchos::RCP<stk::mesh::BulkData> bulkData = mesh.getBulkData();
440 Teuchos::RCP<stk::mesh::MetaData> metaData = mesh.getMetaData();
441
442 stk::mesh::Part * edge_block = mesh.getEdgeBlock(edgeBlockName_);
443
444 stk::mesh::Selector owned_block = metaData->locally_owned_part();
445
446 std::vector<stk::mesh::Entity> edges;
447 bulkData->get_entities(mesh.getEdgeRank(), owned_block, edges);
448 mesh.addEntitiesToEdgeBlock(edges, edge_block);
449
450 mesh.endModification();
451}
452
454{
455 // Vector of pointers to field data
456 const auto& fields = lin_mesh.getMetaData()->get_fields(lin_mesh.getElementRank());
457
458 // loop over fields and add the data to the new mesh.
459 for (const auto& field : fields) {
460
461 if (print_debug_)
462 std::cout << "Adding field values for \"" << *field << "\" to the linear mesh!\n";
463
464 auto field_name = field->name();
465
466 // Divide into scalar and vector fields, ignoring any other types
467 // for now.
468 if (field->type_is<double>() &&
469 field_name != "coordinates" &&
470 field_name != "PROC_ID" &&
471 field_name != "LOAD_BAL") {
472
473 // Loop over element blocks
474 std::vector<std::string> block_names;
475 quadMesh_->getElementBlockNames(block_names);
476 for (const auto& block : block_names) {
477
478 auto* lin_field = lin_mesh.getCellField(field_name,block);
479 TEUCHOS_ASSERT(lin_field != nullptr);
480 // The q mesh doesn't have the field names set up, so a query
481 // fails. Go to stk directly in this case.
482 auto* q_field = quadMesh_->getMetaData()->get_field(quadMesh_->getElementRank(),field_name);
483#ifdef PANZER_DEBUG
484 TEUCHOS_ASSERT(q_field != nullptr);
485#endif
486
487 // Get the elements for this block.
488 std::vector<stk::mesh::Entity> lin_elements;
489 lin_mesh.getMyElements(block,lin_elements);
490 std::vector<stk::mesh::Entity> q_elements;
491 quadMesh_->getMyElements(block,q_elements);
492 TEUCHOS_ASSERT(lin_elements.size() == q_elements.size());
493
494 for (size_t i=0; i < lin_elements.size(); ++i) {
495#ifdef PANZER_DEBUG
496 TEUCHOS_ASSERT(lin_mesh.getBulkData()->identifier(lin_elements[i]) ==
497 quadMesh_->getBulkData()->identifier(q_elements[i]));
498#endif
499
500 double* const lin_val = static_cast<double*>(stk::mesh::field_data(*lin_field,lin_elements[i]));
501 const double* const q_val = static_cast<double*>(stk::mesh::field_data(*q_field,q_elements[i]));
502 *lin_val = *q_val;
503
504 if (print_debug_) {
505 std::cout << "field=" << field_name << ", block=" << block
506 << ", lin_e(" << lin_mesh.getBulkData()->identifier(lin_elements[i]) << ") = " << *lin_val
507 << ", q_e(" << quadMesh_->getBulkData()->identifier(q_elements[i]) << ") = " << *q_val
508 << std::endl;
509 }
510
511 }
512 }
513 }
514 }
515}
516
517} // 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
void buildMetaData(stk::ParallelMachine parallelMach, STK_Interface &mesh) const
std::map< const std::string, const CellTopologyData * > outputTopoMap_
QuadraticToLinearMeshFactory(const std::string &quadMeshFileName, stk::ParallelMachine mpi_comm=MPI_COMM_WORLD, const bool print_debug=false)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Derived from ParameterListAcceptor.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
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
virtual Teuchos::RCP< STK_Interface > buildUncommitedMesh(stk::ParallelMachine parallelMach) const
std::vector< shards::CellTopology > supportedInputTopos_
Nodes in one element of the linear basis.
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)
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::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)