Sierra Toolkit  Version of the Day
UnitTestTopology.cpp
1 /*------------------------------------------------------------------------*/
2 /* _ Copyright 2010 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 <sstream>
10 #include <stdexcept>
11 
12 #include <stk_util/unit_test_support/stk_utest_macros.hpp>
13 
14 #include <stk_util/parallel/Parallel.hpp>
15 
16 #include <stk_mesh/base/MetaData.hpp>
17 #include <stk_mesh/base/BulkData.hpp>
18 #include <stk_mesh/base/GetEntities.hpp>
19 #include <stk_mesh/base/Field.hpp>
20 #include <stk_mesh/base/FieldData.hpp>
21 #include <stk_mesh/base/Comm.hpp>
22 #include <stk_mesh/base/EntityComm.hpp>
23 #include <stk_mesh/base/Part.hpp>
24 #include <stk_mesh/base/Entity.hpp>
25 #include <stk_mesh/base/GetBuckets.hpp>
26 #include <stk_mesh/base/Bucket.hpp>
27 #include <stk_mesh/base/BulkModification.hpp>
28 #include <stk_mesh/base/Entity.hpp>
29 #include <stk_mesh/base/Bucket.hpp>
30 #include <stk_mesh/base/Ghosting.hpp>
31 
32 #include <stk_mesh/baseImpl/BucketImpl.hpp>
33 
34 #include <stk_mesh/fem/FEMHelpers.hpp>
35 #include <stk_mesh/fem/FEMMetaData.hpp>
36 #include <stk_mesh/fem/BoundaryAnalysis.hpp>
37 
38 #include <Shards_BasicTopologies.hpp>
39 
45 using stk_classic::mesh::EntityRank;
46 using stk_classic::mesh::EntityId;
47 using stk_classic::mesh::EntitySideComponent;
50 using stk_classic::mesh::EntityRank;
52 
53 class TopologyHelpersTestingFixture
54 {
55  public:
56  TopologyHelpersTestingFixture(ParallelMachine pm);
57  ~TopologyHelpersTestingFixture() {}
58 
59  const int spatial_dimension;
60  FEMMetaData meta;
61  BulkData bulk;
62  const EntityRank element_rank;
63  const EntityRank side_rank;
64  Part & generic_element_part;
65  Part & element_tet_part;
66  Part & element_wedge_part;
67  Part & generic_face_part;
68  Part & another_generic_face_part;
69  Part & face_quad_part;
70  Part & another_generic_element_part;
71 
72  EntityId nextEntityId()
73  { return psize*(++entity_id)+prank; }
74 
75  Entity & create_entity( EntityRank rank, Part& part_membership)
76  {
77  PartVector part_intersection;
78  part_intersection.push_back ( &part_membership );
79  return bulk.declare_entity(rank, nextEntityId(), part_intersection);
80  }
81 
82  private:
83  EntityId entity_id;
84  const int psize;
85  const int prank;
86 };
87 
88 TopologyHelpersTestingFixture::TopologyHelpersTestingFixture(ParallelMachine pm)
89  : spatial_dimension( 3 )
90  , meta( spatial_dimension )
91  , bulk( FEMMetaData::get_meta_data(meta), pm, 100 )
92  , element_rank( meta.element_rank())
93  , side_rank( meta.side_rank())
94  , generic_element_part( meta.declare_part("another part", element_rank ) )
95  , element_tet_part( stk_classic::mesh::fem::declare_part<shards::Tetrahedron<4> >( meta, "block_left_1" ) )
96  , element_wedge_part( stk_classic::mesh::fem::declare_part<shards::Wedge<15> >(meta, "block_left_2" ) )
97  , generic_face_part( stk_classic::mesh::fem::declare_part<shards::Quadrilateral<4> >(meta, "A_1" ) )
98  , another_generic_face_part( meta.declare_part("A_2", side_rank ) )
99  , face_quad_part( meta.declare_part("A_3", side_rank ) )
100  , another_generic_element_part( meta.declare_part("B_3", element_rank ) )
101  , entity_id(0u)
102  , psize(bulk.parallel_size())
103  , prank(bulk.parallel_rank())
104 {
105  meta.commit();
106 }
107 
108 namespace {
109 
110 const EntityRank NODE_RANK = FEMMetaData::NODE_RANK;
111 
112 STKUNIT_UNIT_TEST( testTopologyHelpers, get_cell_topology_based_on_part)
113 {
114  TopologyHelpersTestingFixture fix(MPI_COMM_WORLD);
115  fix.bulk.modification_begin();
116  Entity & elem1 = fix.create_entity( fix.side_rank, fix.generic_face_part );
117 
118  PartVector tmp(1);
119  tmp[0] = & fix.face_quad_part;
120  fix.bulk.change_entity_parts ( elem1 , tmp );
121  STKUNIT_ASSERT_EQUAL( stk_classic::mesh::fem::get_cell_topology(elem1).getCellTopologyData(), shards::getCellTopologyData< shards::Quadrilateral<4> >() );
122  fix.bulk.change_entity_parts ( elem1 , tmp );
123  STKUNIT_ASSERT_EQUAL( stk_classic::mesh::fem::get_cell_topology(elem1).getCellTopologyData(), shards::getCellTopologyData< shards::Quadrilateral<4> >() );
124  tmp[0] = & fix.another_generic_face_part;
125  fix.bulk.change_entity_parts ( elem1 , tmp );
126  STKUNIT_ASSERT_EQUAL( stk_classic::mesh::fem::get_cell_topology(elem1).getCellTopologyData(), shards::getCellTopologyData< shards::Quadrilateral<4> >() );
127  STKUNIT_ASSERT_NE( stk_classic::mesh::fem::get_cell_topology( elem1).getCellTopologyData() , shards::getCellTopologyData< shards::Wedge<15> >() );
128 
129  fix.bulk.modification_end();
130 }
131 
132 STKUNIT_UNIT_TEST( testTopologyHelpers, get_cell_topology_multiple_topologies )
133 {
134  // Coverage for get_cell_topology in TopologyHelpers.cpp; (FAILED WITH MULTIPLE LOCAL TOPOLOGIES)
135  TopologyHelpersTestingFixture fix(MPI_COMM_WORLD);
136 
137  fix.bulk.modification_begin();
138  Entity & elem = fix.create_entity( fix.element_rank, fix.generic_element_part );
139  PartVector add_parts;
140  add_parts.push_back( &fix.element_tet_part );
141  add_parts.push_back( &fix.element_wedge_part );
142  fix.bulk.change_entity_parts( elem, add_parts );
143  fix.bulk.modification_end();
144  STKUNIT_ASSERT_THROW( stk_classic::mesh::fem::get_cell_topology( elem ).getCellTopologyData(), std::runtime_error );
145 }
146 
147 // No longer in the public API
148 // STKUNIT_UNIT_TEST( testTopologyHelpers, get_adjacent_entities_trivial )
149 // {
150 // // Element, elem2, has NULL topology
151 // TopologyHelpersTestingFixture fix(MPI_COMM_WORLD);
152 //
153 // if ( 1 == fix.bulk.parallel_size() ) {
154 //
155 // fix.bulk.modification_begin();
156 // Entity & elem2 = fix.create_entity( fix.element_rank, fix.generic_element_part );
157 // fix.bulk.modification_end();
158 //
159 // std::vector<EntitySideComponent> adjacent_entities;
160 // const EntityRank subcell_rank = fix.element_rank;
161 // const EntityId subcell_identifier = 1;
162 // get_adjacent_entities( elem2 , subcell_rank, subcell_identifier, adjacent_entities);
163 // STKUNIT_ASSERT_TRUE( true );
164 // }
165 // }
166 //
167 // STKUNIT_UNIT_TEST( testTopologyHelpers, get_adjacent_entities_invalid )
168 // {
169 // TopologyHelpersTestingFixture fix(MPI_COMM_WORLD);
170 // fix.bulk.modification_begin();
171 // Entity & elem3 = fix.create_entity( fix.element_rank , fix.generic_element_part );
172 //
173 // PartVector add_parts;
174 // add_parts.push_back( & fix.element_tet_part );
175 // fix.bulk.change_entity_parts ( elem3 , add_parts );
176 // fix.bulk.modification_end();
177 // std::vector<EntitySideComponent> adjacent_entities2;
178 // {
179 // const EntityRank invalid_subcell_rank = 4;
180 // const EntityId valid_subcell_identifier = 0;
181 // STKUNIT_ASSERT_THROW(
182 // get_adjacent_entities( elem3 , invalid_subcell_rank, valid_subcell_identifier, adjacent_entities2),
183 // std::invalid_argument
184 // );
185 // }
186 // {
187 // const EntityRank valid_subcell_rank = 1;
188 // const EntityId invalid_subcell_identifier = 8;
189 // STKUNIT_ASSERT_THROW(
190 // get_adjacent_entities( elem3 , valid_subcell_rank, invalid_subcell_identifier, adjacent_entities2),
191 // std::invalid_argument
192 // );
193 // }
194 // }
195 
196 STKUNIT_UNIT_TEST( testTopologyHelpers, declare_element_side_no_topology )
197 {
198  // Coverage for declare_element_side - TopologyHelpers.cpp - "Cannot discern element topology"
199  TopologyHelpersTestingFixture fix(MPI_COMM_WORLD);
200 
201  fix.bulk.modification_begin();
202  Entity & elem4 = fix.create_entity( fix.element_rank , fix.generic_element_part );
203  STKUNIT_ASSERT_THROW(
204  stk_classic::mesh::fem::declare_element_side( fix.bulk, fix.element_rank, elem4, fix.nextEntityId(), &fix.element_wedge_part ),
205  std::runtime_error
206  );
207  fix.bulk.modification_end();
208 
209 
210  {
211  EntityId elem_node[4];
212  elem_node[0] = 1;
213  elem_node[1] = 2;
214  elem_node[2] = 3;
215  elem_node[3] = 4;
216  fix.bulk.modification_begin();
217  // Cannot declare an element without a topology defined
218  STKUNIT_ASSERT_THROW(
219  stk_classic::mesh::fem::declare_element(fix.bulk, fix.generic_element_part, fix.nextEntityId(), elem_node),
220  std::runtime_error
221  );
222  fix.bulk.modification_end();
223  }
224 }
225 
226 STKUNIT_UNIT_TEST( testTopologyHelpers, declare_element_side_wrong_bulk_data)
227 {
228  // Coverage for verify_declare_element_side - in TopologyHelpers.cpp - "BulkData for 'elem' and 'side' are different"
229  TopologyHelpersTestingFixture fix1(MPI_COMM_WORLD);
230 
231  fix1.bulk.modification_begin();
232 
233  TopologyHelpersTestingFixture fix2(MPI_COMM_WORLD);
234  fix2.bulk.modification_begin();
235  Entity & elem4_2 = fix2.create_entity( fix2.element_rank , fix2.generic_element_part );
236  fix2.bulk.modification_end();
237 
238  STKUNIT_ASSERT_THROW(
239  stk_classic::mesh::fem::declare_element_side( fix1.bulk, fix1.element_rank, elem4_2, fix1.nextEntityId(), &fix1.element_wedge_part),
240  std::runtime_error
241  );
242  fix1.bulk.modification_end();
243 }
244 
245 STKUNIT_UNIT_TEST( testTopologyHelpers, declare_element_side_no_topology_2 )
246 {
247  // Coverage for verify_declare_element_side - in TopologyHelpers.cpp - "No element topology found and cell side id exceeds..."
248  TopologyHelpersTestingFixture fix(MPI_COMM_WORLD);
249  fix.bulk.modification_begin();
250 
251  EntityId elem_node[4];
252  elem_node[0] = 1;
253  elem_node[1] = 2;
254  elem_node[2] = 3;
255  elem_node[3] = 4;
256  Entity & element = stk_classic::mesh::fem::declare_element(fix.bulk, fix.element_tet_part, fix.nextEntityId(), elem_node);
257  const CellTopologyData * const elem_top = stk_classic::mesh::fem::get_cell_topology( element ).getCellTopologyData();
258  const EntityId nSideCount = elem_top->side_count + 10 ;
259  STKUNIT_ASSERT_THROW(
260  stk_classic::mesh::fem::declare_element_side( fix.bulk, fix.nextEntityId(), element, nSideCount, &fix.element_tet_part ),
261  std::runtime_error
262  );
263  fix.bulk.modification_end();
264 }
265 
266 STKUNIT_UNIT_TEST( testTopologyHelpers, declare_element_side_full )
267 {
268  // Go all way the through declare_element_side - use new element
269  TopologyHelpersTestingFixture fix(MPI_COMM_WORLD);
270 
271  fix.bulk.modification_begin();
272 
273  EntityId elem_node[4];
274  elem_node[0] = 1;
275  elem_node[1] = 2;
276  elem_node[2] = 3;
277  elem_node[3] = 4;
278 
279  Entity& element = stk_classic::mesh::fem::declare_element(fix.bulk, fix.element_tet_part, fix.nextEntityId(), elem_node );
280 
281  const EntityId zero_side_count = 0;
282  Entity& face2 = stk_classic::mesh::fem::declare_element_side( fix.bulk, fix.nextEntityId(), element, zero_side_count);
283  fix.bulk.modification_end();
284 
285  PairIterRelation rel2 = face2.relations(NODE_RANK);
286 
287  STKUNIT_ASSERT_TRUE( true );
288 }
289 
290 STKUNIT_UNIT_TEST( testTopologyHelpers, element_side_polarity_valid )
291 {
292  // Coverage of element_side_polarity in TopologyHelpers.cpp 168-181 and 200-215
293  TopologyHelpersTestingFixture fix(MPI_COMM_WORLD);
294  EntityId elem_node[4];
295  elem_node[0] = 1;
296  elem_node[1] = 2;
297  elem_node[2] = 3;
298  elem_node[3] = 4;
299 
300  fix.bulk.modification_begin();
301  Entity & element = stk_classic::mesh::fem::declare_element(fix.bulk, fix.element_tet_part, fix.nextEntityId(), elem_node );
302  const EntityId zero_side_count = 0;
303  Entity& face2 = stk_classic::mesh::fem::declare_element_side( fix.bulk, fix.nextEntityId(), element, zero_side_count);
304  fix.bulk.modification_end();
305 
306  const int local_side_id = 0;
307  STKUNIT_ASSERT_TRUE( stk_classic::mesh::fem::element_side_polarity( element, face2, local_side_id) );
308 
309 }
310 
311 STKUNIT_UNIT_TEST( testTopologyHelpers, element_side_polarity_invalid_1 )
312 {
313  TopologyHelpersTestingFixture fix(MPI_COMM_WORLD);
314  EntityId elem_node[4];
315  elem_node[0] = 1;
316  elem_node[1] = 2;
317  elem_node[2] = 3;
318  elem_node[3] = 4;
319 
320  // Coverage of element_side_polarity in TopologyHelpers.cpp
321  {
322  fix.bulk.modification_begin();
323  Entity & element = stk_classic::mesh::fem::declare_element(fix.bulk, fix.element_tet_part, fix.nextEntityId(), elem_node );
324  const EntityId zero_side_count = 0;
325  Entity& face = stk_classic::mesh::fem::declare_element_side( fix.bulk, fix.nextEntityId(), element, zero_side_count);
326  fix.bulk.modification_end();
327 
328  const int invalid_local_side_id = -1;
329  // Hits "Unsuported local_side_id" error condition:
330  STKUNIT_ASSERT_THROW(
331  stk_classic::mesh::fem::element_side_polarity( element, face, invalid_local_side_id),
332  std::runtime_error
333  );
334  }
335 }
336 
337 STKUNIT_UNIT_TEST( testTopologyHelpers, element_side_polarity_invalid_2 )
338 {
339  TopologyHelpersTestingFixture fix(MPI_COMM_WORLD);
340  EntityId elem_node[4];
341  elem_node[0] = 1;
342  elem_node[1] = 2;
343  elem_node[2] = 3;
344  elem_node[3] = 4;
345 
346  // Coverage of element_side_polarity in TopologyHelpers.cpp - NULL = elem_top
347  fix.bulk.modification_begin();
348 
349  PartVector part_intersection;
350  part_intersection.push_back ( &fix.generic_element_part);
351  Entity & element = fix.bulk.declare_entity(fix.element_rank, fix.nextEntityId(), part_intersection);
352  STKUNIT_ASSERT_TRUE( stk_classic::mesh::fem::get_cell_topology( element ).getCellTopologyData() == NULL );
353 
354  Entity & element_with_top = stk_classic::mesh::fem::declare_element(fix.bulk, fix.element_tet_part, fix.nextEntityId(), elem_node );
355  STKUNIT_ASSERT_TRUE( stk_classic::mesh::fem::get_cell_topology( element_with_top ).getCellTopologyData() != NULL );
356 
357  const EntityId zero_side_count = 0;
358  Entity& face_with_top = stk_classic::mesh::fem::declare_element_side( fix.bulk, fix.nextEntityId(), element_with_top, zero_side_count);
359 
360  fix.bulk.modification_end();
361 
362  const int valid_local_side_id = 0;
363  // Hits "Element has no defined topology" error condition:
364  STKUNIT_ASSERT_TRUE( stk_classic::mesh::fem::get_cell_topology( element ).getCellTopologyData() == NULL );
365  STKUNIT_ASSERT_THROW(
366  stk_classic::mesh::fem::element_side_polarity( element, face_with_top, valid_local_side_id),
367  std::runtime_error
368  );
369 
370 }
371 
372 //----------------------------------------------------------------------
373 //----------------------------------------------------------------------
374 
375 }
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
Entity & declare_element(BulkData &mesh, Part &part, const EntityId elem_id, const EntityId node_id[])
Declare an element member of a Part with a CellTopology and nodes conformal to that topology...
Definition: FEMHelpers.cpp:72
int parallel_rank()
function parallel_rank returns the rank of this processor in the current mpi communicator.
Definition: Env.cpp:318
An application-defined subset of a problem domain.
Definition: Part.hpp:49
PairIterRelation relations() const
All Entity relations for which this entity is a member. The relations are ordered from lowest entity-...
Definition: Entity.hpp:161
Manager for an integrated collection of entities, entity relations, and buckets of field data...
Definition: BulkData.hpp:49
EntityId entity_id(const EntityKey &key)
Given an entity key, return the identifier for the entity.
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
Part & declare_part(FEMMetaData &meta_data, const std::string &name)
Declare a part with a given cell topology. This is just a convenient function that wraps FEMMetaData&#39;...
Definition: FEMHelpers.hpp:93
Sierra Toolkit.
MPI_Comm ParallelMachine
Definition: Parallel.hpp:32
int parallel_size()
function parallel_size returns the number of processors in the current mpi communicator.
Definition: Env.cpp:314
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
bool element_side_polarity(const Entity &elem, const Entity &side, int local_side_id)
Determine the polarity of the local side, more efficient if the local_side_id is known.
Definition: FEMHelpers.cpp:401
Entity & declare_element_side(Entity &elem, Entity &side, const unsigned local_side_id, Part *part)
Create (or find) an element side.
Definition: FEMHelpers.cpp:104