Sierra Toolkit  Version of the Day
LinsysFunctions.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 
10 #include <stk_linsys/FeiBaseIncludes.hpp>
11 #include <stk_linsys/LinsysFunctions.hpp>
12 #include <stk_linsys/ImplDetails.hpp>
13 
14 #include <stk_mesh/base/GetBuckets.hpp>
15 #include <stk_mesh/base/FieldData.hpp>
16 
17 #include <fei_trilinos_macros.hpp>
18 #include <fei_Solver_AztecOO.hpp>
19 #include <fei_Trilinos_Helpers.hpp>
20 
21 
22 namespace stk_classic {
23 namespace linsys {
24 
25 void add_connectivities(stk_classic::linsys::LinearSystemInterface& ls,
26  stk_classic::mesh::EntityRank entity_rank,
27  stk_classic::mesh::EntityRank connected_entity_rank,
28  const stk_classic::mesh::FieldBase& field,
29  const stk_classic::mesh::Selector& selector,
30  const stk_classic::mesh::BulkData& mesh_bulk)
31 {
32  const std::vector<mesh::Bucket*>& mesh_buckets = mesh_bulk.buckets(entity_rank);
33  std::vector<mesh::Bucket*> part_buckets;
34  stk_classic::mesh::get_buckets(selector, mesh_buckets, part_buckets);
35 
36  if (part_buckets.empty()) return;
37 
38  DofMapper& dof_mapper = ls.get_DofMapper();
39 
40  dof_mapper.add_dof_mappings(mesh_bulk, selector, connected_entity_rank, field);
41 
42  int field_id = dof_mapper.get_field_id(field);
43 
44  stk_classic::mesh::Entity& first_entity = *(part_buckets[0]->begin());
45  stk_classic::mesh::PairIterRelation rel = first_entity.relations(connected_entity_rank);
46  int num_connected = rel.second - rel.first;
47 
48  fei::SharedPtr<fei::MatrixGraph> matgraph = ls.get_fei_MatrixGraph();
49 
50  int pattern_id = matgraph->definePattern(num_connected, connected_entity_rank, field_id);
51 
52  int num_entities = 0;
53  for(size_t i=0; i<part_buckets.size(); ++i) {
54  num_entities += part_buckets[i]->size();
55  }
56 
57  int block_id = matgraph->initConnectivityBlock(num_entities, pattern_id);
58 
59  std::vector<int> connected_ids(num_connected);
60 
61  for(size_t i=0; i<part_buckets.size(); ++i) {
62  stk_classic::mesh::Bucket::iterator
63  b_iter = part_buckets[i]->begin(),
64  b_end = part_buckets[i]->end();
65  for(; b_iter != b_end; ++b_iter) {
66  stk_classic::mesh::Entity& entity = *b_iter;
67  rel = entity.relations(connected_entity_rank);
68  for(int j=0; rel.first != rel.second; ++rel.first, ++j) {
69  connected_ids[j] = rel.first->entity()->identifier();
70  }
71  int conn_id = entity.identifier();
72  matgraph->initConnectivity(block_id, conn_id, &connected_ids[0]);
73  }
74  }
75 }
76 
77 void dirichlet_bc(stk_classic::linsys::LinearSystemInterface& ls,
78  const stk_classic::mesh::BulkData& mesh,
79  const stk_classic::mesh::Part& bcpart,
80  stk_classic::mesh::EntityRank entity_rank,
81  const stk_classic::mesh::FieldBase& field,
82  unsigned field_component,
83  double prescribed_value)
84 {
85  std::vector<stk_classic::mesh::Bucket*> buckets;
86  stk_classic::mesh::Selector selector(bcpart);
87  stk_classic::mesh::get_buckets(selector, mesh.buckets(entity_rank), buckets);
88 
89  int field_id = ls.get_DofMapper().get_field_id(field);
90  std::vector<int> entity_ids;
91 
92  for(size_t i=0; i<buckets.size(); ++i) {
93  stk_classic::mesh::Bucket::iterator
94  iter = buckets[i]->begin(), iend = buckets[i]->end();
95  for(; iter!=iend; ++iter) {
96  const stk_classic::mesh::Entity& entity = *iter;
97  entity_ids.push_back(stk_classic::linsys::impl::entityid_to_int(entity.identifier()));
98  }
99  }
100 
101  std::vector<double> prescribed_values(entity_ids.size(), prescribed_value);
102 
103  ls.get_fei_LinearSystem()->loadEssentialBCs(entity_ids.size(),
104  &entity_ids[0],
106  field_id, field_component,
107  &prescribed_values[0]);
108 }
109 
110 int fei_solve(int & status, fei::LinearSystem &fei_ls, const Teuchos::ParameterList & params )
111 {
112 //Note: hard-coded to Trilinos/Aztec!!!
113  Solver_AztecOO solver_az;
114 
115  fei::ParameterSet param_set;
116 
117  Trilinos_Helpers::copy_parameterlist(params, param_set);
118 
119  int iterationsTaken = 0;
120 
121  return solver_az.solve( & fei_ls,
122  NULL,
123  param_set,
124  iterationsTaken,
125  status
126  );
127 }
128 
129 double compute_residual_norm2(fei::LinearSystem& fei_ls, fei::Vector& r)
130 {
131  fei::SharedPtr<fei::Matrix> A = fei_ls.getMatrix();
132  fei::SharedPtr<fei::Vector> x = fei_ls.getSolutionVector();
133  fei::SharedPtr<fei::Vector> b = fei_ls.getRHS();
134 
135  //form r = A*x
136  A->multiply(x.get(), &r);
137  //form r = b - r (i.e., r = b - A*x)
138  r.update(1, b.get(), -1);
139 
141 //terrible data copy. fei::Vector should provide a norm operation instead
142 //of making me roll my own here...
143 
144  std::vector<int> indices;
145  r.getVectorSpace()->getIndices_Owned(indices);
146  std::vector<double> coefs(indices.size());
147  r.copyOut(indices.size(), &indices[0], &coefs[0]);
148  double local_sum = 0;
149  for(size_t i=0; i<indices.size(); ++i) {
150  local_sum += coefs[i]*coefs[i];
151  }
152 #ifdef HAVE_MPI
153  MPI_Comm comm = r.getVectorSpace()->getCommunicator();
154  double global_sum = 0;
155  int num_doubles = 1;
156  MPI_Allreduce(&local_sum, &global_sum, num_doubles, MPI_DOUBLE, MPI_SUM, comm);
157 #else
158  double global_sum = local_sum;
159 #endif
160  return std::sqrt(global_sum);
161 }
162 
163 void copy_vector_to_mesh( fei::Vector & vec,
164  const DofMapper & dof,
165  stk_classic::mesh::BulkData & mesh_bulk_data
166  )
167 {
168  vec.scatterToOverlap();
169 
170  std::vector<int> shared_and_owned_indices;
171 
172  vec.getVectorSpace()->getIndices_SharedAndOwned(shared_and_owned_indices);
173 
174  size_t num_values = shared_and_owned_indices.size();
175 
176  if(num_values == 0) {
177  return;
178  }
179 
180  std::vector<double> values(num_values);
181  vec.copyOut(num_values,&shared_and_owned_indices[0],&values[0]);
182 
183  stk_classic::mesh::EntityRank ent_type;
184  stk_classic::mesh::EntityId ent_id;
185  const stk_classic::mesh::FieldBase * field;
186  int offset_into_field;
187 
188  for(size_t i = 0; i < num_values; ++i)
189  {
190 
191  dof.get_dof( shared_and_owned_indices[i],
192  ent_type,
193  ent_id,
194  field,
195  offset_into_field
196  );
197 
198  stk_classic::mesh::Entity & entity = *mesh_bulk_data.get_entity(ent_type, ent_id);
199 
200  void * data = stk_classic::mesh::field_data(*field,entity);
201 
202  if(!(field->type_is<double>()) || data == NULL) {
203  std::ostringstream oss;
204  oss << "stk_classic::linsys::copy_vector_to_mesh ERROR, bad data type, or ";
205  oss << " field (" << field->name() << ") not found on entity with type " << entity.entity_rank();
206  oss << " and ID " << entity.identifier();
207  std::string str = oss.str();
208  throw std::runtime_error(str.c_str());
209  }
210 
211  double * double_data = reinterpret_cast<double *>(data);
212 
213  double_data[offset_into_field] = values[i];
214 
215  }
216 }
217 
218 void scale_matrix(double scalar, fei::Matrix& matrix)
219 {
220  fei::SharedPtr<fei::VectorSpace> vspace = matrix.getMatrixGraph()->getRowSpace();
221 
222  int numRows = vspace->getNumIndices_Owned();
223  std::vector<int> rows(numRows);
224  vspace->getIndices_Owned(numRows, &rows[0], numRows);
225 
226  std::vector<int> indices;
227  std::vector<double> coefs;
228 
229  for(size_t i=0; i<rows.size(); ++i) {
230  int rowlen = 0;
231  matrix.getRowLength(rows[i], rowlen);
232 
233  if ((int)indices.size() < rowlen) {
234  indices.resize(rowlen);
235  coefs.resize(rowlen);
236  }
237 
238  matrix.copyOutRow(rows[i], rowlen, &coefs[0], &indices[0]);
239 
240  for(int j=0; j<rowlen; ++j) {
241  coefs[j] *= scalar;
242  }
243 
244  double* coefPtr = &coefs[0];
245  matrix.copyIn(1, &rows[i], rowlen, &indices[0], &coefPtr);
246  }
247 }
248 
249 void add_matrix_to_matrix(double scalar,
250  const fei::Matrix& src_matrix,
251  fei::Matrix& dest_matrix)
252 {
253  fei::SharedPtr<fei::VectorSpace> vspace = src_matrix.getMatrixGraph()->getRowSpace();
254 
255  int numRows = vspace->getNumIndices_Owned();
256  std::vector<int> rows(numRows);
257  vspace->getIndices_Owned(numRows, &rows[0], numRows);
258 
259  std::vector<int> indices;
260  std::vector<double> coefs;
261 
262  for(size_t i=0; i<rows.size(); ++i) {
263  int rowlen = 0;
264  src_matrix.getRowLength(rows[i], rowlen);
265 
266  if ((int)indices.size() < rowlen) {
267  indices.resize(rowlen);
268  coefs.resize(rowlen);
269  }
270 
271  src_matrix.copyOutRow(rows[i], rowlen, &coefs[0], &indices[0]);
272 
273  for(int j=0; j<rowlen; ++j) {
274  coefs[j] *= scalar;
275  }
276 
277  double* coefPtr = &coefs[0];
278  dest_matrix.sumIn(1, &rows[i], rowlen, &indices[0], &coefPtr);
279  }
280 }
281 
282 void scale_vector(double scalar,
283  fei::Vector& vec)
284 {
285  fei::SharedPtr<fei::VectorSpace> vspace = vec.getVectorSpace();
286 
287  int numIndices = vspace->getNumIndices_Owned();
288  std::vector<int> indices(numIndices);
289  vspace->getIndices_Owned(numIndices, &indices[0], numIndices);
290 
291  std::vector<double> coefs(numIndices);
292 
293  vec.copyOut(numIndices, &indices[0], &coefs[0]);
294 
295  for(size_t j=0; j<coefs.size(); ++j) {
296  coefs[j] *= scalar;
297  }
298 
299  vec.copyIn(numIndices, &indices[0], &coefs[0]);
300 }
301 
302 void add_vector_to_vector(double scalar,
303  const fei::Vector& src_vector,
304  fei::Vector& dest_vector)
305 {
306  fei::SharedPtr<fei::VectorSpace> vspace = src_vector.getVectorSpace();
307 
308  int numIndices = vspace->getNumIndices_Owned();
309  std::vector<int> indices(numIndices);
310  vspace->getIndices_Owned(numIndices, &indices[0], numIndices);
311 
312  std::vector<double> coefs(numIndices);
313 
314  src_vector.copyOut(numIndices, &indices[0], &coefs[0]);
315 
316  for(size_t j=0; j<coefs.size(); ++j) {
317  coefs[j] *= scalar;
318  }
319 
320  dest_vector.sumIn(numIndices, &indices[0], &coefs[0]);
321 }
322 
323 }//namespace linsys
324 }//namespace stk_classic
325 
void scale_vector(double scalar, fei::Vector &vec)
void add_connectivities(stk_classic::linsys::LinearSystemInterface &ls, stk_classic::mesh::EntityRank entity_rank, stk_classic::mesh::EntityRank connected_entity_rank, const stk_classic::mesh::FieldBase &field, const stk_classic::mesh::Selector &selector, const stk_classic::mesh::BulkData &mesh_bulk)
int entitytype_to_int(stk_classic::mesh::EntityRank entity_rank)
Definition: ImplDetails.cpp:70
int fei_solve(int &status, fei::LinearSystem &fei_ls, const Teuchos::ParameterList &params)
Field base class with an anonymous data type and anonymous multi-dimension.
Definition: FieldBase.hpp:53
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
double compute_residual_norm2(fei::LinearSystem &fei_ls, fei::Vector &r)
void get_dof(int global_index, stk_classic::mesh::EntityRank &ent_type, stk_classic::mesh::EntityId &ent_id, const stk_classic::mesh::FieldBase *&field, int &offset_into_field) const
Definition: DofMapper.cpp:169
This is a class for selecting buckets based on a set of meshparts and set logic.
Definition: Selector.hpp:112
const std::vector< Bucket * > & buckets(EntityRank rank) const
Query all buckets of a given entity rank.
Definition: BulkData.hpp:195
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
void dirichlet_bc(stk_classic::linsys::LinearSystemInterface &ls, const stk_classic::mesh::BulkData &mesh, const stk_classic::mesh::Part &bcpart, stk_classic::mesh::EntityRank entity_rank, const stk_classic::mesh::FieldBase &field, unsigned field_component, double prescribed_value)
int entityid_to_int(stk_classic::mesh::EntityId id)
Definition: ImplDetails.cpp:79
void copy_vector_to_mesh(fei::Vector &vec, const DofMapper &dof, stk_classic::mesh::BulkData &mesh_bulk_data)
PairIterRelation relations() const
All Entity relations for which this entity is a member. The relations are ordered from lowest entity-...
Definition: Entity.hpp:161
int get_field_id(const stk_classic::mesh::FieldBase &field) const
Definition: DofMapper.cpp:136
Manager for an integrated collection of entities, entity relations, and buckets of field data...
Definition: BulkData.hpp:49
void add_dof_mappings(const stk_classic::mesh::BulkData &mesh_bulk, const stk_classic::mesh::Selector &selector, stk_classic::mesh::EntityRank ent_type, const stk_classic::mesh::FieldBase &field)
Definition: DofMapper.cpp:50
void add_matrix_to_matrix(double scalar, const fei::Matrix &src_matrix, fei::Matrix &dest_matrix)
A fundamental unit within the discretization of a problem domain, including but not limited to nodes...
Definition: Entity.hpp:120
Sierra Toolkit.
void add_vector_to_vector(double scalar, const fei::Vector &src_vector, fei::Vector &dest_vector)
AllSelectedBucketsRange get_buckets(const Selector &selector, const BulkData &mesh)
Definition: GetBuckets.cpp:26
void scale_matrix(double scalar, fei::Matrix &matrix)
EntityRank entity_rank(const EntityKey &key)
Given an entity key, return an entity type (rank).