Sierra Toolkit  Version of the Day
UnitTestBoundaryAnalysis.cpp
1 /*------------------------------------------------------------------------*/
2 /* Copyright 2010, 2011 Sandia Corporation. */
3 /* Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive */
4 /* license for use of this work by or on behalf of the U.S. Government. */
5 /* Export of this program may require a license from the */
6 /* United States Government. */
7 /*------------------------------------------------------------------------*/
8 
9 #include <stk_util/unit_test_support/stk_utest_macros.hpp>
10 #include <Shards_BasicTopologies.hpp>
11 
12 #include <stk_util/parallel/Parallel.hpp>
13 
14 #include <stk_mesh/base/MetaData.hpp>
15 #include <stk_mesh/base/BulkData.hpp>
16 #include <stk_mesh/base/Entity.hpp>
17 #include <stk_mesh/base/GetEntities.hpp>
18 #include <stk_mesh/base/Selector.hpp>
19 #include <stk_mesh/base/GetBuckets.hpp>
20 
21 #include <stk_mesh/fem/BoundaryAnalysis.hpp>
22 #include <stk_mesh/fem/FEMHelpers.hpp>
23 
24 #include <stk_mesh/fixtures/GridFixture.hpp>
25 
26 #include <iomanip>
27 #include <algorithm>
28 
29 static const size_t NODE_RANK = stk_classic::mesh::fem::FEMMetaData::NODE_RANK;
30 
31 class UnitTestStkMeshBoundaryAnalysis {
32 public:
33  UnitTestStkMeshBoundaryAnalysis(stk_classic::ParallelMachine pm) : m_comm(pm), m_num_procs(0), m_rank(0)
34  {
35  m_num_procs = stk_classic::parallel_machine_size( m_comm );
36  m_rank = stk_classic::parallel_machine_rank( m_comm );
37  }
38 
39  void test_boundary_analysis();
40  void test_boundary_analysis_null_topology();
41 
43  int m_num_procs;
44  int m_rank;
45 };
46 
47 namespace {
48 
49 STKUNIT_UNIT_TEST( UnitTestStkMeshBoundaryAnalysis , testUnit )
50 {
51  UnitTestStkMeshBoundaryAnalysis unit(MPI_COMM_WORLD);
52  unit.test_boundary_analysis();
53 }
54 
55 STKUNIT_UNIT_TEST( UnitTestStkMeshBoundaryAnalysis , testNullTopology )
56 {
57  UnitTestStkMeshBoundaryAnalysis unit(MPI_COMM_WORLD);
58  unit.test_boundary_analysis_null_topology();
59 }
60 
61 } //end namespace
62 
63 void UnitTestStkMeshBoundaryAnalysis::test_boundary_analysis()
64 {
65  // Test the boundary_analysis algorithm in stk_mesh/fem/BoundaryAnalysis.hpp
66  // with a boundary that is both shelled and non-shelled.
67  //
68  // We will be testing the algorithm on the following 2D mesh:
69  //
70  // 17---18---19---20---21
71  // | 1 | 2 | 3 || 4 |
72  // 22---23---24---25---26
73  // | 5 | 6 | 7 || 8 |
74  // 27---28---29---30---31
75  // | 9 | 10 | 11 ||12 |
76  // 32---33---34---35---36
77  // | 13 | 14 | 15 ||16 |
78  // 37---38---39---40---41
79  //
80  // Note the shells along nodes 20-40.
81  //
82  // We will be computing the boundary of the closure of
83  // elements: 6, 7, 10, 11, 14, 15
84 
85  // This test will only work for np=1
86  if (m_num_procs > 1) {
87  return;
88  }
89 
90  // set up grid_mesh
91  stk_classic::mesh::fixtures::GridFixture grid_mesh(MPI_COMM_WORLD);
92 
93  stk_classic::mesh::fem::FEMMetaData& fem_meta = grid_mesh.fem_meta();
94  stk_classic::mesh::BulkData& bulk_data = grid_mesh.bulk_data();
95 
96  const stk_classic::mesh::EntityRank element_rank = fem_meta.element_rank();
97 
98  // make shell part
99  stk_classic::mesh::fem::CellTopology line_top(shards::getCellTopologyData<shards::ShellLine<2> >());
100  stk_classic::mesh::Part& shell_part = fem_meta.declare_part("shell_part", line_top);
101 
102  fem_meta.commit();
103 
104  bulk_data.modification_begin();
105  grid_mesh.generate_grid();
106 
107  // Add some shells
108  const unsigned num_shells = 4;
109 
110  // get a count of entities that have already been created
111  std::vector<unsigned> count;
112  stk_classic::mesh::Selector locally_owned(fem_meta.locally_owned_part());
113  stk_classic::mesh::count_entities(locally_owned, bulk_data, count);
114  const unsigned num_entities = count[NODE_RANK] + count[element_rank];
115 
116  // Declare the shell entities, placing them in the shell part
117  std::vector<stk_classic::mesh::Entity*> shells;
118  stk_classic::mesh::PartVector shell_parts;
119  shell_parts.push_back(&shell_part);
120  for (unsigned i = 1; i <= num_shells; ++i) {
121  stk_classic::mesh::Entity& new_shell = bulk_data.declare_entity(element_rank,
122  num_entities + i,
123  shell_parts);
124  shells.push_back(&new_shell);
125  }
126 
127  // declare shell relationships
128  unsigned node_list[5] = {20, 25, 30, 35, 40};
129  for (unsigned i = 0; i < num_shells; ++i) {
130  stk_classic::mesh::Entity& shell = *(shells[i]);
131  stk_classic::mesh::Entity& node1 = *(bulk_data.get_entity(NODE_RANK, node_list[i]));
132  stk_classic::mesh::Entity& node2 = *(bulk_data.get_entity(NODE_RANK, node_list[i+1]));
133  bulk_data.declare_relation(shell, node1, 0);
134  bulk_data.declare_relation(shell, node2, 1);
135  }
136 
137  bulk_data.modification_end();
138 
139  // create the closure we want to analyze
140  std::vector<stk_classic::mesh::Entity*> closure;
141  unsigned num_elems_in_closure = 6;
142  stk_classic::mesh::EntityId ids_of_entities_in_closure[] =
143  {6, 7, 10, 11, 14, 15, 23, 24, 25, 28, 29, 30, 33, 34, 35, 38, 39, 40};
144  for (unsigned i = 0;
145  i < sizeof(ids_of_entities_in_closure)/sizeof(stk_classic::mesh::EntityId);
146  ++i) {
147  stk_classic::mesh::EntityRank rank_of_entity;
148  if (i < num_elems_in_closure) {
149  rank_of_entity = element_rank;
150  }
151  else {
152  rank_of_entity = NODE_RANK;
153  }
154  stk_classic::mesh::Entity* closure_entity =
155  bulk_data.get_entity(rank_of_entity, ids_of_entities_in_closure[i]);
156  closure.push_back(closure_entity);
157  }
158  // sort the closure (boundary analysis expects it this way)
159  std::sort(closure.begin(), closure.end(), stk_classic::mesh::EntityLess());
160 
161  // Run the bounary analysis!
162  stk_classic::mesh::EntitySideVector boundary;
163  stk_classic::mesh::boundary_analysis(bulk_data, closure, element_rank, boundary);
164  STKUNIT_EXPECT_TRUE(!boundary.empty());
165 
166  // Prepare the expected-results as a vector of pairs of pairs representing
167  // ( inside, outside ) where
168  // inside is ( element-id, side-ordinal ) of element inside the closure
169  // outside is ( element-id, side-ordinal ) of element outside the closure
170  // and inside and outside border each other
171 
172  typedef std::pair<stk_classic::mesh::EntityId, stk_classic::mesh::Ordinal> BoundaryItem;
173  typedef std::pair<BoundaryItem, BoundaryItem> BoundaryPair;
174 
175  // Note that certain sides of elements 7, 11, 15 have a boundary with
176  // a shell AND the adjacent element outside the closure.
177 
178  BoundaryPair results[] = {
179  BoundaryPair(BoundaryItem(6, 0), BoundaryItem(5, 2)),
180 
181  BoundaryPair(BoundaryItem(6, 3), BoundaryItem(2, 1)),
182 
183  BoundaryPair(BoundaryItem(7, 2), BoundaryItem(8, 0)),
184  BoundaryPair(BoundaryItem(7, 2), BoundaryItem(43, 0)),
185 
186  BoundaryPair(BoundaryItem(7, 3), BoundaryItem(3, 1)),
187 
188  BoundaryPair(BoundaryItem(10, 0), BoundaryItem(9, 2)),
189 
190  BoundaryPair(BoundaryItem(11, 2), BoundaryItem(12, 0)),
191  BoundaryPair(BoundaryItem(11, 2), BoundaryItem(44, 0)),
192 
193  BoundaryPair(BoundaryItem(14, 0), BoundaryItem(13, 2)),
194 
195  BoundaryPair(BoundaryItem(14, 1), BoundaryItem(0, 0)),
196 
197  BoundaryPair(BoundaryItem(15, 1), BoundaryItem(0, 0)),
198 
199  BoundaryPair(BoundaryItem(15, 2), BoundaryItem(16, 0)),
200  BoundaryPair(BoundaryItem(15, 2), BoundaryItem(45, 0))
201  };
202 
203  // Convert the boundary returned by boundary_analysis into a data-structure
204  // comparable to expected_results
205 
206  BoundaryPair expected_results[sizeof(results)/sizeof(BoundaryPair)];
207 
208  unsigned i = 0;
209  stk_classic::mesh::EntitySideVector::iterator itr = boundary.begin();
210 
211  for (; itr != boundary.end(); ++itr, ++i)
212  {
213  stk_classic::mesh::EntitySide& side = *itr;
214  stk_classic::mesh::EntitySideComponent& inside_closure = side.inside;
215  stk_classic::mesh::EntityId inside_id = inside_closure.entity != NULL ? inside_closure.entity->identifier() : 0;
216  stk_classic::mesh::EntityId inside_side = inside_closure.entity != NULL ? inside_closure.side_ordinal : 0;
217  stk_classic::mesh::EntitySideComponent& outside_closure = side.outside;
218  stk_classic::mesh::EntityId outside_id = outside_closure.entity != NULL ? outside_closure.entity->identifier() : 0;
219  stk_classic::mesh::EntityId outside_side = outside_closure.entity != NULL ? outside_closure.side_ordinal : 0;
220 
221  expected_results[i] = BoundaryPair(BoundaryItem(inside_id, inside_side),
222  BoundaryItem(outside_id, outside_side));
223  }
224 
225  // Check that results match expected results
226 
227  STKUNIT_EXPECT_EQ(sizeof(results), sizeof(expected_results));
228 
229  for (i = 0; i < sizeof(results)/sizeof(BoundaryPair); ++i) {
230  STKUNIT_EXPECT_TRUE(results[i] == expected_results[i]);
231  }
232 }
233 
234 void UnitTestStkMeshBoundaryAnalysis::test_boundary_analysis_null_topology()
235 {
236  //test on boundary_analysis for closure with a NULL topology - coverage of lines 39-40 of BoundaryAnalysis.cpp
237 
238  //create new fem_meta, bulk and boundary for this test
239  const int spatial_dimension = 3;
241  fem_meta.FEM_initialize(spatial_dimension, stk_classic::mesh::fem::entity_rank_names(spatial_dimension));
242 
243  const stk_classic::mesh::EntityRank side_rank = fem_meta.side_rank();
244 
245  //declare part with topology = NULL
246  stk_classic::mesh::Part & quad_part = fem_meta.declare_part("quad_part", side_rank);
247  fem_meta.commit();
248 
249  stk_classic::ParallelMachine comm(MPI_COMM_WORLD);
251  stk_classic::mesh::BulkData bulk ( meta , comm , 100 );
252 
253  stk_classic::mesh::EntitySideVector boundary;
254  std::vector<stk_classic::mesh::Entity*> newclosure;
255 
257  face_parts.push_back(&quad_part);
258 
259  bulk.modification_begin();
260  if (m_rank == 0) {
261  stk_classic::mesh::Entity & new_face = bulk.declare_entity(side_rank, 1, face_parts);
262  newclosure.push_back(&new_face);
263  }
264 
265  stk_classic::mesh::boundary_analysis(bulk, newclosure, side_rank, boundary);
266  /*
267  STKUNIT_EXPECT_TRUE(!boundary.empty());
268  */
269 
270  bulk.modification_end();
271 }
Comparison operator for entities compares the entities&#39; keys.
Definition: Entity.hpp:478
Part & locally_owned_part() const
Subset for the problem domain that is owned by the local process. Ghost entities are not members of t...
void declare_relation(Entity &e_from, Entity &e_to, const RelationIdentifier local_id)
Declare a relation and its converse between entities in the same mesh.
FEMMetaData is a class that implements a Finite Element Method skin on top of the Sierra Tool Kit Met...
Definition: FEMMetaData.hpp:54
The manager of an integrated collection of parts and fields.
Definition: MetaData.hpp:56
EntityRank side_rank() const
Returns the side rank which changes depending on spatial dimension.
EntityRank element_rank() const
Returns the element rank which is always equal to spatial dimension.
This is a class for selecting buckets based on a set of meshparts and set logic.
Definition: Selector.hpp:112
Entity * get_entity(EntityRank entity_rank, EntityId entity_id) const
Get entity with a given key.
Definition: BulkData.hpp:211
An application-defined subset of a problem domain.
Definition: Part.hpp:49
unsigned parallel_machine_rank(ParallelMachine parallel_machine)
Member function parallel_machine_rank ...
Definition: Parallel.cpp:29
Part & declare_part(const std::string &name, fem::CellTopology cell_topology)
Declare a part with a given cell topology.
bool modification_end()
Parallel synchronization of modifications and transition to the guaranteed parallel consistent state...
bool modification_begin()
Begin a modification phase during which the mesh bulk data could become parallel inconsistent. This is a parallel synchronous call. The first time this method is called the mesh meta data is verified to be committed and parallel consistent. An exception is thrown if this verification fails.
Definition: BulkData.cpp:172
unsigned parallel_machine_size(ParallelMachine parallel_machine)
Member function parallel_machine_size ...
Definition: Parallel.cpp:18
Manager for an integrated collection of entities, entity relations, and buckets of field data...
Definition: BulkData.hpp:49
static MetaData & get_meta_data(FEMMetaData &fem_meta)
Getter for MetaData off of a FEMMetaData object.
void commit()
Commit the part and field declarations so that the meta data manager can be used to create mesh bulk ...
A fundamental unit within the discretization of a problem domain, including but not limited to nodes...
Definition: Entity.hpp:120
void count_entities(const Selector &selector, const BulkData &mesh, std::vector< EntityRank > &count)
Local count selected entities of each type.
MPI_Comm ParallelMachine
Definition: Parallel.hpp:32
Entity & declare_entity(EntityRank ent_rank, EntityId ent_id, const PartVector &parts)
Create or retrieve a locally owned entity of a given rank and id.
Definition: BulkData.cpp:215
std::vector< Part *> PartVector
Collections of parts are frequently maintained as a vector of Part pointers.
Definition: Types.hpp:31
eastl::iterator_traits< InputIterator >::difference_type count(InputIterator first, InputIterator last, const T &value)
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