adaptive_pitchfork.cc File Reference

#include "generic.h"
#include "navier_stokes.h"
#include "meshes/rectangular_quadmesh.h"

Classes
class	CompareNodes
class	RefineablePorousChannelProblem< ELEMENT >
	Porous channel flow on a refineable mesh. More...

Namespaces
	Global_Physical_Variables
	Global variables.

Functions
int	main ()
	Driver for RefineablePorousChannel test problem using a. More...

Variables
double	Global_Physical_Variables::L = 8.0
	Length of beam. More...

Function Documentation

main()

int main

(

)

Driver for RefineablePorousChannel test problem using a.

///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////

{
 // Build the problem with QTaylorHoodElements
 RefineablePorousChannelProblem<RefineableQTaylorHoodElement<2> > 
  problem;
 
 problem.set_boundary_conditions();
 problem.steady_newton_solve();
 
 
 //Assign memory for the eigenvalues and eigenvectors
 Vector<std::complex<double> > eigenvalues;
 Vector<DoubleVector> eigenvectors;
 //Set the eigen solver to the LAPACK version which is
 //in every distributioin of oomph-lib
 problem.eigen_solver_pt() = new LAPACK_QZ;
 //Solve the eigenproblem
 problem.solve_eigenproblem(4,eigenvalues,eigenvectors);  
 
 //Find the eigenvalue with greatest real part
 unsigned ev_crit_index=0;  double ev_crit_value=0.0;
 
 //Loop over the values until we get a finite one
 for(unsigned ev=0;ev<eigenvalues.size();ev++)
  {
   double ev_crit_value = eigenvalues[ev].real();
   if(std::isfinite(ev_crit_value)) {ev_crit_index = ev; break;}
  }
 
 //Now loop over the other values and find the eigenvalue with 
 //greatest real part
 for(unsigned ev=ev_crit_index+1;ev<eigenvalues.size();ev++)
  {
   //Get the current value
   double real_part = eigenvalues[ev].real();
   //If it's finite do the comparison
   if(std::isfinite(real_part))
    {
     if(real_part > ev_crit_value)
      {
       ev_crit_value = real_part;
       ev_crit_index = ev;
      }
    }
  }
 
 problem.activate_pitchfork_tracking(&Global_Physical_Variables::Re,
                                     eigenvectors[ev_crit_index]);
 
 //Solve with two rounds of adaptivity
 problem.steady_newton_solve(2);
 
 //Output the solution at the bifurcation point
 problem.mesh_pt()->output("bif_soln.dat",5);
 
 //Report the final answer
 std::cout << "Pitchfork found at " << Global_Physical_Variables::Re  
           << "\n";
 std::cout << "The slack parameter is " << problem.dof(problem.ndof()-1) 
           << "\n";
 
 //Write the final answer into a file
 std::ofstream trace("trace.dat");
 trace << Global_Physical_Variables::Re << " "
       << problem.dof(problem.ndof()-1) << "\n";
 trace.close();
 
 //Deactivate the bifurcation tracking (for neatness)
 problem.deactivate_bifurcation_tracking();
 
 //Delete the eigensolver that we allocated
 delete problem.eigen_solver_pt();
 
} // end_of_main

References isfinite, problem, and oomph::Global_Physical_Variables::Re.