bifurcation_tracking/pitchfork.cc File Reference

#include <iostream>
#include <fstream>
#include <cmath>
#include <typeinfo>
#include <algorithm>
#include <cstdio>
#include "generic.h"

Classes
class	Mesh1D< ELEMENT >
class	SSPorousChannelEquations
class	SSPorousChannelElement
class	UniformTranspiration< ELEMENT >

Namespaces
	Global_Physical_Variables
	Global variables.

Functions
int	main ()

Function Documentation

main()

int main

(

)

{
 //Load the namespace
 using namespace Global_Physical_Variables;
 
 //Set the number of elements in each section of the mesh.
 unsigned Nx1=50, Nx2=50;
 //Set up the problem
 UniformTranspiration<SSPorousChannelElement> problem(Nx1,Nx2);
 
 //Open output file
 ofstream trace("trace.dat");
 
 //Track the node in the middle
 //The element Nx1 contains this node
 SSPorousChannelElement *Test_pt = 
  dynamic_cast<SSPorousChannelElement *>(problem.mesh_pt()->element_pt(Nx1));
  //Set the position of the node within the element
 Vector<double> s(1); s[0] = 1.0;
 
 //Step up in Reynolds number
 for(unsigned i=0;i<6;i++)
  {
   //Increase the Reynolds number by 0.5 each time
   Re = i*0.5;
   problem.steady_newton_solve();
   //Output to the trace file
   trace << Re << " " 
        << Test_pt->interpolated_f(s,1) << " " 
        << Test_pt->interpolated_f(s,0) << std::endl;
  }
 
 //Specify the symmetry the hard way
  unsigned n_dof = problem.ndof();
  Vector<double> backup(n_dof);
  for(unsigned n=0;n<n_dof;n++)
   {
    backup[n] = problem.dof(n);
   }
  //Now sort out the problem
  unsigned n_node=problem.mesh_pt()->nnode();
  for(unsigned n=0;n<n_node;n++)
   {
    Node* nod_pt = problem.mesh_pt()->node_pt(n);
    double x = nod_pt->x(0);
    if(!nod_pt->is_pinned(1))
     {
      nod_pt->set_value(1,x*x);
     }
    if(!nod_pt->is_pinned(0))
     {
      nod_pt->set_value(0,0.0);
     }
   }
 
  LinearAlgebraDistribution dist(problem.communicator_pt(),n_dof,false);
  DoubleVector symm(&dist,0.0);
  for(unsigned n=0;n<n_dof;n++)
   {
    symm[n] = problem.dof(n);
    problem.dof(n) = backup[n];
   }
   
 
 //Let's try to find the pitchfork it
 problem.activate_pitchfork_tracking(&Re,symm);
 
 problem.steady_newton_solve();
 
 std::cout << "Pitchfork detected at " << Re  << std::endl;
 std::cout << "The slack parameter is " << problem.dof(2*n_dof+1) << std::endl;
 
 //Output to the trace file
 trace << Re << " " 
       << Test_pt->interpolated_f(s,1) << " " 
       << Test_pt->interpolated_f(s,0) << std::endl;
 
 
 }

References i, SSPorousChannelEquations::interpolated_f(), oomph::Data::is_pinned(), n, problem, GlobalPhysicalVariables::Re, s, oomph::Data::set_value(), symm(), plotDoE::x, and oomph::Node::x().