Albany: a Trilinos-based PDE code: Albany_MultiSTKFieldContainer

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 
00007 #include <iostream>
00008 #include <string>
00009 
00010 #include "Albany_MultiSTKFieldContainer.hpp"
00011 
00012 #include "Albany_Utils.hpp"
00013 #include <stk_mesh/base/GetBuckets.hpp>
00014 
00015 #ifdef ALBANY_SEACAS
00016 #include <stk_io/IossBridge.hpp>
00017 #endif
00018 
00019 template<bool Interleaved>
00020 Albany::MultiSTKFieldContainer<Interleaved>::MultiSTKFieldContainer(
00021   const Teuchos::RCP<Teuchos::ParameterList>& params_,
00022   stk::mesh::fem::FEMMetaData* metaData_,
00023   stk::mesh::BulkData* bulkData_,
00024   const int neq_,
00025   const AbstractFieldContainer::FieldContainerRequirements& req,
00026   const int numDim_,
00027   const Teuchos::RCP<Albany::StateInfoStruct>& sis,
00028   const Teuchos::Array<std::string>& solution_vector,
00029   const Teuchos::Array<std::string>& residual_vector)
00030   : GenericSTKFieldContainer<Interleaved>(params_, metaData_, bulkData_, neq_, numDim_),
00031     haveResidual(false) {
00032 
00033   typedef typename AbstractSTKFieldContainer::VectorFieldType VFT;
00034   typedef typename AbstractSTKFieldContainer::ScalarFieldType SFT;
00035 
00036   // Check the input
00037 
00038   if(solution_vector.size() == 0) { // Do the default solution vector
00039 
00040     std::string name = params_->get<std::string>("Exodus Solution Name", "solution");
00041     VFT* solution = & metaData_->declare_field< VFT >(name);
00042     stk::mesh::put_field(*solution, metaData_->node_rank() , metaData_->universal_part(), neq_);
00043 #ifdef ALBANY_SEACAS
00044     stk::io::set_field_role(*solution, Ioss::Field::TRANSIENT);
00045 #endif
00046 
00047     sol_vector_name.push_back(name);
00048     sol_index.push_back(this->neq);
00049 
00050   }
00051 
00052   else if(solution_vector.size() == 1) { // User is just renaming the entire solution vector
00053 
00054     VFT* solution = & metaData_->declare_field< VFT >(solution_vector[0]);
00055     stk::mesh::put_field(*solution, metaData_->node_rank() , metaData_->universal_part(), neq_);
00056 #ifdef ALBANY_SEACAS
00057     stk::io::set_field_role(*solution, Ioss::Field::TRANSIENT);
00058 #endif
00059 
00060     sol_vector_name.push_back(solution_vector[0]);
00061     sol_index.push_back(neq_);
00062 
00063   }
00064 
00065   else { // user is breaking up the solution into multiple fields
00066 
00067     // make sure the number of entries is even
00068 
00069     TEUCHOS_TEST_FOR_EXCEPTION((solution_vector.size() % 2), std::logic_error,
00070                                "Error in input file: specification of solution vector layout is incorrect." << std::endl);
00071 
00072     int len, accum = 0;
00073 
00074     for(int i = 0; i < solution_vector.size(); i += 2) {
00075 
00076       if(solution_vector[i + 1] == "V") {
00077 
00078         len = numDim_; // vector
00079         accum += len;
00080         VFT* solution = & metaData_->declare_field< VFT >(solution_vector[i]);
00081         stk::mesh::put_field(*solution, metaData_->node_rank() , metaData_->universal_part(), len);
00082 #ifdef ALBANY_SEACAS
00083         stk::io::set_field_role(*solution, Ioss::Field::TRANSIENT);
00084 #endif
00085         sol_vector_name.push_back(solution_vector[i]);
00086         sol_index.push_back(len);
00087 
00088       }
00089 
00090       else if(solution_vector[i + 1] == "S") {
00091 
00092         len = 1; // scalar
00093         accum += len;
00094         SFT* solution = & metaData_->declare_field< SFT >(solution_vector[i]);
00095         stk::mesh::put_field(*solution, metaData_->node_rank() , metaData_->universal_part());
00096 #ifdef ALBANY_SEACAS
00097         stk::io::set_field_role(*solution, Ioss::Field::TRANSIENT);
00098 #endif
00099         sol_vector_name.push_back(solution_vector[i]);
00100         sol_index.push_back(len);
00101 
00102       }
00103 
00104       else
00105 
00106         TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00107                                    "Error in input file: specification of solution vector layout is incorrect." << std::endl);
00108 
00109     }
00110 
00111     TEUCHOS_TEST_FOR_EXCEPTION(accum != neq_, std::logic_error,
00112                                "Error in input file: specification of solution vector layout is incorrect." << std::endl);
00113 
00114   }
00115 
00116 #ifdef ALBANY_LCM
00117   // do the residual next
00118 
00119   if(residual_vector.size() == 0) { // Do the default residual vector
00120 
00121     std::string name = params_->get<std::string>("Exodus Residual Name", "residual");
00122     VFT* residual = & metaData_->declare_field< VFT >(name);
00123     stk::mesh::put_field(*residual, metaData_->node_rank() , metaData_->universal_part(), neq_);
00124 #ifdef ALBANY_SEACAS
00125     stk::io::set_field_role(*residual, Ioss::Field::TRANSIENT);
00126 #endif
00127 
00128     res_vector_name.push_back(name);
00129     res_index.push_back(neq_);
00130 
00131   }
00132 
00133   else if(residual_vector.size() == 1) { // User is just renaming the entire residual vector
00134 
00135     VFT* residual = & metaData_->declare_field< VFT >(residual_vector[0]);
00136     stk::mesh::put_field(*residual, metaData_->node_rank() , metaData_->universal_part(), neq_);
00137 #ifdef ALBANY_SEACAS
00138     stk::io::set_field_role(*residual, Ioss::Field::TRANSIENT);
00139 #endif
00140 
00141     res_vector_name.push_back(residual_vector[0]);
00142     res_index.push_back(neq_);
00143 
00144   }
00145 
00146   else { // user is breaking up the residual into multiple fields
00147 
00148     // make sure the number of entries is even
00149 
00150     TEUCHOS_TEST_FOR_EXCEPTION((residual_vector.size() % 2), std::logic_error,
00151                                "Error in input file: specification of residual vector layout is incorrect." << std::endl);
00152 
00153     int len, accum = 0;
00154 
00155     for(int i = 0; i < residual_vector.size(); i += 2) {
00156 
00157       if(residual_vector[i + 1] == "V") {
00158 
00159         len = numDim_; // vector
00160         accum += len;
00161         VFT* residual = & metaData_->declare_field< VFT >(residual_vector[i]);
00162         stk::mesh::put_field(*residual, metaData_->node_rank() , metaData_->universal_part(), len);
00163 #ifdef ALBANY_SEACAS
00164         stk::io::set_field_role(*residual, Ioss::Field::TRANSIENT);
00165 #endif
00166         res_vector_name.push_back(residual_vector[i]);
00167         res_index.push_back(len);
00168 
00169       }
00170 
00171       else if(residual_vector[i + 1] == "S") {
00172 
00173         len = 1; // scalar
00174         accum += len;
00175         SFT* residual = & metaData_->declare_field< SFT >(residual_vector[i]);
00176         stk::mesh::put_field(*residual, metaData_->node_rank() , metaData_->universal_part());
00177 #ifdef ALBANY_SEACAS
00178         stk::io::set_field_role(*residual, Ioss::Field::TRANSIENT);
00179 #endif
00180         res_vector_name.push_back(residual_vector[i]);
00181         res_index.push_back(len);
00182 
00183       }
00184 
00185       else
00186 
00187         TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00188                                    "Error in input file: specification of residual vector layout is incorrect." << std::endl);
00189 
00190     }
00191 
00192     TEUCHOS_TEST_FOR_EXCEPTION(accum != neq_, std::logic_error,
00193                                "Error in input file: specification of residual vector layout is incorrect." << std::endl);
00194 
00195   }
00196 
00197   haveResidual = true;
00198 
00199 #endif
00200 
00201   //Do the coordinates
00202   this->coordinates_field = & metaData_->declare_field< VFT >("coordinates");
00203   stk::mesh::put_field(*this->coordinates_field , metaData_->node_rank() , metaData_->universal_part(), numDim_);
00204 #ifdef ALBANY_SEACAS
00205   stk::io::set_field_role(*this->coordinates_field, Ioss::Field::MESH);
00206 #endif
00207 
00208   this->buildStateStructs(sis);
00209 
00210   initializeSTKAdaptation();
00211 
00212 }
00213 
00214 template<bool Interleaved>
00215 Albany::MultiSTKFieldContainer<Interleaved>::~MultiSTKFieldContainer() {
00216 }
00217 
00218 template<bool Interleaved>
00219 void Albany::MultiSTKFieldContainer<Interleaved>::initializeSTKAdaptation() {
00220 
00221   typedef typename AbstractSTKFieldContainer::IntScalarFieldType ISFT;
00222 
00223   this->proc_rank_field =
00224       & this->metaData->template declare_field< ISFT >("proc_rank");
00225 
00226   this->refine_field =
00227       & this->metaData->template declare_field< ISFT >("refine_field");
00228 
00229   // Processor rank field, a scalar
00230   stk::mesh::put_field(
00231       *this->proc_rank_field,
00232       this->metaData->element_rank(),
00233       this->metaData->universal_part());
00234 
00235   stk::mesh::put_field(
00236       *this->refine_field,
00237       this->metaData->element_rank(),
00238       this->metaData->universal_part());
00239 
00240 #ifdef ALBANY_LCM
00241   // Fracture state used for adaptive insertion.
00242   // It exists for all entities except cells (elements).
00243   this->fracture_state =
00244       & this->metaData->template declare_field< ISFT >("fracture_state");
00245 
00246   stk::mesh::EntityRank const
00247   cell_rank = this->metaData->element_rank();
00248 
00249   for (stk::mesh::EntityRank rank = 0; rank < cell_rank; ++rank) {
00250     stk::mesh::put_field(
00251         *this->fracture_state,
00252         rank,
00253         this->metaData->universal_part());
00254 
00255   }
00256 #endif // ALBANY_LCM
00257 
00258 
00259 #ifdef ALBANY_SEACAS
00260   stk::io::set_field_role(*this->proc_rank_field, Ioss::Field::MESH);
00261   stk::io::set_field_role(*this->refine_field, Ioss::Field::MESH);
00262 #ifdef ALBANY_LCM
00263   stk::io::set_field_role(*this->fracture_state, Ioss::Field::MESH);
00264 #endif // ALBANY_LCM
00265 #endif
00266 
00267 }
00268 
00269 template<bool Interleaved>
00270 void Albany::MultiSTKFieldContainer<Interleaved>::fillSolnVector(Epetra_Vector& soln,
00271     stk::mesh::Selector& sel, const Teuchos::RCP<Epetra_Map>& node_map) {
00272 
00273   typedef typename AbstractSTKFieldContainer::VectorFieldType VFT;
00274   typedef typename AbstractSTKFieldContainer::ScalarFieldType SFT;
00275 
00276   // Iterate over the on-processor nodes by getting node buckets and iterating over each bucket.
00277   stk::mesh::BucketVector all_elements;
00278   stk::mesh::get_buckets(sel, this->bulkData->buckets(this->metaData->node_rank()), all_elements);
00279   this->numNodes = node_map->NumMyElements(); // Needed for the getDOF function to work correctly
00280   // This is either numOwnedNodes or numOverlapNodes, depending on
00281   // which map is passed in
00282 
00283   for(stk::mesh::BucketVector::const_iterator it = all_elements.begin() ; it != all_elements.end() ; ++it) {
00284 
00285     const stk::mesh::Bucket& bucket = **it;
00286 
00287     int offset = 0;
00288 
00289     for(int k = 0; k < sol_index.size(); k++) {
00290 
00291       if(sol_index[k] == 1) { // Scalar
00292 
00293         SFT* field = this->metaData->template get_field<SFT>(sol_vector_name[k]);
00294         this->fillVectorHelper(soln, field, node_map, bucket, offset);
00295 
00296       }
00297 
00298       else {
00299 
00300         VFT* field = this->metaData->template get_field<VFT>(sol_vector_name[k]);
00301         this->fillVectorHelper(soln, field, node_map, bucket, offset);
00302 
00303       }
00304 
00305       offset += sol_index[k];
00306 
00307     }
00308 
00309   }
00310 }
00311 
00312 template<bool Interleaved>
00313 void Albany::MultiSTKFieldContainer<Interleaved>::saveSolnVector(const Epetra_Vector& soln,
00314     stk::mesh::Selector& sel, const Teuchos::RCP<Epetra_Map>& node_map) {
00315 
00316   typedef typename AbstractSTKFieldContainer::VectorFieldType VFT;
00317   typedef typename AbstractSTKFieldContainer::ScalarFieldType SFT;
00318 
00319   // Iterate over the on-processor nodes by getting node buckets and iterating over each bucket.
00320   stk::mesh::BucketVector all_elements;
00321   stk::mesh::get_buckets(sel, this->bulkData->buckets(this->metaData->node_rank()), all_elements);
00322   this->numNodes = node_map->NumMyElements(); // Needed for the getDOF function to work correctly
00323   // This is either numOwnedNodes or numOverlapNodes, depending on
00324   // which map is passed in
00325 
00326   for(stk::mesh::BucketVector::const_iterator it = all_elements.begin() ; it != all_elements.end() ; ++it) {
00327 
00328     const stk::mesh::Bucket& bucket = **it;
00329 
00330     int offset = 0;
00331 
00332     for(int k = 0; k < sol_index.size(); k++) {
00333 
00334       if(sol_index[k] == 1) { // Scalar
00335 
00336         SFT* field = this->metaData->template get_field<SFT>(sol_vector_name[k]);
00337         this->saveVectorHelper(soln, field, node_map, bucket, offset);
00338 
00339       }
00340 
00341       else {
00342 
00343         VFT* field = this->metaData->template get_field<VFT>(sol_vector_name[k]);
00344         this->saveVectorHelper(soln, field, node_map, bucket, offset);
00345 
00346       }
00347 
00348       offset += sol_index[k];
00349 
00350     }
00351 
00352   }
00353 }
00354 
00355 template<bool Interleaved>
00356 void Albany::MultiSTKFieldContainer<Interleaved>::saveResVector(const Epetra_Vector& res,
00357     stk::mesh::Selector& sel, const Teuchos::RCP<Epetra_Map>& node_map) {
00358 
00359   typedef typename AbstractSTKFieldContainer::VectorFieldType VFT;
00360   typedef typename AbstractSTKFieldContainer::ScalarFieldType SFT;
00361 
00362   // Iterate over the on-processor nodes by getting node buckets and iterating over each bucket.
00363   stk::mesh::BucketVector all_elements;
00364   stk::mesh::get_buckets(sel, this->bulkData->buckets(this->metaData->node_rank()), all_elements);
00365   this->numNodes = node_map->NumMyElements(); // Needed for the getDOF function to work correctly
00366   // This is either numOwnedNodes or numOverlapNodes, depending on
00367   // which map is passed in
00368 
00369   for(stk::mesh::BucketVector::const_iterator it = all_elements.begin() ; it != all_elements.end() ; ++it) {
00370 
00371     const stk::mesh::Bucket& bucket = **it;
00372 
00373     int offset = 0;
00374 
00375     for(int k = 0; k < res_index.size(); k++) {
00376 
00377       if(res_index[k] == 1) { // Scalar
00378 
00379         SFT* field = this->metaData->template get_field<SFT>(res_vector_name[k]);
00380         this->saveVectorHelper(res, field, node_map, bucket, offset);
00381 
00382       }
00383 
00384       else {
00385 
00386         VFT* field = this->metaData->template get_field<VFT>(res_vector_name[k]);
00387         this->saveVectorHelper(res, field, node_map, bucket, offset);
00388 
00389       }
00390 
00391       offset += res_index[k];
00392 
00393     }
00394 
00395   }
00396 }
00397 
00398 template<bool Interleaved>
00399 void Albany::MultiSTKFieldContainer<Interleaved>::transferSolutionToCoords() {
00400 
00401   const bool MultiSTKFieldContainer_transferSolutionToCoords_not_implemented = true;
00402   TEUCHOS_TEST_FOR_EXCEPT(MultiSTKFieldContainer_transferSolutionToCoords_not_implemented);
00403   //     this->copySTKField(solution_field, this->coordinates_field);
00404 
00405 }
Albany_MultiSTKFieldContainer_Def.hpp
