Sierra Toolkit  Version of the Day
UseCase_Rebal_3.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 <use_cases/UseCase_Rebal_3.hpp>
10 
11 #include <stk_util/parallel/Parallel.hpp>
12 #include <stk_util/parallel/ParallelReduce.hpp>
13 
14 #include <stk_mesh/base/FieldData.hpp>
15 #include <stk_mesh/base/GetEntities.hpp>
16 #include <stk_mesh/base/MetaData.hpp>
17 #include <stk_mesh/base/BulkData.hpp>
18 
19 #include <stk_mesh/fem/FEMMetaData.hpp>
20 #include <stk_mesh/fem/FEMHelpers.hpp>
21 #include <stk_mesh/fem/CoordinateSystems.hpp>
22 #include <stk_mesh/fem/CreateAdjacentEntities.hpp>
23 
24 #include <stk_mesh/fixtures/HexFixture.hpp>
25 
29 
30 #include <stk_rebalance_utils/RebalanceUtils.hpp>
31 
32 //----------------------------------------------------------------------
33 
34 using namespace stk_classic::mesh::fixtures;
35 
37 
38 namespace stk_classic {
39 namespace rebalance {
40 namespace use_cases {
41 
42 class Test_Case_3_Partition : public stk_classic::rebalance::Zoltan {
43  public :
44  explicit Test_Case_3_Partition(ParallelMachine pm,
45  const unsigned ndim,
46  Teuchos::ParameterList & rebal_region_parameters,
47  std::string parameters_name=default_parameters_name());
48  virtual ~Test_Case_3_Partition();
49  virtual bool partition_dependents_needed() const;
50 };
51 
52 Test_Case_3_Partition::Test_Case_3_Partition(ParallelMachine pm,
53  const unsigned ndim,
54  Teuchos::ParameterList & rebal_region_parameters,
55  std::string parameters_name) :
56  stk_classic::rebalance::Zoltan(pm, ndim, rebal_region_parameters, parameters_name) {}
57 Test_Case_3_Partition::~Test_Case_3_Partition() {}
58 bool Test_Case_3_Partition::partition_dependents_needed() const
59 {return false;} //Do NOT move dependent entities for this case
60 
61 enum { nx = 4, ny = 4 };
62 
63 bool test_contact_surfaces( stk_classic::ParallelMachine comm )
64 {
65  unsigned spatial_dimension = 2;
66  std::vector<std::string> rank_names = stk_classic::mesh::fem::entity_rank_names(spatial_dimension);
67  const stk_classic::mesh::EntityRank constraint_rank = rank_names.size();
68  rank_names.push_back("Constraint");
69 
71  fem_meta.FEM_initialize(spatial_dimension, rank_names);
72 
74  stk_classic::mesh::BulkData bulk_data( meta_data , comm , 100 );
75  const stk_classic::mesh::EntityRank element_rank = fem_meta.element_rank();
76  const stk_classic::mesh::EntityRank side_rank = fem_meta.side_rank();
77  const stk_classic::mesh::EntityRank node_rank = fem_meta.node_rank();
78 
79  stk_classic::mesh::fem::CellTopology quad_top(shards::getCellTopologyData<shards::Quadrilateral<4> >());
80  stk_classic::mesh::Part & quad_part( fem_meta.declare_part("quad", quad_top ) );
81  stk_classic::mesh::fem::CellTopology side_top(shards::getCellTopologyData<shards::Line<2> >());
82  stk_classic::mesh::Part & side_part( fem_meta.declare_part("line", side_top ) );
83  VectorField & coord_field( fem_meta.declare_field< VectorField >( "coordinates" ) );
84  ScalarField & weight_field( fem_meta.declare_field< ScalarField >( "element_weights" ) );
85 
86  stk_classic::mesh::put_field( coord_field , node_rank , fem_meta.universal_part() );
87  stk_classic::mesh::put_field(weight_field , element_rank , fem_meta.universal_part() );
88 
89  fem_meta.commit();
90 
91  //const unsigned p_size = bulk_data.parallel_size();
92  const unsigned p_rank = bulk_data.parallel_rank();
93 
94  bulk_data.modification_begin();
95 
96  if ( !p_rank ) {
97 
98  std::vector<stk_classic::mesh::EntityVector> quads(nx);
99  for ( unsigned ix = 0 ; ix < nx ; ++ix ) quads[ix].resize(ny);
100 
101  const unsigned nnx = nx + 1 ;
102  const unsigned nny = ny + 1 ;
103  unsigned face_id = 1;
104  for ( unsigned iy = 0 ; iy < ny ; ++iy ) {
105  for ( unsigned ix = 0 ; ix < nx ; ++ix ) {
106  stk_classic::mesh::EntityId elem = 1 + ix + iy * nx ;
107  stk_classic::mesh::EntityId nodes[4] ;
108  nodes[0] = 1 + ix + iy * nnx ;
109  nodes[1] = 2 + ix + iy * nnx ;
110  nodes[2] = 2 + ix + ( iy + 1 ) * nnx ;
111  nodes[3] = 1 + ix + ( iy + 1 ) * nnx ;
112 
113  stk_classic::mesh::Entity &q = stk_classic::mesh::fem::declare_element( bulk_data , quad_part , elem , nodes );
114  if (0==ix) {
115  stk_classic::mesh::fem::declare_element_side( bulk_data, face_id, q, 0 /*local side id*/, &side_part);
116  ++face_id;
117  }
118  quads[ix][iy] = &q;
119  }
120  }
121 
122  for ( unsigned iy = 0 ; iy < ny ; ++iy ) {
123  for ( unsigned ix = 0 ; ix < nx ; ++ix ) {
124  stk_classic::mesh::EntityId elem = 1 + ix + iy * nx ;
125  stk_classic::mesh::Entity * e = bulk_data.get_entity( element_rank, elem );
126  double * const e_weight = stk_classic::mesh::field_data( weight_field , *e );
127  *e_weight = 1.0;
128  }
129  }
130  for ( unsigned iy = 0 ; iy <= ny ; ++iy ) {
131  for ( unsigned ix = 0 ; ix <= nx ; ++ix ) {
132  stk_classic::mesh::EntityId nid = 1 + ix + iy * nnx ;
133  stk_classic::mesh::Entity * n = bulk_data.get_entity( node_rank, nid );
134  double * const coord = stk_classic::mesh::field_data( coord_field , *n );
135  coord[0] = .1*ix;
136  coord[1] = .1*iy;
137  coord[2] = 0;
138  }
139  }
140 
141  // Assign constraint relations between nodes at top and bottom of mesh
142  {
143  const unsigned iy_bottom = 0;
144  const unsigned iy_top = ny;
146  for ( unsigned ix = 0 ; ix <= nx ; ++ix ) {
147  stk_classic::mesh::EntityId nid_bottom = 1 + ix + iy_bottom * nnx ;
148  stk_classic::mesh::EntityId nid_top = 1 + ix + iy_top * nnx ;
149  stk_classic::mesh::Entity * n_bottom = bulk_data.get_entity( node_rank, nid_bottom );
150  stk_classic::mesh::Entity * n_top = bulk_data.get_entity( node_rank, nid_top );
151  const stk_classic::mesh::EntityId constraint_entity_id = 1 + ix + nny * nnx;
152  stk_classic::mesh::Entity & c = bulk_data.declare_entity( constraint_rank, constraint_entity_id, add );
153  bulk_data.declare_relation( c , *n_bottom , 0 );
154  bulk_data.declare_relation( c , *n_top , 1 );
155  }
156  } // end snippet
157 
158  }
159 
160  bulk_data.modification_end();
161 
162  stk_classic::mesh::Selector selector(side_part);
163  selector &= fem_meta.locally_owned_part();
164  // Coordinates are passed to support geometric-based load balancing algorithms
165  // Zoltan partition is specialized form a virtual base class, stk_classic::rebalance::Partition.
166  // Other specializations are possible.
167  Teuchos::ParameterList emptyList;
168  stk_classic::rebalance::use_cases::Test_Case_3_Partition zoltan_partition(comm, spatial_dimension, emptyList);
169  stk_classic::rebalance::rebalance(bulk_data, selector, &coord_field, NULL, zoltan_partition, side_rank);
170 
171  const double imbalance_threshold = 1.5;
172  const bool do_rebal = imbalance_threshold < stk_classic::rebalance::check_balance(bulk_data, NULL, side_rank, &selector);
173 
174  if( !p_rank )
175  std::cerr << std::endl
176  << "imbalance_threshold after rebalance = " << imbalance_threshold <<", "<<do_rebal << std::endl;
177 
178  // Check that we satisfy our threshhold
179  bool result = !do_rebal ;
180 
181  return result;
182 }
183 
184 } //namespace use_cases
185 } //namespace rebalance
186 } //namespace stk_classic
187 
188 
Part & locally_owned_part() const
Subset for the problem domain that is owned by the local process. Ghost entities are not members of t...
FEMMetaData is a class that implements a Finite Element Method skin on top of the Sierra Tool Kit Met...
Definition: FEMMetaData.hpp:54
bool rebalance(mesh::BulkData &bulk_data, const mesh::Selector &selector, const VectorField *coord_ref, const ScalarField *elem_weight_ref, Partition &partition, const stk_classic::mesh::EntityRank rank=stk_classic::mesh::InvalidEntityRank)
Rebalance with a Partition object.
Definition: Rebalance.cpp:164
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.
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...
This is a class for selecting buckets based on a set of meshparts and set logic.
Definition: Selector.hpp:112
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.
An application-defined subset of a problem domain.
Definition: Part.hpp:49
Part & declare_part(const std::string &name, fem::CellTopology cell_topology)
Declare a part with a given cell topology.
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 ...
field_type & declare_field(const std::string &name, unsigned number_of_states=1)
Declare a field of the given field_type, test name, and number of states.
A fundamental unit within the discretization of a problem domain, including but not limited to nodes...
Definition: Entity.hpp:120
Sierra Toolkit.
MPI_Comm ParallelMachine
Definition: Parallel.hpp:32
For partitioning of mesh entities over a processing grid.
Static functions for dynamic load balancing.
EntityRank node_rank() const
Returns the node rank, which is always zero.
std::vector< Part *> PartVector
Collections of parts are frequently maintained as a vector of Part pointers.
Definition: Types.hpp:31
Class for implementing Zoltan based rebalancing.
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
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