Sierra Toolkit  Version of the Day
HexFixture.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/HexFixture.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/FEMHelpers.hpp>
22 #include <stk_mesh/fem/BoundaryAnalysis.hpp>
23 
24 namespace stk_classic {
25 namespace mesh {
26 namespace fixtures {
27 
28 HexFixture::HexFixture(stk_classic::ParallelMachine pm, unsigned nx, unsigned ny, unsigned nz)
29  : m_spatial_dimension(3),
30  m_nx(nx),
31  m_ny(ny),
32  m_nz(nz),
33  m_fem_meta( m_spatial_dimension, fem::entity_rank_names(m_spatial_dimension) ),
34  m_bulk_data( stk_classic::mesh::fem::FEMMetaData::get_meta_data(m_fem_meta) , pm ),
35  m_hex_part( fem::declare_part<shards::Hexahedron<8> >(m_fem_meta, "hex_part") ),
36  m_node_part( m_fem_meta.declare_part("node_part", m_fem_meta.node_rank() ) ),
37  m_coord_field( m_fem_meta.declare_field<CoordFieldType>("Coordinates") ),
38  m_coord_gather_field(m_fem_meta.declare_field<CoordGatherFieldType>("GatherCoordinates"))
39 {
40  typedef shards::Hexahedron<8> Hex8 ;
41  const unsigned nodes_per_elem = Hex8::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  //put coord-gather-field on all elements:
51  put_field(
52  m_coord_gather_field,
53  m_fem_meta.element_rank(),
54  m_fem_meta.universal_part(),
55  nodes_per_elem);
56 
57  // Field relation so coord-gather-field on elements points
58  // to coord-field of the element's nodes
59  m_fem_meta.declare_field_relation( m_coord_gather_field,
60  fem::element_node_stencil<Hex8, 3>,
61  m_coord_field);
62 }
63 
65 {
66  std::vector<EntityId> element_ids_on_this_processor;
67 
68  const unsigned p_size = m_bulk_data.parallel_size();
69  const unsigned p_rank = m_bulk_data.parallel_rank();
70  const unsigned num_elems = m_nx * m_ny * m_nz ;
71 
72  const EntityId beg_elem = 1 + ( num_elems * p_rank ) / p_size ;
73  const EntityId end_elem = 1 + ( num_elems * ( p_rank + 1 ) ) / p_size ;
74 
75  for ( EntityId i = beg_elem; i != end_elem; ++i) {
76  element_ids_on_this_processor.push_back(i);
77  }
78 
79  generate_mesh(element_ids_on_this_processor);
80 }
81 
82 void HexFixture::node_x_y_z( EntityId entity_id, unsigned &x , unsigned &y , unsigned &z ) const
83 {
84  entity_id -= 1;
85 
86  x = entity_id % (m_nx+1);
87  entity_id /= (m_nx+1);
88 
89  y = entity_id % (m_ny+1);
90  entity_id /= (m_ny+1);
91 
92  z = entity_id;
93 }
94 
95 void HexFixture::elem_x_y_z( EntityId entity_id, unsigned &x , unsigned &y , unsigned &z ) const
96 {
97  entity_id -= 1;
98 
99  x = entity_id % m_nx;
100  entity_id /= m_nx;
101 
102  y = entity_id % m_ny;
103  entity_id /= m_ny;
104 
105  z = entity_id;
106 }
107 
108 void HexFixture::generate_mesh(std::vector<EntityId> & element_ids_on_this_processor)
109 {
110  {
111  //sort and unique the input elements
112  std::vector<EntityId>::iterator ib = element_ids_on_this_processor.begin();
113  std::vector<EntityId>::iterator ie = element_ids_on_this_processor.end();
114 
115  std::sort( ib, ie);
116  ib = std::unique( ib, ie);
117  element_ids_on_this_processor.erase(ib, ie);
118  }
119 
120  m_bulk_data.modification_begin();
121 
122  {
123  // Declare the elements that belong on this process
124 
125  PartVector add_parts;
126  add_parts.push_back(&m_node_part);
127 
128  std::vector<EntityId>::iterator ib = element_ids_on_this_processor.begin();
129  const std::vector<EntityId>::iterator ie = element_ids_on_this_processor.end();
130  for (; ib != ie; ++ib) {
131  EntityId entity_id = *ib;
132  unsigned ix = 0, iy = 0, iz = 0;
133  elem_x_y_z(entity_id, ix, iy, iz);
134 
135  stk_classic::mesh::EntityId elem_node[8] ;
136 
137  elem_node[0] = node_id( ix , iy , iz );
138  elem_node[1] = node_id( ix+1 , iy , iz );
139  elem_node[2] = node_id( ix+1 , iy , iz+1 );
140  elem_node[3] = node_id( ix , iy , iz+1 );
141  elem_node[4] = node_id( ix , iy+1 , iz );
142  elem_node[5] = node_id( ix+1 , iy+1 , iz );
143  elem_node[6] = node_id( ix+1 , iy+1 , iz+1 );
144  elem_node[7] = node_id( ix , iy+1 , iz+1 );
145 
146  stk_classic::mesh::fem::declare_element( m_bulk_data, m_hex_part, elem_id( ix , iy , iz ) , elem_node);
147 
148  for (unsigned i = 0; i<8; ++i) {
149  stk_classic::mesh::Entity * const node = m_bulk_data.get_entity( fem::FEMMetaData::NODE_RANK , elem_node[i] );
150  m_bulk_data.change_entity_parts(*node, add_parts);
151 
152  ThrowRequireMsg( node != NULL,
153  "This process should know about the nodes that make up its element");
154 
155  // Compute and assign coordinates to the node
156  unsigned nx = 0, ny = 0, nz = 0;
157  node_x_y_z(elem_node[i], nx, ny, nz);
158 
159  Scalar * data = stk_classic::mesh::field_data( m_coord_field , *node );
160 
161  data[0] = nx ;
162  data[1] = ny ;
163  data[2] = -(Scalar)nz ;
164  }
165  }
166  }
167 
168  m_bulk_data.modification_end();
169 }
170 
171 } // fixtures
172 } // mesh
173 } // 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.
EntityId elem_id(unsigned x, unsigned y, unsigned z) const
Definition: HexFixture.hpp:72
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
void change_entity_parts(Entity &entity, const PartVector &add_parts, const PartVector &remove_parts=PartVector())
Change the parallel-locally-owned entity&#39;s part membership by adding and/or removing parts...
Definition: BulkData.hpp:249
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
void elem_x_y_z(EntityId entity_id, unsigned &x, unsigned &y, unsigned &z) const
Definition: HexFixture.cpp:95
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.
MPI_Comm ParallelMachine
Definition: Parallel.hpp:32
EntityId node_id(unsigned x, unsigned y, unsigned z) const
Definition: HexFixture.hpp:64
void node_x_y_z(EntityId entity_id, unsigned &x, unsigned &y, unsigned &z) const
Definition: HexFixture.cpp:82
unsigned parallel_rank() const
Rank of the parallel machine&#39;s local processor.
Definition: BulkData.hpp:85
std::vector< Part *> PartVector
Collections of parts are frequently maintained as a vector of Part pointers.
Definition: Types.hpp:31
Entity * node(unsigned x, unsigned y, unsigned z) const
Definition: HexFixture.hpp:80
void declare_field_relation(PointerFieldType &pointer_field, relation_stencil_ptr stencil, ReferencedFieldType &referenced_field)
Declare a field relation.