LinsysFunctions.cpp

00001 /*------------------------------------------------------------------------*/
00002 /*                 Copyright 2010 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 <stk_mesh/fem/TopologyHelpers.hpp>
00018 
00019 #include <fei_trilinos_macros.hpp>
00020 #include <fei_Solver_AztecOO.hpp>
00021 #include <fei_Trilinos_Helpers.hpp>
00022 
00023 
00024 namespace stk {
00025 namespace linsys {
00026 
00027 void add_connectivities(stk::linsys::LinearSystem& ls,
00028                         stk::mesh::EntityType entity_type,
00029                         stk::mesh::EntityType connected_entity_type,
00030                         const stk::mesh::FieldBase& field,
00031                         const stk::mesh::PartVector& part_intersection,
00032                         const stk::mesh::BulkData& mesh_bulk)
00033 {
00034   const std::vector<mesh::Bucket*>& mesh_buckets = mesh_bulk.buckets(entity_type);
00035   stk::mesh::Selector selector = selectIntersection(part_intersection);
00036   std::vector<mesh::Bucket*> part_buckets;
00037   stk::mesh::get_buckets(selector, mesh_buckets, part_buckets);
00038 
00039   if (part_buckets.empty()) return;
00040 
00041   DofMapper& dof_mapper = ls.get_DofMapper();
00042 
00043   dof_mapper.add_dof_mappings(mesh_bulk, part_intersection, connected_entity_type, field);
00044 
00045   int field_id = dof_mapper.get_field_id(field);
00046 
00047   stk::mesh::Entity& first_entity = *(part_buckets[0]->begin());
00048   stk::mesh::PairIterRelation rel = first_entity.relations(connected_entity_type);
00049   int num_connected = rel.second - rel.first;
00050 
00051   fei::SharedPtr<fei::MatrixGraph> matgraph = ls.get_fei_MatrixGraph();
00052 
00053   int pattern_id = matgraph->definePattern(num_connected, connected_entity_type, field_id);
00054 
00055   int num_entities = 0;
00056   for(size_t i=0; i<part_buckets.size(); ++i) {
00057     num_entities += part_buckets[i]->size();
00058   }
00059 
00060   int part_ordinal_sum = 0;
00061   for(size_t i=0; i<part_intersection.size(); ++i) {
00062     part_ordinal_sum += part_intersection[i]->mesh_meta_data_ordinal();
00063   }
00064 
00065   int block_id = 1000*(part_ordinal_sum+1) + entity_type;
00066 
00067   matgraph->initConnectivityBlock(block_id, num_entities, pattern_id);
00068 
00069   std::vector<int> connected_ids(num_connected);
00070 
00071   for(size_t i=0; i<part_buckets.size(); ++i) {
00072     stk::mesh::Bucket::iterator
00073       b_iter = part_buckets[i]->begin(),
00074       b_end  = part_buckets[i]->end();
00075     for(; b_iter != b_end; ++b_iter) {
00076       stk::mesh::Entity& entity = *b_iter;
00077       rel = entity.relations(connected_entity_type);
00078       for(int j=0; rel.first != rel.second; ++rel.first, ++j) {
00079         connected_ids[j] = rel.first->entity()->identifier();
00080       }
00081       int conn_id = entity.identifier();
00082       matgraph->initConnectivity(block_id, conn_id, &connected_ids[0]);
00083     }
00084   }
00085 }
00086 
00087 void dirichlet_bc(stk::linsys::LinearSystem& ls,
00088                   const stk::mesh::BulkData& mesh,
00089                   const stk::mesh::Part& bcpart,
00090                   stk::mesh::EntityType entity_type,
00091                   const stk::mesh::FieldBase& field,
00092                   unsigned field_component,
00093                   double prescribed_value)
00094 {
00095   std::vector<stk::mesh::Bucket*> buckets;
00096   stk::mesh::Selector selector(bcpart);
00097   stk::mesh::get_buckets(selector, mesh.buckets(entity_type), buckets);
00098 
00099   int field_id = ls.get_DofMapper().get_field_id(field);
00100   std::vector<int> entity_ids;
00101 
00102   for(size_t i=0; i<buckets.size(); ++i) {
00103     stk::mesh::Bucket::iterator
00104       iter = buckets[i]->begin(), iend = buckets[i]->end();
00105     for(; iter!=iend; ++iter) {
00106       const stk::mesh::Entity& entity = *iter;
00107       entity_ids.push_back(stk::linsys::impl::entityid_to_int(entity.identifier()));
00108     }
00109   }
00110 
00111   std::vector<double> prescribed_values(entity_ids.size(), prescribed_value);
00112 
00113   ls.get_fei_LinearSystem()->loadEssentialBCs(entity_ids.size(),
00114                                        &entity_ids[0],
00115                                        stk::linsys::impl::entitytype_to_int(entity_type),
00116                                        field_id, field_component,
00117                                        &prescribed_values[0]);
00118 }
00119 
00120 int fei_solve(int & status, fei::LinearSystem &fei_ls, const Teuchos::ParameterList & params )
00121 {
00122 //Note: hard-coded to Trilinos/Aztec!!!
00123   Solver_AztecOO solver_az;
00124 
00125   fei::ParameterSet param_set;
00126 
00127   Trilinos_Helpers::copy_parameterlist(params, param_set);
00128 
00129   int iterationsTaken = 0;
00130 
00131   return solver_az.solve( & fei_ls,
00132                      NULL,
00133                      param_set,
00134                      iterationsTaken,
00135                      status
00136                   );
00137 }
00138 
00139 void copy_vector_to_mesh( fei::Vector & vec,
00140                           const DofMapper & dof,
00141                           stk::mesh::BulkData & mesh_bulk_data
00142                         )
00143 {
00144   vec.scatterToOverlap();
00145 
00146   std::vector<int> shared_and_owned_indices;
00147 
00148   vec.getVectorSpace()->getIndices_SharedAndOwned(shared_and_owned_indices);
00149 
00150   size_t num_values = shared_and_owned_indices.size();
00151 
00152   if(num_values == 0) {
00153     return;
00154   }
00155 
00156   std::vector<double> values(num_values);
00157   vec.copyOut(num_values,&shared_and_owned_indices[0],&values[0]);
00158 
00159   stk::mesh::EntityType ent_type;
00160   stk::mesh::EntityId ent_id;
00161   const stk::mesh::FieldBase * field;
00162   int offset_into_field;
00163 
00164   for(size_t i = 0; i < num_values; ++i)
00165   {
00166 
00167     dof.get_dof( shared_and_owned_indices[i],
00168                  ent_type,
00169                  ent_id,
00170                  field,
00171                  offset_into_field
00172                 );
00173 
00174     stk::mesh::Entity & entity = *mesh_bulk_data.get_entity(ent_type, ent_id);
00175 
00176     void * data = stk::mesh::field_data(*field,entity);
00177 
00178     if(!(field->type_is<double>()) || data == NULL) {
00179       std::ostringstream oss;
00180       oss << "stk::linsys::copy_vector_to_mesh ERROR, bad data type, or ";
00181       oss << " field (" << field->name() << ") not found on entity with type " << entity.entity_type();
00182       oss << " and ID " << entity.identifier();
00183       std::string str = oss.str();
00184       throw std::runtime_error(str.c_str());
00185     }
00186 
00187     double * double_data = reinterpret_cast<double *>(data);
00188 
00189     double_data[offset_into_field] = values[i];
00190 
00191   }
00192 }
00193 
00194 }//namespace linsys
00195 }//namespace stk
00196 

Generated on Tue Jul 13 09:27:32 2010 for Sierra Toolkit by  doxygen 1.4.7