Albany: a Trilinos-based PDE code: AAdapt_STKAdapt

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 "AAdapt_STKAdapt.hpp"
00008 #include "Teuchos_TimeMonitor.hpp"
00009 #include "Albany_AbstractSTKFieldContainer.hpp"
00010 
00011 #include <stk_adapt/IEdgeBasedAdapterPredicate.hpp>
00012 #include <stk_adapt/IElementBasedAdapterPredicate.hpp>
00013 #include <stk_adapt/PredicateBasedElementAdapter.hpp>
00014 #include <stk_adapt/PredicateBasedEdgeAdapter.hpp>
00015 
00016 
00017 template<class SizeField>
00018 AAdapt::STKAdapt<SizeField>::
00019 STKAdapt(const Teuchos::RCP<Teuchos::ParameterList>& params_,
00020          const Teuchos::RCP<ParamLib>& paramLib_,
00021          Albany::StateManager& StateMgr_,
00022          const Teuchos::RCP<const Epetra_Comm>& comm_) :
00023   AAdapt::AbstractAdapter(params_, paramLib_, StateMgr_, comm_),
00024   remeshFileIndex(1) {
00025 
00026   disc = StateMgr_.getDiscretization();
00027 
00028   stk_discretization = static_cast<Albany::STKDiscretization*>(disc.get());
00029 
00030   genericMeshStruct = Teuchos::rcp_dynamic_cast<Albany::GenericSTKMeshStruct>(stk_discretization->getSTKMeshStruct());
00031 
00032   eMesh = genericMeshStruct->getPerceptMesh();
00033   TEUCHOS_TEST_FOR_EXCEPT(eMesh.is_null());
00034 
00035   refinerPattern = genericMeshStruct->getRefinerPattern();
00036   TEUCHOS_TEST_FOR_EXCEPT(refinerPattern.is_null());
00037 
00038   num_iterations = adapt_params_->get<int>("Max Number of STK Adapt Iterations", 1);
00039 
00040   // Save the initial output file name
00041   base_exo_filename = stk_discretization->getSTKMeshStruct()->exoOutFile;
00042 
00043 
00044 }
00045 
00046 template<class SizeField>
00047 AAdapt::STKAdapt<SizeField>::
00048 ~STKAdapt() {
00049 }
00050 
00051 template<class SizeField>
00052 bool
00053 AAdapt::STKAdapt<SizeField>::queryAdaptationCriteria() {
00054 
00055   if(adapt_params_->get<std::string>("Remesh Strategy", "None").compare("Continuous") == 0){
00056 
00057     if(iter > 1)
00058 
00059       return true;
00060 
00061     else
00062 
00063       return false;
00064 
00065   }
00066 
00067   Teuchos::Array<int> remesh_iter = adapt_params_->get<Teuchos::Array<int> >("Remesh Step Number");
00068 
00069   for(int i = 0; i < remesh_iter.size(); i++)
00070 
00071     if(iter == remesh_iter[i])
00072 
00073       return true;
00074 
00075   return false;
00076 
00077 }
00078 
00079 template<class SizeField>
00080 void
00081 AAdapt::STKAdapt<SizeField>::printElementData() {
00082 
00083   Albany::StateArrays& sa = disc->getStateArrays();
00084   Albany::StateArrayVec& esa = sa.elemStateArrays;
00085   int numElemWorksets = esa.size();
00086   Teuchos::RCP<Albany::StateInfoStruct> stateInfo = state_mgr_.getStateInfoStruct();
00087 
00088   std::cout << "Num Worksets = " << numElemWorksets << std::endl;
00089 
00090   for(unsigned int i = 0; i < stateInfo->size(); i++) {
00091 
00092     const std::string stateName = (*stateInfo)[i]->name;
00093     const std::string init_type = (*stateInfo)[i]->initType;
00094     std::vector<int> dims;
00095     esa[0][stateName].dimensions(dims);
00096     int size = dims.size();
00097 
00098     std::cout << "Meshadapt: have element field \"" << stateName << "\" of type \"" << init_type << "\"" << std::endl;
00099 
00100     if(init_type == "scalar") {
00101 
00102 
00103       switch(size) {
00104 
00105         case 1:
00106           std::cout << "esa[ws][stateName](0)" << std::endl;
00107           std::cout << "Size = " << dims[0] << std::endl;
00108           break;
00109 
00110         case 2:
00111           std::cout << "esa[ws][stateName](cell, qp)" << std::endl;
00112           std::cout << "Size = " << dims[0] << " , " << dims[1] << std::endl;
00113           break;
00114 
00115         case 3:
00116           std::cout << "esa[ws][stateName](cell, qp, i)" << std::endl;
00117           std::cout << "Size = " << dims[0] << " , " << dims[1] << " , " << dims[2] << std::endl;
00118           break;
00119 
00120         case 4:
00121           std::cout << "esa[ws][stateName](cell, qp, i, j)" << std::endl;
00122           std::cout << "Size = " << dims[0] << " , " << dims[1] << " , " << dims[2] << " , " << dims[3] << std::endl;
00123           break;
00124 
00125       }
00126     }
00127 
00128     else if(init_type == "identity") {
00129       std::cout << "Have an identity matrix: " << "esa[ws][stateName](cell, qp, i, j)" << std::endl;
00130     }
00131   }
00132 }
00133 
00134 template<class SizeField>
00135 bool
00136 AAdapt::STKAdapt<SizeField>::adaptMesh(const Epetra_Vector& sol, const Epetra_Vector& ovlp_sol) {
00137 
00138   *output_stream_ << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
00139   *output_stream_ << "Adapting mesh using AAdapt::STKAdapt method        " << std::endl;
00140   *output_stream_ << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
00141 
00142   Albany::AbstractSTKFieldContainer::IntScalarFieldType* proc_rank_field 
00143       = genericMeshStruct->getFieldContainer()->getProcRankField();
00144   Albany::AbstractSTKFieldContainer::IntScalarFieldType* refine_field 
00145       = genericMeshStruct->getFieldContainer()->getRefineField();
00146 
00147   // Save the current results and close the exodus file
00148 
00149   // Create a remeshed output file naming convention by adding the remesh_file_index_ ahead of the period
00150   std::ostringstream ss;
00151   std::string str = base_exo_filename;
00152   ss << "_" << remeshFileIndex << ".";
00153   str.replace(str.find('.'), 1, ss.str());
00154 
00155   *output_stream_ << "Remeshing: renaming output file to - " << str << std::endl;
00156 
00157   // Open the new exodus file for results
00158   stk_discretization->reNameExodusOutput(str);
00159 
00160   remeshFileIndex++;
00161 
00162 //  printElementData();
00163 
00164   SizeField set_ref_field(*eMesh);
00165   eMesh->elementOpLoop(set_ref_field, refine_field);
00166 
00167   //    SetUnrefineField set_unref_field(*eMesh);
00168   //eMesh.elementOpLoop(set_ref_field, refine_field);
00169 
00170 //  eMesh->save_as("local_tet_N_5_ElementBased_0_.e");
00171 
00172   stk::adapt::ElementRefinePredicate erp(0, refine_field, 0.0);
00173 
00174   stk::adapt::PredicateBasedElementAdapter<stk::adapt::ElementRefinePredicate>
00175   breaker(erp, *eMesh, *refinerPattern, proc_rank_field);
00176 
00177   breaker.setRemoveOldElements(false);
00178   breaker.setAlwaysInitializeNodeRegistry(false);
00179 
00180   for(int ipass = 0; ipass < 3; ipass++) {
00181 
00182     eMesh->elementOpLoop(set_ref_field, refine_field);
00183 
00184 #if 0
00185     std::vector<stk::mesh::Entity*> elems;
00186     const std::vector<stk::mesh::Bucket*>& buckets = eMesh->get_bulk_data()->buckets(eMesh->element_rank());
00187 
00188     for(std::vector<stk::mesh::Bucket*>::const_iterator k = buckets.begin() ; k != buckets.end() ; ++k) {
00189       stk::mesh::Bucket& bucket = **k ;
00190 
00191       const unsigned num_elements_in_bucket = bucket.size();
00192 
00193       for(unsigned i_element = 0; i_element < num_elements_in_bucket; i_element++) {
00194         stk::mesh::Entity& element = bucket[i_element];
00195         double* f_data = stk::percept::PerceptMesh::field_data_entity(refine_field, element);
00196 
00197         std::cout << "Element: " << element.identifier() << "Refine field: " << f_data[0] << std::endl;
00198       }
00199     }
00200 
00201 #endif
00202 
00203 
00204     //     std::cout << "P[" << eMesh->get_rank() << "] ipass= " << ipass << std::endl;
00205     breaker.doBreak();
00206     //     std::cout << "P[" << eMesh->get_rank() << "] done... ipass= " << ipass << std::endl;
00207     //     eMesh->save_as("local_tet_N_5_ElementBased_1_ipass_"+Teuchos::toString(ipass)+"_.e");
00208   }
00209 
00210   breaker.deleteParentElements();
00211 //  eMesh->save_as("local_tet_N_5_ElementBased_1_.e");
00212 
00213   // Throw away all the Albany data structures and re-build them from the mesh
00214 
00215   if(adapt_params_->get<bool>("Rebalance", false))
00216 
00217     genericMeshStruct->rebalanceAdaptedMesh(adapt_params_, epetra_comm_);
00218     
00219   stk_discretization->updateMesh();
00220 //  printElementData();
00221 
00222   return true;
00223 
00224 }
00225 
00227 template<class SizeField>
00228 void
00229 AAdapt::STKAdapt<SizeField>::
00230 solutionTransfer(const Epetra_Vector& oldSolution,
00231                  Epetra_Vector& newSolution) {
00232 #if 0
00233 
00234   // Just copy across for now!
00235 
00236   std::cout << "WARNING: solution transfer not implemented yet!!!" << std::endl;
00237 
00238 
00239   std::cout << "AAdapt<> will now throw an exception from line #156" << std::endl;
00240 
00241   newSolution = oldSolution;
00242 
00243 #endif
00244 
00245 }
00246 
00247 template<class SizeField>
00248 Teuchos::RCP<const Teuchos::ParameterList>
00249 AAdapt::STKAdapt<SizeField>::getValidAdapterParameters() const {
00250   Teuchos::RCP<Teuchos::ParameterList> validPL =
00251     this->getGenericAdapterParams("ValidSTKAdaptParams");
00252 
00253   Teuchos::Array<int> defaultArgs;
00254 
00255   validPL->set<Teuchos::Array<int> >("Remesh Step Number", defaultArgs, "Iteration step at which to remesh the problem");
00256   validPL->set<std::string>("Remesh Strategy", "", "Strategy to use when remeshing: Continuous - remesh every step.");
00257   validPL->set<int>("Max Number of STK Adapt Iterations", 1, "Number of iterations to limit stk_adapt to");
00258   validPL->set<std::string>("Refiner Pattern", "", "Element pattern to use for refinement");
00259   validPL->set<double>("Target Element Size", 0.1, "Seek this element size when isotropically adapting");
00260   validPL->set<bool>("Rebalance", "1", "Rebalance mesh after each refinement operation");
00261 
00262   return validPL;
00263 }
00264 
00265
AAdapt_STKAdapt_Def.hpp
