GlobalFct Namespace Reference

Functions
void	reset (Vector< double > &unknowns)
	Reset unknowns to some nontrivial garbage. More...
void	get_residuals (const Vector< double > &param, const Vector< double > &unknowns, Vector< double > &residuals)
	Global residual fct. More...
void	set_spurious_solution (Vector< double > &unknowns)
	Set unknowns to spurious solution. More...
void	plot_it (const std::string filename)
	Plot "landscape" of residuals (only for 2D problems!) More...

Variables
double	Force =0.0
	Force. More...
double	Stiffness =1.0
	Spring stiffness. More...
double	U_max =5.5
	Max. displacement. More...
unsigned	N_primary =1

Function Documentation

get_residuals()

void GlobalFct::get_residuals	(	const Vector< double > &	param,
		const Vector< double > &	unknowns,
		Vector< double > &	residuals
	)

Global residual fct.

 {
  // Loop over primary unknowns
  for (unsigned j=0;j<N_primary;j++)
   {    
    
    // Mash 'em up
    double alpha=1.0/double(5*N_primary);
    double u=unknowns[j];
    double u_combined=u;
    for (unsigned k=0;k<N_primary;k++)
     {    
      if (k!=j)
       {
        u_combined+=unknowns[k]*alpha;
       }
     }
    
    // Contact force
    double fc=unknowns[N_primary+j];
    
    // Eqn numbers
    unsigned displ_eqn=j;
    unsigned contact_eqn=N_primary+j;
    
    // Equilibrium equation
    //=====================
    double N=double(N_primary);
    residuals[displ_eqn]=Stiffness*u_combined-
     Force*(2.0*j+2.0+alpha*N+alpha*N*N-2.0*alpha*j-2.0*alpha)/2.0
     +fc;
    
    // Contact equation
    //=================
 
    // Old version
    //------------
    if (CommandLineArgs::command_line_flag_has_been_set
        ("--old_version"))
     {
      
      // Equation for "melt" rate
      if (u<U_max)
       {
        // "Lower left" quadrant of (u-U_max,fc) space:
        // Linear variation with m forces things back to m=0
        residuals[contact_eqn]=fc;
       }
      else
       {
        // Lower right quadrant of (u-U_max,fc) space
        if (fc<0.0)
         {
          residuals[contact_eqn]=fc;
         }
        // Upper right quadrant of (u-U_max,fc) space
        else
         {
          // Linear variation forces things back to U_max and maintains
          // continuity with upper left and lower right quadrants
          residuals[contact_eqn]=u-U_max; 
         }
       }
     }
    // Completely smooth surface for melt rate residual
    //-------------------------------------------------
    else if (CommandLineArgs::command_line_flag_has_been_set
             ("--completely_smooth"))
     {
      double phi=atan2(fc,(U_max-u));
      if (phi<0) phi+=2.0*MathematicalConstants::Pi;
      
      // Quadratic fit to periodic polynomial that's only
      // above zero between phi=90 and 180 degrees.
      double t0 = 768.0/77.0/0.3141592653589793E1*phi-2176.0/77.0/(
       0.3141592653589793E1*0.3141592653589793E1)*phi*phi+128.0/7.0/(
        0.3141592653589793E1*0.3141592653589793E1*0.3141592653589793E1)*
       phi*phi*phi
       -256.0/77.0/(0.3141592653589793E1*0.3141592653589793E1*
                    0.3141592653589793E1*
                    0.3141592653589793E1)*phi*phi*phi*phi;
      
      residuals[contact_eqn]=(fc*fc+(U_max-u)*(U_max-u))*t0;
     }
    // Single kink
    //------------
    else if (CommandLineArgs::command_line_flag_has_been_set
             ("--single_kink"))
     {
      // Piecewise linear variation with a single kink
      if ((U_max-u)>fc)
       {
        residuals[contact_eqn]=fc;
       }
      else
       {
        residuals[contact_eqn]=(U_max-u);
       }
     }
 
    // Kuhn Tucker
    //------------
    else if (CommandLineArgs::command_line_flag_has_been_set
             ("--kuhn_tucker"))
     {
      residuals[contact_eqn]=fc*(U_max-u);
     }
    else
     {
      oomph_info << "Please specify which method you want to use for\n"
                 << "the enforcement of the contact condition: \n";
      CommandLineArgs::doc_available_flags();
      abort();
     }
   }
 }

References alpha, atan2(), oomph::CommandLineArgs::command_line_flag_has_been_set(), oomph::CommandLineArgs::doc_available_flags(), Force, j, k, N, N_primary, oomph::oomph_info, BiharmonicTestFunctions2::Pi, Stiffness, and U_max.

