Sierra Toolkit Version of the Day
LinsysFunctions.cpp
00001 /*------------------------------------------------------------------------*/
00002 /*                 Copyright 2010, 2011 Sandia Corporation.                     */
00003 /*  Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive   */
00004 /*  license for use of this work by or on behalf of the U.S. Government.  */
00005 /*  Export of this program may require a license from the                 */
00006 /*  United States Government.                                             */
00007 /*------------------------------------------------------------------------*/
00008 
00009 
00010 #include <stk_linsys/FeiBaseIncludes.hpp>
00011 #include <stk_linsys/LinsysFunctions.hpp>
00012 #include <stk_linsys/ImplDetails.hpp>
00013 
00014 #include <stk_mesh/base/GetBuckets.hpp>
00015 #include <stk_mesh/base/FieldData.hpp>
00016 
00017 #include <fei_trilinos_macros.hpp>
00018 #include <fei_Solver_AztecOO.hpp>
00019 #include <fei_Trilinos_Helpers.hpp>
00020 
00021 
00022 namespace stk {
00023 namespace linsys {
00024 
00025 void add_connectivities(stk::linsys::LinearSystemInterface& ls,
00026                         stk::mesh::EntityRank entity_rank,
00027                         stk::mesh::EntityRank connected_entity_rank,
00028                         const stk::mesh::FieldBase& field,
00029                         const stk::mesh::Selector& selector,
00030                         const stk::mesh::BulkData& mesh_bulk)
00031 {
00032   const std::vector<mesh::Bucket*>& mesh_buckets = mesh_bulk.buckets(entity_rank);
00033   std::vector<mesh::Bucket*> part_buckets;
00034   stk::mesh::get_buckets(selector, mesh_buckets, part_buckets);
00035 
00036   if (part_buckets.empty()) return;
00037 
00038   DofMapper& dof_mapper = ls.get_DofMapper();
00039 
00040   dof_mapper.add_dof_mappings(mesh_bulk, selector, connected_entity_rank, field);
00041 
00042   int field_id = dof_mapper.get_field_id(field);
00043 
00044   stk::mesh::Entity& first_entity = *(part_buckets[0]->begin());
00045   stk::mesh::PairIterRelation rel = first_entity.relations(connected_entity_rank);
00046   int num_connected = rel.second - rel.first;
00047 
00048   fei::SharedPtr<fei::MatrixGraph> matgraph = ls.get_fei_MatrixGraph();
00049 
00050   int pattern_id = matgraph->definePattern(num_connected, connected_entity_rank, field_id);
00051 
00052   int num_entities = 0;
00053   for(size_t i=0; i<part_buckets.size(); ++i) {
00054     num_entities += part_buckets[i]->size();
00055   }
00056 
00057   int block_id = matgraph->initConnectivityBlock(num_entities, pattern_id);
00058 
00059   std::vector<int> connected_ids(num_connected);
00060 
00061   for(size_t i=0; i<part_buckets.size(); ++i) {
00062     stk::mesh::Bucket::iterator
00063       b_iter = part_buckets[i]->begin(),
00064       b_end  = part_buckets[i]->end();
00065     for(; b_iter != b_end; ++b_iter) {
00066       stk::mesh::Entity& entity = *b_iter;
00067       rel = entity.relations(connected_entity_rank);
00068       for(int j=0; rel.first != rel.second; ++rel.first, ++j) {
00069         connected_ids[j] = rel.first->entity()->identifier();
00070       }
00071       int conn_id = entity.identifier();
00072       matgraph->initConnectivity(block_id, conn_id, &connected_ids[0]);
00073     }
00074   }
00075 }
00076 
00077 void dirichlet_bc(stk::linsys::LinearSystemInterface& ls,
00078                   const stk::mesh::BulkData& mesh,
00079                   const stk::mesh::Part& bcpart,
00080                   stk::mesh::EntityRank entity_rank,
00081                   const stk::mesh::FieldBase& field,
00082                   unsigned field_component,
00083                   double prescribed_value)
00084 {
00085   std::vector<stk::mesh::Bucket*> buckets;
00086   stk::mesh::Selector selector(bcpart);
00087   stk::mesh::get_buckets(selector, mesh.buckets(entity_rank), buckets);
00088 
00089   int field_id = ls.get_DofMapper().get_field_id(field);
00090   std::vector<int> entity_ids;
00091 
00092   for(size_t i=0; i<buckets.size(); ++i) {
00093     stk::mesh::Bucket::iterator
00094       iter = buckets[i]->begin(), iend = buckets[i]->end();
00095     for(; iter!=iend; ++iter) {
00096       const stk::mesh::Entity& entity = *iter;
00097       entity_ids.push_back(stk::linsys::impl::entityid_to_int(entity.identifier()));
00098     }
00099   }
00100 
00101   std::vector<double> prescribed_values(entity_ids.size(), prescribed_value);
00102 
00103   ls.get_fei_LinearSystem()->loadEssentialBCs(entity_ids.size(),
00104                                        &entity_ids[0],
00105                                        stk::linsys::impl::entitytype_to_int(entity_rank),
00106                                        field_id, field_component,
00107                                        &prescribed_values[0]);
00108 }
00109 
00110 int fei_solve(int & status, fei::LinearSystem &fei_ls, const Teuchos::ParameterList & params )
00111 {
00112 //Note: hard-coded to Trilinos/Aztec!!!
00113   Solver_AztecOO solver_az;
00114 
00115   fei::ParameterSet param_set;
00116 
00117   Trilinos_Helpers::copy_parameterlist(params, param_set);
00118 
00119   int iterationsTaken = 0;
00120 
00121   return solver_az.solve( & fei_ls,
00122                      NULL,
00123                      param_set,
00124                      iterationsTaken,
00125                      status
00126                   );
00127 }
00128 
00129 double compute_residual_norm2(fei::LinearSystem& fei_ls, fei::Vector& r)
00130 {
00131   fei::SharedPtr<fei::Matrix> A = fei_ls.getMatrix();
00132   fei::SharedPtr<fei::Vector> x = fei_ls.getSolutionVector();
00133   fei::SharedPtr<fei::Vector> b = fei_ls.getRHS();
00134 
00135   //form r = A*x
00136   A->multiply(x.get(), &r);
00137   //form r = b - r   (i.e., r = b - A*x)
00138   r.update(1, b.get(), -1);
00139 
00141 //terrible data copy. fei::Vector should provide a norm operation instead
00142 //of making me roll my own here...
00143 
00144   std::vector<int> indices;
00145   r.getVectorSpace()->getIndices_Owned(indices);
00146   std::vector<double> coefs(indices.size());
00147   r.copyOut(indices.size(), &indices[0], &coefs[0]);
00148   double local_sum = 0;
00149   for(size_t i=0; i<indices.size(); ++i) {
00150     local_sum += coefs[i]*coefs[i];
00151   }
00152 #ifdef HAVE_MPI
00153   MPI_Comm comm = r.getVectorSpace()->getCommunicator();
00154   double global_sum = 0;
00155   int num_doubles = 1;
00156   MPI_Allreduce(&local_sum, &global_sum, num_doubles, MPI_DOUBLE, MPI_SUM, comm);
00157 #else
00158   double global_sum = local_sum;
00159 #endif
00160   return std::sqrt(global_sum);
00161 }
00162 
00163 void copy_vector_to_mesh( fei::Vector & vec,
00164                           const DofMapper & dof,
00165                           stk::mesh::BulkData & mesh_bulk_data
00166                         )
00167 {
00168   vec.scatterToOverlap();
00169 
00170   std::vector<int> shared_and_owned_indices;
00171 
00172   vec.getVectorSpace()->getIndices_SharedAndOwned(shared_and_owned_indices);
00173 
00174   size_t num_values = shared_and_owned_indices.size();
00175 
00176   if(num_values == 0) {
00177     return;
00178   }
00179 
00180   std::vector<double> values(num_values);
00181   vec.copyOut(num_values,&shared_and_owned_indices[0],&values[0]);
00182 
00183   stk::mesh::EntityRank ent_type;
00184   stk::mesh::EntityId ent_id;
00185   const stk::mesh::FieldBase * field;
00186   int offset_into_field;
00187 
00188   for(size_t i = 0; i < num_values; ++i)
00189   {
00190 
00191     dof.get_dof( shared_and_owned_indices[i],
00192                  ent_type,
00193                  ent_id,
00194                  field,
00195                  offset_into_field
00196                 );
00197 
00198     stk::mesh::Entity & entity = *mesh_bulk_data.get_entity(ent_type, ent_id);
00199 
00200     void * data = stk::mesh::field_data(*field,entity);
00201 
00202     if(!(field->type_is<double>()) || data == NULL) {
00203       std::ostringstream oss;
00204       oss << "stk::linsys::copy_vector_to_mesh ERROR, bad data type, or ";
00205       oss << " field (" << field->name() << ") not found on entity with type " << entity.entity_rank();
00206       oss << " and ID " << entity.identifier();
00207       std::string str = oss.str();
00208       throw std::runtime_error(str.c_str());
00209     }
00210 
00211     double * double_data = reinterpret_cast<double *>(data);
00212 
00213     double_data[offset_into_field] = values[i];
00214 
00215   }
00216 }
00217 
00218 void scale_matrix(double scalar, fei::Matrix& matrix)
00219 {
00220   fei::SharedPtr<fei::VectorSpace> vspace = matrix.getMatrixGraph()->getRowSpace();
00221 
00222   int numRows = vspace->getNumIndices_Owned();
00223   std::vector<int> rows(numRows);
00224   vspace->getIndices_Owned(numRows, &rows[0], numRows);
00225 
00226   std::vector<int> indices;
00227   std::vector<double> coefs;
00228 
00229   for(size_t i=0; i<rows.size(); ++i) {
00230     int rowlen = 0;
00231     matrix.getRowLength(rows[i], rowlen);
00232 
00233     if ((int)indices.size() < rowlen) {
00234       indices.resize(rowlen);
00235       coefs.resize(rowlen);
00236     }
00237 
00238     matrix.copyOutRow(rows[i], rowlen, &coefs[0], &indices[0]);
00239 
00240     for(int j=0; j<rowlen; ++j) {
00241       coefs[j] *= scalar;
00242     }
00243 
00244     double* coefPtr = &coefs[0];
00245     matrix.copyIn(1, &rows[i], rowlen, &indices[0], &coefPtr);
00246   }
00247 }
00248 
00249 void add_matrix_to_matrix(double scalar,
00250                           const fei::Matrix& src_matrix,
00251                           fei::Matrix& dest_matrix)
00252 {
00253   fei::SharedPtr<fei::VectorSpace> vspace = src_matrix.getMatrixGraph()->getRowSpace();
00254 
00255   int numRows = vspace->getNumIndices_Owned();
00256   std::vector<int> rows(numRows);
00257   vspace->getIndices_Owned(numRows, &rows[0], numRows);
00258 
00259   std::vector<int> indices;
00260   std::vector<double> coefs;
00261 
00262   for(size_t i=0; i<rows.size(); ++i) {
00263     int rowlen = 0;
00264     src_matrix.getRowLength(rows[i], rowlen);
00265 
00266     if ((int)indices.size() < rowlen) {
00267       indices.resize(rowlen);
00268       coefs.resize(rowlen);
00269     }
00270 
00271     src_matrix.copyOutRow(rows[i], rowlen, &coefs[0], &indices[0]);
00272 
00273     for(int j=0; j<rowlen; ++j) {
00274       coefs[j] *= scalar;
00275     }
00276 
00277     double* coefPtr = &coefs[0];
00278     dest_matrix.sumIn(1, &rows[i], rowlen, &indices[0], &coefPtr);
00279   }
00280 }
00281 
00282 void scale_vector(double scalar,
00283                   fei::Vector& vec)
00284 {
00285   fei::SharedPtr<fei::VectorSpace> vspace = vec.getVectorSpace();
00286 
00287   int numIndices = vspace->getNumIndices_Owned();
00288   std::vector<int> indices(numIndices);
00289   vspace->getIndices_Owned(numIndices, &indices[0], numIndices);
00290 
00291   std::vector<double> coefs(numIndices);
00292 
00293   vec.copyOut(numIndices, &indices[0], &coefs[0]);
00294 
00295   for(size_t j=0; j<coefs.size(); ++j) {
00296     coefs[j] *= scalar;
00297   }
00298 
00299   vec.copyIn(numIndices, &indices[0], &coefs[0]);
00300 }
00301 
00302 void add_vector_to_vector(double scalar,
00303                           const fei::Vector& src_vector,
00304                           fei::Vector& dest_vector)
00305 {
00306   fei::SharedPtr<fei::VectorSpace> vspace = src_vector.getVectorSpace();
00307 
00308   int numIndices = vspace->getNumIndices_Owned();
00309   std::vector<int> indices(numIndices);
00310   vspace->getIndices_Owned(numIndices, &indices[0], numIndices);
00311 
00312   std::vector<double> coefs(numIndices);
00313 
00314   src_vector.copyOut(numIndices, &indices[0], &coefs[0]);
00315 
00316   for(size_t j=0; j<coefs.size(); ++j) {
00317     coefs[j] *= scalar;
00318   }
00319 
00320   dest_vector.sumIn(numIndices, &indices[0], &coefs[0]);
00321 }
00322 
00323 }//namespace linsys
00324 }//namespace stk
00325 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends