Sierra Toolkit  Version of the Day
QuadFixture.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 <algorithm>
10 
11 #include <stk_util/environment/ReportHandler.hpp>
12 
13 #include <stk_mesh/fixtures/QuadFixture.hpp>
14 
15 #include <stk_mesh/base/FieldData.hpp>
16 #include <stk_mesh/base/Types.hpp>
17 #include <stk_mesh/base/Entity.hpp>
18 #include <stk_mesh/base/BulkModification.hpp>
19 
20 #include <stk_mesh/fem/Stencils.hpp>
21 #include <stk_mesh/fem/BoundaryAnalysis.hpp>
22 #include <stk_mesh/fem/FEMMetaData.hpp>
23 #include <stk_mesh/fem/FEMHelpers.hpp>
24 
25 namespace stk_classic {
26 namespace mesh {
27 namespace fixtures {
28 
29 QuadFixture::QuadFixture( stk_classic::ParallelMachine pm ,
30  unsigned nx , unsigned ny )
31  : m_spatial_dimension(2),
32  m_fem_meta( m_spatial_dimension, fem::entity_rank_names(m_spatial_dimension) ),
33  m_bulk_data( stk_classic::mesh::fem::FEMMetaData::get_meta_data(m_fem_meta) , pm ),
34  m_quad_part( fem::declare_part<shards::Quadrilateral<4> >(m_fem_meta, "quad_part" ) ),
35  m_coord_field( m_fem_meta.declare_field<CoordFieldType>("Coordinates") ),
36  m_coord_gather_field( m_fem_meta.declare_field<CoordGatherFieldType>("GatherCoordinates") ),
37  m_nx( nx ),
38  m_ny( ny )
39 {
40  typedef shards::Quadrilateral<4> Quad4 ;
41  const unsigned nodes_per_elem = Quad4::node_count;
42 
43  //put coord-field on all nodes:
44  put_field(
45  m_coord_field,
46  fem::FEMMetaData::NODE_RANK,
47  m_fem_meta.universal_part(),
48  m_spatial_dimension
49  );
50 
51  //put coord-gather-field on all elements:
52  put_field(
53  m_coord_gather_field,
54  m_fem_meta.element_rank(),
55  m_fem_meta.universal_part(),
56  nodes_per_elem
57  );
58 
59  // Field relation so coord-gather-field on elements points
60  // to coord-field of the element's nodes
61  m_fem_meta.declare_field_relation(
62  m_coord_gather_field,
63  fem::element_node_stencil<Quad4, 2>,
64  m_coord_field
65  );
66 }
67 
68 void QuadFixture::node_x_y( EntityId entity_id, unsigned &x , unsigned &y ) const
69 {
70  entity_id -= 1;
71 
72  x = entity_id % (m_nx+1);
73  entity_id /= (m_nx+1);
74 
75  y = entity_id;
76 }
77 
78 void QuadFixture::elem_x_y( EntityId entity_id, unsigned &x , unsigned &y ) const
79 {
80  entity_id -= 1;
81 
82  x = entity_id % m_nx;
83  entity_id /= m_nx;
84 
85  y = entity_id;
86 }
87 
88 
90 {
91  std::vector<EntityId> element_ids_on_this_processor;
92 
93  const unsigned p_size = m_bulk_data.parallel_size();
94  const unsigned p_rank = m_bulk_data.parallel_rank();
95  const unsigned num_elems = m_nx * m_ny;
96 
97  const EntityId beg_elem = 1 + ( num_elems * p_rank ) / p_size ;
98  const EntityId end_elem = 1 + ( num_elems * ( p_rank + 1 ) ) / p_size ;
99 
100  for ( EntityId i = beg_elem; i != end_elem; ++i) {
101  element_ids_on_this_processor.push_back(i);
102  }
103 
104  generate_mesh(element_ids_on_this_processor);
105 }
106 
107 void QuadFixture::generate_mesh(std::vector<EntityId> & element_ids_on_this_processor)
108 {
109  {
110  //sort and unique the input elements
111  std::vector<EntityId>::iterator ib = element_ids_on_this_processor.begin();
112  std::vector<EntityId>::iterator ie = element_ids_on_this_processor.end();
113 
114  std::sort( ib, ie);
115  ib = std::unique( ib, ie);
116  element_ids_on_this_processor.erase(ib, ie);
117  }
118 
119  m_bulk_data.modification_begin();
120 
121  {
122  // Declare the elements that belong on this process
123 
124  std::vector<EntityId>::const_iterator ib = element_ids_on_this_processor.begin();
125  const std::vector<EntityId>::const_iterator ie = element_ids_on_this_processor.end();
126  for (; ib != ie; ++ib) {
127  EntityId entity_id = *ib;
128  unsigned ix = 0, iy = 0;
129  elem_x_y(entity_id, ix, iy);
130 
131  stk_classic::mesh::EntityId elem_nodes[4] ;
132 
133  elem_nodes[0] = node_id( ix , iy );
134  elem_nodes[1] = node_id( ix+1 , iy );
135  elem_nodes[2] = node_id( ix+1 , iy+1 );
136  elem_nodes[3] = node_id( ix , iy+1 );
137 
138  stk_classic::mesh::fem::declare_element( m_bulk_data, m_quad_part, elem_id( ix , iy ) , elem_nodes);
139  for (unsigned i = 0; i<4; ++i) {
140  stk_classic::mesh::Entity * const node = m_bulk_data.get_entity( fem::FEMMetaData::NODE_RANK , elem_nodes[i] );
141 
142  ThrowRequireMsg( node != NULL,
143  "This process should know about the nodes that make up its element");
144 
145  // Compute and assign coordinates to the node
146  unsigned nx = 0, ny = 0;
147  node_x_y(elem_nodes[i], nx, ny);
148 
149  Scalar * data = stk_classic::mesh::field_data( m_coord_field , *node );
150 
151  data[0] = nx ;
152  data[1] = ny ;
153  }
154  }
155  }
156 
157  m_bulk_data.modification_end();
158 }
159 
160 
161 } // fixtures
162 } // mesh
163 } // stk
FieldTraits< field_type >::data_type * field_data(const field_type &f, const Bucket::iterator i)
Pointer to the field data array.
Definition: FieldData.hpp:116
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
EntityRank element_rank() const
Returns the element rank which is always equal to spatial dimension.
Part & universal_part() const
Universal subset for the problem domain. All other parts are a subset of the universal part...
field_type & put_field(field_type &field, EntityRank entity_rank, const Part &part, const void *init_value=NULL)
Declare a field to exist for a given entity type and Part.
Entity * get_entity(EntityRank entity_rank, EntityId entity_id) const
Get entity with a given key.
Definition: BulkData.hpp:211
Field with defined data type and multi-dimensions (if any)
Definition: Field.hpp:118
EntityId elem_id(unsigned x, unsigned y) const
Definition: QuadFixture.hpp:71
Entity * node(unsigned x, unsigned y) const
Definition: QuadFixture.hpp:79
bool modification_end()
Parallel synchronization of modifications and transition to the guaranteed parallel consistent state...
unsigned parallel_size() const
Size of the parallel machine.
Definition: BulkData.hpp:82
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
EntityId entity_id(const EntityKey &key)
Given an entity key, return the identifier for the entity.
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.
void node_x_y(EntityId entity_id, unsigned &x, unsigned &y) const
Definition: QuadFixture.cpp:68
MPI_Comm ParallelMachine
Definition: Parallel.hpp:32
unsigned parallel_rank() const
Rank of the parallel machine&#39;s local processor.
Definition: BulkData.hpp:85
void elem_x_y(EntityId entity_id, unsigned &x, unsigned &y) const
Definition: QuadFixture.cpp:78
void declare_field_relation(PointerFieldType &pointer_field, relation_stencil_ptr stencil, ReferencedFieldType &referenced_field)
Declare a field relation.
EntityId node_id(unsigned x, unsigned y) const
Definition: QuadFixture.hpp:64