Referenced by SegregatedFSICollapsibleChannelProblem< ELEMENT >::actions_before_newton_convergence_check(), SegregatedFSICollapsibleChannelProblem< ELEMENT >::actions_before_segregated_convergence_check(), BaseCoupling< M, O >::checkResidual(), oomph::SpineLineMarangoniFluidInterfaceElement< ELEMENT >::fill_in_contribution_to_jacobian(), oomph::FixedVolumeSpineLineMarangoniFluidInterfaceElement< ELEMENT >::fill_in_contribution_to_jacobian(), oomph::SpineAxisymmetricMarangoniSurfactantFluidInterfaceElement< ELEMENT >::fill_in_contribution_to_jacobian(), oomph::SpineLineMarangoniSurfactantFluidInterfaceElement< ELEMENT >::fill_in_contribution_to_jacobian(), BrethertonElement< ELEMENT >::get_jacobian(), TwoNodePoissonElement::get_jacobian(), oomph::Hijacked< ELEMENT >::get_residuals(), main(), BrethertonProblem< ELEMENT >::parameter_study(), plot_it(), set_spurious_solution(), and HomogenisationProblem< ELEMENT >::sub_solve().

plot_it()

void GlobalFct::plot_it

(

const std::string

filename

)

Plot "landscape" of residuals (only for 2D problems!)

 {
  if (N_primary!=1)
   {
    oomph_info << "Skipping plot_it()\n";
    return;
   }
 
  ofstream outfile;
  outfile.open(filename.c_str());
  Vector<double> x(2);
  Vector<double> f(2);
  Vector<double> params(0);
  double x_min=0.0;
  double x_max=10.0;
  double y_min=-5.0;
  double y_max=5.0;
 
  unsigned nplot=100;
  outfile << "ZONE I=" << nplot << ", J=" << nplot << std::endl;
  for (unsigned i=0;i<nplot;i++)
   {
    x[0]=x_min+(x_max-x_min)*double(i)/double(nplot-1);
    for (unsigned j=0;j<nplot;j++)
     {
      x[1]=y_min+(y_max-y_min)*double(j)/double(nplot-1);
      get_residuals(params,x,f);
      outfile << x[0] << " " 
              << x[1] << " " 
              << f[0] << " " 
              << f[1] << "\n" ;
     }
   }
  outfile.close();
 
 }

References f(), MergeRestartFiles::filename, get_residuals(), i, j, N_primary, oomph::oomph_info, Particles2023AnalysisHung::params, and plotDoE::x.

Referenced by oomph::StefanBoltzmannHelper::bin_helper(), oomph::FishMesh< ELEMENT >::build_mesh(), main(), and oomph::StefanBoltzmannHelper::setup_stefan_boltzmann_visibility().

reset()

void GlobalFct::reset

(

Vector< double > &

unknowns

)

Reset unknowns to some nontrivial garbage.

 {
  unsigned n=N_primary*2;
  unknowns.resize(n);
  for (unsigned j=0;j<n;j++)
   {
    unknowns[j]=double(j)+20.5;
   }
 }

References j, n, and N_primary.

Referenced by main().

set_spurious_solution()

void GlobalFct::set_spurious_solution

(

Vector< double > &

unknowns

)

Set unknowns to spurious solution.

 {
  Force=100.0;
  N_primary=1;
  
  Vector<double> param; 
  unknowns.resize(2); 
  Vector<double> residuals(2);
  
  // Displacement
  unknowns[0]=Force/Stiffness+0.5*BlackBoxFDNewtonSolver::Tol;
  
  // Contact force 
  unknowns[1]=-0.5*BlackBoxFDNewtonSolver::Tol;
  
  // Get residuals
  get_residuals(param,unknowns, residuals);
  oomph_info << "Residual for spurious solution: " 
             << residuals[0] << " " 
             << residuals[1] << "\n"; 
 }

References Force, get_residuals(), N_primary, oomph::oomph_info, Stiffness, and oomph::ToleranceForFourierDecomposedHelmholtzOuterBoundary::Tol.

Referenced by main().

Variable Documentation

Force

double GlobalFct::Force =0.0

Force.

Referenced by get_residuals(), main(), and set_spurious_solution().

N_primary

unsigned GlobalFct::N_primary =1

Number of primary variables. Equal to number of raw "displacement" variables for contact problem; total number of degrees of freedom in that case is twice as big because of the unknown contact forces.

Referenced by get_residuals(), main(), plot_it(), reset(), and set_spurious_solution().

Stiffness

double GlobalFct::Stiffness =1.0

Spring stiffness.

Referenced by get_residuals(), and set_spurious_solution().

U_max

double GlobalFct::U_max =5.5

Max. displacement.

Referenced by get_residuals(), and main().