10 #include <stk_linsys/FeiBaseIncludes.hpp> 11 #include <stk_linsys/LinsysFunctions.hpp> 12 #include <stk_linsys/ImplDetails.hpp> 14 #include <stk_mesh/base/GetBuckets.hpp> 15 #include <stk_mesh/base/FieldData.hpp> 17 #include <fei_trilinos_macros.hpp> 18 #include <fei_Solver_AztecOO.hpp> 19 #include <fei_Trilinos_Helpers.hpp> 27 stk_classic::mesh::EntityRank connected_entity_rank,
33 std::vector<mesh::Bucket*> part_buckets;
36 if (part_buckets.empty())
return;
38 DofMapper& dof_mapper = ls.get_DofMapper();
40 dof_mapper.
add_dof_mappings(mesh_bulk, selector, connected_entity_rank, field);
46 int num_connected = rel.second - rel.first;
48 fei::SharedPtr<fei::MatrixGraph> matgraph = ls.get_fei_MatrixGraph();
50 int pattern_id = matgraph->definePattern(num_connected, connected_entity_rank, field_id);
53 for(
size_t i=0; i<part_buckets.size(); ++i) {
54 num_entities += part_buckets[i]->size();
57 int block_id = matgraph->initConnectivityBlock(num_entities, pattern_id);
59 std::vector<int> connected_ids(num_connected);
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) {
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();
71 int conn_id = entity.identifier();
72 matgraph->initConnectivity(block_id, conn_id, &connected_ids[0]);
82 unsigned field_component,
83 double prescribed_value)
85 std::vector<stk_classic::mesh::Bucket*> buckets;
89 int field_id = ls.get_DofMapper().get_field_id(field);
90 std::vector<int> entity_ids;
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) {
101 std::vector<double> prescribed_values(entity_ids.size(), prescribed_value);
103 ls.get_fei_LinearSystem()->loadEssentialBCs(entity_ids.size(),
106 field_id, field_component,
107 &prescribed_values[0]);
110 int fei_solve(
int & status, fei::LinearSystem &fei_ls,
const Teuchos::ParameterList & params )
113 Solver_AztecOO solver_az;
115 fei::ParameterSet param_set;
117 Trilinos_Helpers::copy_parameterlist(params, param_set);
119 int iterationsTaken = 0;
121 return solver_az.solve( & fei_ls,
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();
136 A->multiply(x.get(), &r);
138 r.update(1, b.get(), -1);
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];
153 MPI_Comm comm = r.getVectorSpace()->getCommunicator();
154 double global_sum = 0;
156 MPI_Allreduce(&local_sum, &global_sum, num_doubles, MPI_DOUBLE, MPI_SUM, comm);
158 double global_sum = local_sum;
160 return std::sqrt(global_sum);
168 vec.scatterToOverlap();
170 std::vector<int> shared_and_owned_indices;
172 vec.getVectorSpace()->getIndices_SharedAndOwned(shared_and_owned_indices);
174 size_t num_values = shared_and_owned_indices.size();
176 if(num_values == 0) {
180 std::vector<double> values(num_values);
181 vec.copyOut(num_values,&shared_and_owned_indices[0],&values[0]);
183 stk_classic::mesh::EntityRank ent_type;
184 stk_classic::mesh::EntityId ent_id;
186 int offset_into_field;
188 for(
size_t i = 0; i < num_values; ++i)
191 dof.
get_dof( shared_and_owned_indices[i],
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());
211 double * double_data =
reinterpret_cast<double *
>(data);
213 double_data[offset_into_field] = values[i];
220 fei::SharedPtr<fei::VectorSpace> vspace = matrix.getMatrixGraph()->getRowSpace();
222 int numRows = vspace->getNumIndices_Owned();
223 std::vector<int> rows(numRows);
224 vspace->getIndices_Owned(numRows, &rows[0], numRows);
226 std::vector<int> indices;
227 std::vector<double> coefs;
229 for(
size_t i=0; i<rows.size(); ++i) {
231 matrix.getRowLength(rows[i], rowlen);
233 if ((
int)indices.size() < rowlen) {
234 indices.resize(rowlen);
235 coefs.resize(rowlen);
238 matrix.copyOutRow(rows[i], rowlen, &coefs[0], &indices[0]);
240 for(
int j=0; j<rowlen; ++j) {
244 double* coefPtr = &coefs[0];
245 matrix.copyIn(1, &rows[i], rowlen, &indices[0], &coefPtr);
250 const fei::Matrix& src_matrix,
251 fei::Matrix& dest_matrix)
253 fei::SharedPtr<fei::VectorSpace> vspace = src_matrix.getMatrixGraph()->getRowSpace();
255 int numRows = vspace->getNumIndices_Owned();
256 std::vector<int> rows(numRows);
257 vspace->getIndices_Owned(numRows, &rows[0], numRows);
259 std::vector<int> indices;
260 std::vector<double> coefs;
262 for(
size_t i=0; i<rows.size(); ++i) {
264 src_matrix.getRowLength(rows[i], rowlen);
266 if ((
int)indices.size() < rowlen) {
267 indices.resize(rowlen);
268 coefs.resize(rowlen);
271 src_matrix.copyOutRow(rows[i], rowlen, &coefs[0], &indices[0]);
273 for(
int j=0; j<rowlen; ++j) {
277 double* coefPtr = &coefs[0];
278 dest_matrix.sumIn(1, &rows[i], rowlen, &indices[0], &coefPtr);
285 fei::SharedPtr<fei::VectorSpace> vspace = vec.getVectorSpace();
287 int numIndices = vspace->getNumIndices_Owned();
288 std::vector<int> indices(numIndices);
289 vspace->getIndices_Owned(numIndices, &indices[0], numIndices);
291 std::vector<double> coefs(numIndices);
293 vec.copyOut(numIndices, &indices[0], &coefs[0]);
295 for(
size_t j=0; j<coefs.size(); ++j) {
299 vec.copyIn(numIndices, &indices[0], &coefs[0]);
303 const fei::Vector& src_vector,
304 fei::Vector& dest_vector)
306 fei::SharedPtr<fei::VectorSpace> vspace = src_vector.getVectorSpace();
308 int numIndices = vspace->getNumIndices_Owned();
309 std::vector<int> indices(numIndices);
310 vspace->getIndices_Owned(numIndices, &indices[0], numIndices);
312 std::vector<double> coefs(numIndices);
314 src_vector.copyOut(numIndices, &indices[0], &coefs[0]);
316 for(
size_t j=0; j<coefs.size(); ++j) {
320 dest_vector.sumIn(numIndices, &indices[0], &coefs[0]);
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)
int fei_solve(int &status, fei::LinearSystem &fei_ls, const Teuchos::ParameterList ¶ms)
Field base class with an anonymous data type and anonymous multi-dimension.
FieldTraits< field_type >::data_type * field_data(const field_type &f, const Bucket::iterator i)
Pointer to the field data array.
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
This is a class for selecting buckets based on a set of meshparts and set logic.
const std::vector< Bucket * > & buckets(EntityRank rank) const
Query all buckets of a given entity rank.
Entity * get_entity(EntityRank entity_rank, EntityId entity_id) const
Get entity with a given key.
An application-defined subset of a problem domain.
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)
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-...
int get_field_id(const stk_classic::mesh::FieldBase &field) const
Manager for an integrated collection of entities, entity relations, and buckets of field data...
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)
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...
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)
void scale_matrix(double scalar, fei::Matrix &matrix)
EntityRank entity_rank(const EntityKey &key)
Given an entity key, return an entity type (rank).