mpi/distribution/fish_poisson/fish_poisson.cc File Reference

#include "generic.h"
#include "poisson.h"
#include "meshes/fish_mesh.h"

Classes
class	RefineableFishPoissonProblem< ELEMENT >

Namespaces
	ConstSourceForPoisson
	Namespace for const source term in Poisson equation.

Functions
void	ConstSourceForPoisson::source_function (const Vector< double > &x, double &source)
	Const source function. More...
void	solve_with_incremental_adaptation ()
void	solve_with_fully_automatic_adaptation ()
void	solve_with_selected_refinement_pattern ()
	Solve 2D fish Poisson problem with a selected refinement pattern. More...
int	main (int argc, char **argv)

Function Documentation

main()

int main	(	int argc	,
		char **	argv
	)

Demonstrate how to solve 2D Poisson problem in fish-shaped domain with mesh adaptation.

{
 // initialise MPI
#ifdef OOMPH_HAS_MPI
 MPI_Helpers::init(argc,argv);
#endif
 
 // Solve using a pre-selected refinement pattern
 solve_with_selected_refinement_pattern();
 
 // Solve with adaptation, docing the intermediate steps
 solve_with_incremental_adaptation();
 
 // Solve directly, with fully automatic adaptation
 solve_with_fully_automatic_adaptation();
 
#ifdef OOMPH_HAS_MPI
 MPI_Helpers::finalize();
#endif
 
} // end of main

References solve_with_fully_automatic_adaptation(), solve_with_incremental_adaptation(), and solve_with_selected_refinement_pattern().

solve_with_fully_automatic_adaptation()

void solve_with_fully_automatic_adaptation

(

)

Demonstrate how to solve 2D Poisson problem in fish-shaped domain with fully automatic mesh adaptation

{
 //Set up the problem with nine-node refineable Poisson elements
 RefineableFishPoissonProblem<RefineableQPoissonElement<2,3> > problem;
  
 // Setup labels for output
 //------------------------
 DocInfo doc_info, mesh_doc_info;
 bool report_stats=true;
  
 // Set output directory
 doc_info.set_directory("RESLT_fully_automatic"); 
 mesh_doc_info.set_directory("RESLT_adapt_mesh"); 
  
 // Step number
 doc_info.number()=0;
 mesh_doc_info.number()=20;
 
 // Doc (default) refinement targets
 //----------------------------------
 problem.mesh_pt()->doc_adaptivity_targets(cout);
 
 // Solve/doc the problem with fully automatic adaptation
 //------------------------------------------------------
 
 // Refine coarse original mesh first
 //----------------------------------
 problem.refine_uniformly();
 problem.refine_uniformly();
 
 oomph_info << "-----------------------------------------" << std::endl;
 oomph_info << "--- Distributing problem (fully auto) ---" << std::endl;
 oomph_info << "-----------------------------------------" << std::endl;
 
 std::ifstream input_file;
 std::ofstream output_file;
 char filename[100];
 
 // Get the partition to be used from file
 unsigned n_partition=problem.mesh_pt()->nelement();
 Vector<unsigned> element_partition(n_partition,0);
 sprintf(filename,"fish_fully_automatic_partition.dat");
 input_file.open(filename);
 std::string input_string;
 for (unsigned e=0;e<n_partition;e++)
  {
   getline(input_file,input_string,'\n');
   element_partition[e]=atoi(input_string.c_str());
  }
 
 // Distribute and check halo schemes
 problem.distribute(element_partition,mesh_doc_info,report_stats);
 problem.check_halo_schemes(mesh_doc_info);
 
 //Maximum number of adaptations:
 unsigned max_adapt=5;
 
 oomph_info << "Solve with max_adapt=" << max_adapt << std::endl;
 
 //Solve the problem; perform up to specified number of adaptations.
 problem.newton_solve(max_adapt);
 
 //Output solution
 oomph_info << "-----------------------" << std::endl;
 oomph_info << "Now output the solution" << std::endl;
 oomph_info << "-----------------------" << std::endl;
 problem.doc_solution(doc_info);   
 
 //increment doc_info number
 doc_info.number()++;
 
} // end black box

References e(), MergeRestartFiles::filename, oomph::DocInfo::number(), oomph::oomph_info, problem, oomph::DocInfo::set_directory(), and oomph::Global_string_for_annotation::string().

Referenced by main().

solve_with_incremental_adaptation()

void solve_with_incremental_adaptation

(

)

Demonstrate how to solve 2D Poisson problem in fish-shaped domain with mesh adaptation. First we solve on the original coarse mesh. Next we do a few uniform refinement steps and re-solve. Finally, we enter into an automatic adapation loop.

{
 
 oomph_info << "Min/max about to start incremental\n";
 
 //Set up the problem with nine-node refineable Poisson elements
 RefineableFishPoissonProblem<RefineableQPoissonElement<2,3> > problem;
 
 // Setup labels for output
 //------------------------
 DocInfo doc_info, mesh_doc_info;
 bool report_stats=true;
 
 // Set output directory
 doc_info.set_directory("RESLT_incremental2"); 
 mesh_doc_info.set_directory("RESLT_incr_mesh"); 
 
 // Step number
 doc_info.number()=1;
 mesh_doc_info.number()=10;
  
 // Increase minimum permitted error from default (10e-5) to 10e-4
 problem.mesh_pt()->min_permitted_error()=1.0e-4;
  
 // Doc (default) refinement targets
 //----------------------------------
 problem.mesh_pt()->doc_adaptivity_targets(cout);
 
 // solve
 problem.newton_solve();
  
 //Output solution
 doc_info.label()="initial";
 problem.doc_solution(doc_info);
 
 //Increment counter for solutions 
 doc_info.number()++;
 
 //Loop for refining uniformly
 for (unsigned mym=0;mym<2;mym++)
  {
   //On the first step refine uniformly twice and then distribute
   if (mym==0) 
    {
     problem.refine_uniformly();
     problem.refine_uniformly();
 
     oomph_info << "Distributing problem" << std::endl;
     std::ifstream input_file;
     std::ofstream output_file;
     char filename[100];
 
     // Get the partition to be used from file
     unsigned n_partition=problem.mesh_pt()->nelement();
     Vector<unsigned> element_partition(n_partition,0);
     sprintf(filename,"fish_incremental_partition.dat");
     input_file.open(filename);
     std::string input_string;
     for (unsigned e=0;e<n_partition;e++)
      {
       getline(input_file,input_string,'\n');
       element_partition[e]=atoi(input_string.c_str());
      }
 
     // Distribute and check halo schemes
     problem.distribute(element_partition,mesh_doc_info,report_stats);
     problem.check_halo_schemes(mesh_doc_info);
    }
  
   // another round of refinement
   oomph_info << "Refine again..." << std::endl;
   problem.refine_uniformly();
 
   // Solve the problem 
   oomph_info << "Solve again..." << std::endl;
 
   problem.newton_solve();
 
   //Output solution
   doc_info.label()="after first solve after distribution and uniform refinement";
   problem.doc_solution(doc_info);
 
   //Increment counter for solutions 
   doc_info.number()++;
  }
 
 oomph_info << "Automatic adaptation starting" << std::endl << std::endl;
 
 // Now do (up to) four rounds of fully automatic adapation in response to 
 //-----------------------------------------------------------------------
 // error estimate
 //---------------
 unsigned max_solve=4;
 for (unsigned isolve=0;isolve<max_solve;isolve++)
  {
   oomph_info << "Adapting problem - isolve=" << isolve << std::endl;
 
   // Adapt problem/mesh
   problem.adapt(); 
 
   // Re-solve the problem if the adaptation has changed anything
#ifdef OOMPH_HAS_MPI
   // Make sure all processors know if refinement is taking place
   // Adaptation only converges if ALL the processes have no
   // refinement or unrefinement to perform
   unsigned total_refined;
   unsigned total_unrefined;
   unsigned n_refined=problem.mesh_pt()->nrefined();
   unsigned n_unrefined=problem.mesh_pt()->nunrefined();
 
   MPI_Allreduce(&n_refined,&total_refined,1,MPI_INT,MPI_SUM,
                 problem.communicator_pt()->mpi_comm());
   n_refined=total_refined;
 
   MPI_Allreduce(&n_unrefined,&total_unrefined,1,MPI_INT,MPI_SUM,
                 problem.communicator_pt()->mpi_comm());
   n_unrefined=total_unrefined;
#endif
 
   oomph_info << "---> " << n_refined << " elements to be refined, and " 
              << n_unrefined << " to be unrefined, in total." << std::endl;
 
   if ((n_refined!=0)||(n_unrefined!=0))
    {
     problem.newton_solve();
    }
   else
    {
     cout << "Mesh wasn't adapted --> we'll stop here" << std::endl;
     break;
    }
   
   //Output solution
   doc_info.label()="solve after adapt";
   problem.doc_solution(doc_info);
   
   //Increment counter for solutions 
   doc_info.number()++;
  }
 
 // Loop for uniform unrefinement
 for(unsigned myn=0;myn<10;myn++)
  {
   // Test the unrefine_uniformly command and re-solve
   //-------------------------------------------------
   problem.unrefine_uniformly();
 
   // Solve the problem 
   problem.newton_solve();
 
   //Output solution
   doc_info.label()="solve after uniform unrefinement";
   problem.doc_solution(doc_info);
 
   //Increment counter for solutions 
   doc_info.number()++;
  }
 
 oomph_info << "Min/max end incremental\n";
 
} // end of incremental

References e(), MergeRestartFiles::filename, oomph::DocInfo::label(), oomph::DocInfo::number(), oomph::oomph_info, problem, oomph::DocInfo::set_directory(), and oomph::Global_string_for_annotation::string().

Referenced by main().

solve_with_selected_refinement_pattern()

void solve_with_selected_refinement_pattern

(

)

Solve 2D fish Poisson problem with a selected refinement pattern.

{
 //Set up the problem with nine-node refineable Poisson elements
 RefineableFishPoissonProblem<RefineableQPoissonElement<2,3> > problem;
  
 // Setup labels for output
 //------------------------
 DocInfo doc_info, mesh_doc_info;
 bool report_stats=true;
  
 // Set output directory
 doc_info.set_directory("RESLT_select_refine"); 
 mesh_doc_info.set_directory("RESLT_select_mesh"); 
  
 // Step number
 doc_info.number()=0;
 mesh_doc_info.number()=20;
 
 // Doc (default) refinement targets
 //----------------------------------
 problem.mesh_pt()->doc_adaptivity_targets(cout);
 
 // Refine coarse original mesh first
 //----------------------------------
 problem.refine_uniformly();
 problem.refine_uniformly();
 
 // Distribute the problem and doc mesh info
 oomph_info << "----------------------------------------" << std::endl;
 oomph_info << "--- Distributing problem (selective) ---" << std::endl;
 oomph_info << "----------------------------------------" << std::endl;
 
 std::ifstream input_file;
 std::ofstream output_file;
 char filename[100];
 
 // Get the partition to be used from file
 unsigned n_partition=problem.mesh_pt()->nelement();
 Vector<unsigned> element_partition(n_partition,0);
 sprintf(filename,"fish_selective_partition.dat");
 input_file.open(filename);
 std::string input_string;
 for (unsigned e=0;e<n_partition;e++)
  {
   getline(input_file,input_string,'\n');
   element_partition[e]=atoi(input_string.c_str());
  }
 
 // Distribute and check halo schemes
 problem.distribute(element_partition,mesh_doc_info,report_stats);
 problem.check_halo_schemes(mesh_doc_info);
 
 // Add some selective refinement on top of this
 unsigned nsoln=3; // number of selective refinements - validate, 3
 
 for (unsigned isoln=0; isoln<nsoln; isoln++)
  {
   // Refine selected elements based upon y-coordinate
   // of the central node of the element (node 4)
   Vector<unsigned> more_middle_els;
   unsigned nels=problem.mesh_pt()->nelement();
   for (unsigned e=0; e<nels; e++)
    {
     FiniteElement* el_pt=dynamic_cast<FiniteElement*>
      (problem.mesh_pt()->element_pt(e));
     if ((el_pt->node_pt(4)->x(1)>=(-0.35)) &&
         (el_pt->node_pt(4)->x(1)<=0.35))
      {
       more_middle_els.push_back(e);
      }
    }
 
   oomph_info << "--------------------------------" << std::endl;
   oomph_info << "This time, refine " << more_middle_els.size() 
              << " selected elements" << std::endl;
   oomph_info << "--------------------------------" << std::endl;
   problem.refine_selected_elements(more_middle_els);
 
   // Now solve the problem
   problem.newton_solve();
 
   // Doc the solution
   problem.doc_solution(doc_info);   
 
   // Increment doc_info number
   doc_info.number()++;
  }
 
}

References e(), MergeRestartFiles::filename, oomph::FiniteElement::node_pt(), oomph::DocInfo::number(), oomph::oomph_info, problem, oomph::DocInfo::set_directory(), oomph::Global_string_for_annotation::string(), and oomph::Node::x().

Referenced by main().