Sierra Toolkit  Version of the Day
FEMMetaData.cpp
1 #include <set>
2 #include <stk_mesh/fem/FEMMetaData.hpp>
3 
4 #include <Shards_CellTopology.hpp>
5 #include <stk_mesh/base/BulkData.hpp>
6 #include <stk_mesh/base/Bucket.hpp>
7 #include <stk_mesh/base/Entity.hpp>
8 #include <stk_mesh/base/Ghosting.hpp>
9 
10 namespace stk_classic {
11 namespace mesh {
12 namespace fem {
13 
14 namespace {
15 
16 void assign_cell_topology(
17  FEMMetaData::PartCellTopologyVector & part_cell_topology_vector,
18  size_t part_ordinal,
19  const fem::CellTopology cell_topology)
20 {
21  if (part_ordinal >= part_cell_topology_vector.size())
22  part_cell_topology_vector.resize(part_ordinal + 1);
23 
24  part_cell_topology_vector[part_ordinal] = cell_topology;
25 
26  if (!cell_topology.getCellTopologyData())
27  {
28  std::cout << "bad topology in FEMMetaData::assign_cell_topology" << std::endl;
29  }
30 
31  ThrowRequireMsg(cell_topology.getCellTopologyData(), "bad topology in FEMMetaData::assign_cell_topology");
32 }
33 
34 } // namespace
35 
36 FEMMetaData::FEMMetaData()
37 :
38  m_fem_initialized(false),
39  m_spatial_dimension(0),
40  m_side_rank(INVALID_RANK),
41  m_element_rank(INVALID_RANK)
42 {
43  // Attach FEMMetaData as attribute on MetaData to enable "get accessors" to FEMMetaData
44  m_meta_data.declare_attribute_no_delete<FEMMetaData>(this);
45 }
46 
47 FEMMetaData::FEMMetaData(size_t spatial_dimension,
48  const std::vector<std::string>& in_entity_rank_names)
49  :
50  m_fem_initialized(false),
51  m_spatial_dimension(0),
52  m_side_rank(INVALID_RANK),
53  m_element_rank(INVALID_RANK)
54 {
55  // Attach FEMMetaData as attribute on MetaData to enable "get accessors" to FEMMetaData
56  m_meta_data.declare_attribute_no_delete<FEMMetaData>(this);
57 
58  FEM_initialize(spatial_dimension, in_entity_rank_names);
59 }
60 
61 void FEMMetaData::FEM_initialize(size_t spatial_dimension, const std::vector<std::string>& rank_names)
62 {
63  ThrowRequireMsg(!m_fem_initialized,"FEM functionality in FEMMetaData can only be initialized once.");
64  if ( rank_names.empty() ) {
65  m_entity_rank_names = fem::entity_rank_names(spatial_dimension);
66  }
67  else {
68  ThrowRequireMsg(rank_names.size() >= spatial_dimension+1,
69  "Entity rank name vector must name every rank");
70  m_entity_rank_names = rank_names;
71  }
72  internal_set_spatial_dimension_and_ranks(spatial_dimension);
73  m_meta_data.set_entity_rank_names(m_entity_rank_names);
74  m_fem_initialized = true;
75  internal_declare_known_cell_topology_parts();
76 }
77 
78 void FEMMetaData::internal_set_spatial_dimension_and_ranks(size_t spatial_dimension)
79 {
80  ThrowRequireMsg( spatial_dimension != 0, "FEMMetaData::internal_set_spatial_dimension_and_ranks: spatial_dimension == 0!");
81  m_spatial_dimension = spatial_dimension;
82  m_meta_data.m_spatial_dimension = spatial_dimension;
83 
84  // TODO: Decide on correct terminology for the FEM Entity Ranks (consider topological vs spatial names and incompatibilities).
85  // spatial_dimension = 1
86  // node = 0, edge = 0, face = 0, side = 0, element = 1
87  // spatial_dimension = 2
88  // node = 0, edge = 1, face = 1, side = 1, element = 2
89  // spatial_dimension = 3
90  // node = 0, edge = 1, face = 2, side = 2, element = 3
91  // spatial_dimension = 4
92  // node = 0, edge = 1, face = 2, side = 3, element = 4
93  m_side_rank = m_spatial_dimension - 1;
94  m_element_rank = m_spatial_dimension;
95 
96 }
97 
98 void FEMMetaData::internal_declare_known_cell_topology_parts()
99 {
100  // Load up appropriate standard cell topologies.
101  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Node >()), NODE_RANK);
102 
103  if (m_spatial_dimension == 1) {
104 
105  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Particle >()), m_element_rank);
106 
107  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Line<2> >()), m_element_rank); // ???
108  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Line<3> >()), m_element_rank); // ???
109 
110  }
111 
112  else if (m_spatial_dimension == 2) {
113 
114  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Line<2> >()), m_side_rank);
115  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Line<3> >()), m_side_rank);
116 
117  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Particle >()), m_element_rank);
118 
119  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Triangle<3> >()), m_element_rank);
120  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Triangle<6> >()), m_element_rank);
121  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Triangle<4> >()), m_element_rank);
122 
123  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<4> >()), m_element_rank);
124  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<8> >()), m_element_rank);
125  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<9> >()), m_element_rank);
126 
127  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Beam<2> >()), m_element_rank);
128  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Beam<3> >()), m_element_rank);
129 
130  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellLine<2> >()), m_element_rank);
131  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellLine<3> >()), m_element_rank);
132  }
133 
134  else if (m_spatial_dimension == 3) {
135 
136  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Line<2> >()), EDGE_RANK);
137  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Line<3> >()), EDGE_RANK);
138 
139  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Triangle<3> >()), m_side_rank);
140  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Triangle<6> >()), m_side_rank);
141  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Triangle<4> >()), m_side_rank);
142 
143  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<4> >()), m_side_rank);
144  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<8> >()), m_side_rank);
145  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<9> >()), m_side_rank);
146 
147  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Particle >()), m_element_rank);
148 
149  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Beam<2> >()), m_element_rank);
150  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Beam<3> >()), m_element_rank);
151 
152  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Tetrahedron<4> >()), m_element_rank);
153  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Tetrahedron<10> >()), m_element_rank);
154  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Tetrahedron<11> >()), m_element_rank);
155  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Tetrahedron<8> >()), m_element_rank);
156 
157  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Pyramid<5> >()), m_element_rank);
158  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Pyramid<13> >()), m_element_rank);
159  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Pyramid<14> >()), m_element_rank);
160 
161  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Wedge<6> >()), m_element_rank);
162  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Wedge<15> >()), m_element_rank);
163  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Wedge<18> >()), m_element_rank);
164 
165  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Hexahedron<8> >()), m_element_rank);
166  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Hexahedron<20> >()), m_element_rank);
167  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Hexahedron<27> >()), m_element_rank);
168 
169  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellTriangle<3> >()), m_element_rank);
170  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellTriangle<6> >()), m_element_rank);
171 
172  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellQuadrilateral<4> >()), m_element_rank);
173  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellQuadrilateral<8> >()), m_element_rank);
174  register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellQuadrilateral<9> >()), m_element_rank);
175  }
176 }
177 
178 void FEMMetaData::register_cell_topology(const fem::CellTopology cell_topology, EntityRank entity_rank)
179 {
180  ThrowRequireMsg(is_FEM_initialized(),"FEMMetaData::register_cell_topology: FEM_initialize() must be called before this function");
181 
182  CellTopologyPartEntityRankMap::const_iterator it = m_cellTopologyPartEntityRankMap.find(cell_topology);
183 
184  const bool duplicate = it != m_cellTopologyPartEntityRankMap.end();
185  const EntityRank existing_rank = duplicate ? (*it).second.second : 0;
186 
187  ThrowInvalidArgMsgIf(m_spatial_dimension < entity_rank,
188  "entity_rank " << entity_rank << ", " <<
189  "exceeds maximum spatial_dimension = " << m_spatial_dimension );
190 
191  ThrowErrorMsgIf(duplicate && existing_rank != entity_rank,
192  "For args: cell_topolgy " << cell_topology.getName() << " and entity_rank " << entity_rank << ", " <<
193  "previously declared rank = " << existing_rank );
194 
195  if (! duplicate) {
196  std::string part_name = std::string("FEM_ROOT_CELL_TOPOLOGY_PART_") + std::string(cell_topology.getName());
197 
198  ThrowErrorMsgIf(get_part(part_name) != 0, "Cannot register topology with same name as existing part '" << cell_topology.getName() << "'" );
199 
200  Part &part = declare_internal_part(part_name, entity_rank);
201  m_cellTopologyPartEntityRankMap[cell_topology] = CellTopologyPartEntityRankMap::mapped_type(&part, entity_rank);
202 
203  assign_cell_topology(m_partCellTopologyVector, part.mesh_meta_data_ordinal(), cell_topology);
204  }
205  //check_topo_db();
206 }
207 
208 
209 fem::CellTopology
211  const std::string & topology_name) const
212 {
213  std::string part_name = convert_to_internal_name(std::string("FEM_ROOT_CELL_TOPOLOGY_PART_") + topology_name);
214 
215  Part *part = get_part(part_name);
216  if (part)
217  return get_cell_topology(*part);
218  else
219  return fem::CellTopology();
220 }
221 
222 
223 Part &FEMMetaData::get_cell_topology_root_part(const fem::CellTopology cell_topology) const
224 {
225  ThrowRequireMsg(is_FEM_initialized(),"FEMMetaData::get_cell_topology_root_part: FEM_initialize() must be called before this function");
226  CellTopologyPartEntityRankMap::const_iterator it = m_cellTopologyPartEntityRankMap.find(cell_topology);
227  ThrowErrorMsgIf(it == m_cellTopologyPartEntityRankMap.end(),
228  "Cell topology " << cell_topology.getName() <<
229  " has not been registered");
230 
231  return *(*it).second.first;
232 }
233 
239 fem::CellTopology FEMMetaData::get_cell_topology( const Part & part) const
240 {
241  ThrowRequireMsg(is_FEM_initialized(),"FEMMetaData::get_cell_topology: FEM_initialize() must be called before this function");
242  fem::CellTopology cell_topology;
243 
244  PartOrdinal part_ordinal = part.mesh_meta_data_ordinal();
245  if (part_ordinal < m_partCellTopologyVector.size())
246  {
247  cell_topology = m_partCellTopologyVector[part_ordinal];
248  }
249 
250  return cell_topology;
251 }
252 
253 #if 0
254  void FEMMetaData::check_topo_db()
255  {
256  std::cout << "FEMMetaData::check_topo_db... m_partCellTopologyVector.size() = " << m_partCellTopologyVector.size() << std::endl;
257 
258  fem::CellTopology cell_topology;
259 
260  for (unsigned i = 0; i < m_partCellTopologyVector.size(); i++)
261  {
262  cell_topology = m_partCellTopologyVector[i];
263  if (!cell_topology.getCellTopologyData())
264  {
265  std::cout << "bad topology in FEMMetaData::check_topo_db" << std::endl;
266  }
267  ThrowRequireMsg(cell_topology.getCellTopologyData(), "bad topology in FEMMetaData::check_topo_db");
268 
269  }
270  std::cout << "FEMMetaData::check_topo_db...done" << std::endl;
271 
272  }
273 #endif
274 
275 namespace {
276 
277 bool root_part_in_subset(stk_classic::mesh::Part & part)
278 {
279  if (is_cell_topology_root_part(part)) {
280  return true;
281  }
282  const PartVector & subsets = part.subsets();
283  for (PartVector::const_iterator it=subsets.begin() ; it != subsets.end() ; ++it) {
284  if (is_cell_topology_root_part( **it )) {
285  return true;
286  }
287  }
288  return false;
289 }
290 
291 void find_cell_topologies_in_part_and_subsets_of_same_rank(const Part & part, EntityRank rank, std::set<fem::CellTopology> & topologies_found)
292 {
293  fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(part);
294  fem::CellTopology top = fem_meta.get_cell_topology(part);
295  if ((top.isValid() && (part.primary_entity_rank() == rank))) {
296  topologies_found.insert(top);
297  }
298  const PartVector & subsets = part.subsets();
299  for (PartVector::const_iterator it=subsets.begin() ; it != subsets.end() ; ++it) {
300  top = fem_meta.get_cell_topology(**it);
301  if (top.isValid() && ( (**it).primary_entity_rank() == rank) ) {
302  topologies_found.insert(top);
303  }
304  }
305 }
306 
307 } // namespace
308 
309 void FEMMetaData::declare_part_subset( Part & superset , Part & subset )
310 {
311  ThrowRequireMsg(is_FEM_initialized(),"FEMMetaData::declare_part_subset: FEM_initialize() must be called before this function");
312  fem::CellTopology superset_top = get_cell_topology(superset);
313 
314  const bool no_superset_topology = !superset_top.isValid();
315  if ( no_superset_topology ) {
316  m_meta_data.declare_part_subset(superset,subset);
317  return;
318  }
319  // Check for cell topology root parts in subset or subset's subsets
320  const bool subset_has_root_part = root_part_in_subset(subset);
321  ThrowErrorMsgIf( subset_has_root_part, "FEMMetaData::declare_part_subset: Error, root cell topology part found in subset or below." );
322 
323  std::set<fem::CellTopology> cell_topologies;
324  find_cell_topologies_in_part_and_subsets_of_same_rank(subset,superset.primary_entity_rank(),cell_topologies);
325 
326  ThrowErrorMsgIf( cell_topologies.size() > 1,
327  "FEMMetaData::declare_part_subset: Error, multiple cell topologies of rank "
328  << superset.primary_entity_rank()
329  << " defined below subset"
330  );
331  const bool non_matching_cell_topology = ((cell_topologies.size() == 1) && (*cell_topologies.begin() != superset_top));
332  ThrowErrorMsgIf( non_matching_cell_topology,
333  "FEMMetaData::declare_part_subset: Error, superset topology = "
334  << superset_top.getName() << " does not match the topology = "
335  << cell_topologies.begin()->getName()
336  << " coming from the subset part"
337  );
338  // Everything is Okay!
339  m_meta_data.declare_part_subset(superset,subset);
340  // Update PartCellTopologyVector for "subset" and same-rank subsets, ad nauseum
341  if (subset.primary_entity_rank() == superset.primary_entity_rank()) {
342  assign_cell_topology(m_partCellTopologyVector, subset.mesh_meta_data_ordinal(), superset_top);
343  const PartVector & subset_parts = subset.subsets();
344  for (PartVector::const_iterator it=subset_parts.begin() ; it != subset_parts.end() ; ++it) {
345  const Part & it_part = **it;
346  if (it_part.primary_entity_rank() == superset.primary_entity_rank()) {
347  assign_cell_topology(m_partCellTopologyVector, it_part.mesh_meta_data_ordinal(), superset_top);
348  }
349  }
350  }
351 }
352 
354  const fem::CellTopology cell_topology) const
355 {
356  CellTopologyPartEntityRankMap::const_iterator it = m_cellTopologyPartEntityRankMap.find(cell_topology);
357  if (it == m_cellTopologyPartEntityRankMap.end())
358  return INVALID_RANK;
359  else
360  return (*it).second.second;
361 }
362 
363 //--------------------------------------------------------------------------------
364 // Free Functions
365 //--------------------------------------------------------------------------------
366 
367 bool is_cell_topology_root_part(const Part & part) {
368  fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(part);
369  fem::CellTopology top = fem_meta.get_cell_topology(part);
370  if (top.isValid()) {
371  const Part & root_part = fem_meta.get_cell_topology_root_part(top);
372  return (root_part == part);
373  }
374  return false;
375 }
376 
382 void set_cell_topology(
383  Part & part,
384  fem::CellTopology cell_topology)
385 {
386  FEMMetaData& fem_meta = FEMMetaData::get(part);
387 
388  ThrowRequireMsg(fem_meta.is_FEM_initialized(),"set_cell_topology: FEM_initialize() must be called before this function");
389 
390  Part &root_part = fem_meta.get_cell_topology_root_part(cell_topology);
391  fem_meta.declare_part_subset(root_part, part);
392 }
393 
394 std::vector<std::string>
395 entity_rank_names( size_t spatial_dimension )
396 {
397  ThrowInvalidArgMsgIf( spatial_dimension < 1 || 3 < spatial_dimension,
398  "Invalid spatial dimension = " << spatial_dimension );
399 
400  std::vector< std::string > names ;
401 
402  names.reserve( spatial_dimension + 1 );
403 
404  names.push_back( std::string( "NODE" ) );
405 
406  if ( 1 < spatial_dimension ) { names.push_back( std::string("EDGE") ); }
407  if ( 2 < spatial_dimension ) { names.push_back( std::string("FACE") ); }
408 
409  names.push_back( std::string("ELEMENT") );
410 
411  return names ;
412 }
413 
414 
415 CellTopology
416 get_cell_topology(
417  const Bucket & bucket)
418 {
419  const BulkData & bulk_data = BulkData::get(bucket);
420  const MetaData & meta_data = MetaData::get(bulk_data);
421  const PartVector & all_parts = meta_data.get_parts();
422 
423  FEMMetaData &fem = FEMMetaData::get(meta_data);
424 
425  CellTopology cell_topology;
426 
427  const std::pair< const unsigned *, const unsigned * > supersets = bucket.superset_part_ordinals();
428 
429  if (supersets.first != supersets.second) {
430  const Part *first_found_part = 0;
431 
432  for ( const unsigned * it = supersets.first ; it != supersets.second ; ++it ) {
433 
434  const Part & part = * all_parts[*it] ;
435 
436  if ( part.primary_entity_rank() == bucket.entity_rank() ) {
437 
438  CellTopology top = fem.get_cell_topology( part );
439 
440  if ( ! cell_topology.getCellTopologyData() ) {
441  cell_topology = top ;
442 
443  if (!first_found_part)
444  first_found_part = &part;
445  }
446  else {
447  ThrowErrorMsgIf( top.getCellTopologyData() && top != cell_topology,
448  "Cell topology is ambiguously defined. It is defined as " << cell_topology.getName() <<
449  " on part " << first_found_part->name() << " and as " << top.getName() << " on its superset part " << part.name() );
450  }
451  }
452  }
453  }
454 
455  return cell_topology ;
456 }
457 
458 } // namespace fem
459 } // namespace mesh
460 } // namespace stk_classic
FEMMetaData is a class that implements a Finite Element Method skin on top of the Sierra Tool Kit Met...
Definition: FEMMetaData.hpp:54
void declare_part_subset(Part &superset, Part &subset)
Declare a superset-subset relationship between parts.
Definition: MetaData.cpp:240
void register_cell_topology(const fem::CellTopology cell_topology, EntityRank in_entity_rank)
This function is used to register new cell topologies and their associated ranks with FEMMetaData...
unsigned primary_entity_rank() const
The primary entity type for this part.
Definition: Part.hpp:64
An application-defined subset of a problem domain.
Definition: Part.hpp:49
bool is_FEM_initialized() const
This function returns whether this class has been initialized or not.
Part & get_cell_topology_root_part(const fem::CellTopology cell_topology) const
Return the root cell topology part associated with the given cell topology. This Part is created in r...
void declare_part_subset(Part &superset, Part &subset)
Declare a superset-subset relationship between parts Note: Cell Topologies are induced through part s...
unsigned mesh_meta_data_ordinal() const
Internally generated ordinal of this part that is unique within the owning meta data manager...
Definition: Part.hpp:72
void set_entity_rank_names(const std::vector< std::string > &entity_rank_names)
entity-rank names
Definition: MetaData.cpp:155
Part * get_part(const std::string &p_name, const char *required_by=NULL) const
Get an existing part by its application-defined text name.
const PartVector & subsets() const
Parts that are subsets of this part.
Definition: Part.hpp:78
Sierra Toolkit.
std::vector< fem::CellTopology > PartCellTopologyVector
PartCellTopologyVector is a fast-lookup vector of size equal to the number of parts.
Definition: FEMMetaData.hpp:60
size_t spatial_dimension() const
Returns the spatial dimension that was passed in through FEM_initialize.
EntityRank get_entity_rank(const fem::CellTopology cell_topology) const
Return the EntityRank that is associated with the given cell topology. In several cases...
std::vector< Part *> PartVector
Collections of parts are frequently maintained as a vector of Part pointers.
Definition: Types.hpp:31
static FEMMetaData & get(const MetaData &meta)
Getter for FEMMetaData off of a MetaData object.
EntityRank entity_rank(const std::string &name) const
Return the rank for the given name that was provided in set_entity_rank_names.
EntityRank entity_rank(const EntityKey &key)
Given an entity key, return an entity type (rank).
fem::CellTopology get_cell_topology(const Part &part) const
Return the cell topology associated with the given part. The cell topology is set on a part through p...
void FEM_initialize(size_t spatial_dimension, const std::vector< std::string > &in_entity_rank_names=std::vector< std::string >())
Initialize the spatial dimension and an optional list of entity rank names associated with each rank...
Definition: FEMMetaData.cpp:61