two_d_biharmonic.cc File Reference

#include "math.h"
#include "generic.h"
#include "biharmonic.h"

Classes
class	BiharmonicTestProblem1
class	BiharmonicTestProblem2

Namespaces
	RayParam
	BiharmonicTestFunctions1
	BiharmonicTestFunctions2

Functions
void	BiharmonicTestFunctions1::u_NE (const double &s, double &u)
void	BiharmonicTestFunctions1::dudn_NE (const double &s, double &dudn)
void	BiharmonicTestFunctions1::u_SW (const double &s, double &u)
void	BiharmonicTestFunctions1::dudn_SW (const double &s, double &dudn)
void	BiharmonicTestFunctions1::surface_load (const Vector< double > &x, double &f)
void	BiharmonicTestFunctions1::flux1_NE (const double &s, double &flux1)
void	BiharmonicTestFunctions1::flux0_NE (const double &s, double &flux0)
void	BiharmonicTestFunctions1::solution (const Vector< double > &x, Vector< double > &u)
void	BiharmonicTestFunctions2::boundary_N (const double &s, Vector< double > &r)
void	BiharmonicTestFunctions2::boundary_E (const double &s, Vector< double > &r)
void	BiharmonicTestFunctions2::boundary_S (const double &s, Vector< double > &r)
void	BiharmonicTestFunctions2::boundary_W (const double &s, Vector< double > &r)
void	BiharmonicTestFunctions2::normal_N (const double &s, Vector< double > &n)
void	BiharmonicTestFunctions2::normal_E (const double &s, Vector< double > &n)
void	BiharmonicTestFunctions2::normal_S (const double &s, Vector< double > &n)
void	BiharmonicTestFunctions2::normal_W (const double &s, Vector< double > &n)
void	BiharmonicTestFunctions2::u_N (const double &s, double &u)
void	BiharmonicTestFunctions2::u_E (const double &s, double &u)
void	BiharmonicTestFunctions2::u_S (const double &s, double &u)
void	BiharmonicTestFunctions2::u_W (const double &s, double &u)
double	BiharmonicTestFunctions2::dudx_0 (const Vector< double > x)
double	BiharmonicTestFunctions2::dudx_1 (const Vector< double > x)
void	BiharmonicTestFunctions2::dudn_N (const double &s, double &dudn)
void	BiharmonicTestFunctions2::dudn_E (const double &s, double &dudn)
void	BiharmonicTestFunctions2::dudn_S (const double &s, double &dudn)
void	BiharmonicTestFunctions2::dudn_W (const double &s, double &dudn)
void	BiharmonicTestFunctions2::surface_load (const Vector< double > &x, double &f)
void	BiharmonicTestFunctions2::solution (const Vector< double > &x, Vector< double > &u)
void	print_elemental_jacobian (const unsigned &element_number, const Problem *const problem_pt)
void	print_connectivity_matrix (const Problem *const problem_pt)
int	main (int argc, char *argv[])
	main More...

Variables
bool	RayParam::Unclamp_right_bound = false
unsigned	RayParam::Noel = 0
int	RayParam::X_min = 0
int	RayParam::X_max = 0
int	RayParam::Y_min = 0
int	RayParam::Y_max = 0
bool	RayParam::Print_connectivity_mat = false
bool	RayParam::Print_natural_jacobian = false
bool	RayParam::Print_subblocks = false
bool	RayParam::Print_elemental_jacobian = false
std::string	RayParam::Matbase_str = ""
std::string	RayParam::Connectivity_dir = ""
std::string	RayParam::Natural_jac_dir = ""
std::string	RayParam::Blocked_dir = ""
std::string	RayParam::Elemental_dir = ""
double	BiharmonicTestFunctions2::Pi = MathematicalConstants::Pi
double	BiharmonicTestFunctions2::theta = Pi/4
double	BiharmonicTestFunctions2::r_min = sqrt(1*Pi)
double	BiharmonicTestFunctions2::r_max = sqrt(3*Pi)

Function Documentation

main()

int main	(	int argc	,
		char *	argv[]
	)

main

{
  // number of element
  unsigned n_element = 20;
  
  // Set up doc info
  DocInfo doc_info;
  doc_info.set_directory("RESLT");
  doc_info.number()=0;
 
 
  // Store commandline arguments
  CommandLineArgs::setup(argc,argv);
 
  CommandLineArgs::specify_command_line_flag("--unclamp_right_bound");
 
  // Need to make a string to reflect the type of problem
  CommandLineArgs::specify_command_line_flag("--noel", 
      &RayParam::Noel);
  CommandLineArgs::specify_command_line_flag("--xmin", 
      &RayParam::X_min);
  CommandLineArgs::specify_command_line_flag("--xmax", 
      &RayParam::X_max);
  CommandLineArgs::specify_command_line_flag("--ymin", 
      &RayParam::Y_min);
  CommandLineArgs::specify_command_line_flag("--ymax", 
      &RayParam::Y_max);
 
  CommandLineArgs::specify_command_line_flag("--connectivity_mat",
      &RayParam::Connectivity_dir);
  CommandLineArgs::specify_command_line_flag("--natural_jacobian",
      &RayParam::Natural_jac_dir);
  CommandLineArgs::specify_command_line_flag("--sub_blocks",
      &RayParam::Blocked_dir);
  CommandLineArgs::specify_command_line_flag("--elemental_jacobian",
      &RayParam::Elemental_dir);
 
  // Parse the above flags.
  CommandLineArgs::parse_and_assign();
  CommandLineArgs::doc_specified_flags();
 
 
 
  // If there is more than one argument, we Milan's tests.
  if(argc > 1)
  {
 
  if(CommandLineArgs::command_line_flag_has_been_set(
        "--unclamp_right_bound"))
  {
    RayParam::Unclamp_right_bound = true;
    RayParam::Matbase_str = "twod_bihar_right_bound_neumann";
  }
  else
  {
    RayParam::Unclamp_right_bound = false;
    RayParam::Matbase_str = "twod_bihar_clamp_all_bound";
  }
 
  if(!CommandLineArgs::command_line_flag_has_been_set("--noel"))
  {
    std::ostringstream err_stream;
    err_stream << "Please set --noel." << std::endl;
    throw OomphLibError(err_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
  else
  {
    std::ostringstream str_ss;
    str_ss << RayParam::Matbase_str << "_Noel" << RayParam::Noel;
    RayParam::Matbase_str = str_ss.str();
  }
 
  if(!CommandLineArgs::command_line_flag_has_been_set("--xmin"))
  {
    std::ostringstream err_stream;
    err_stream << "Please set --xmin." << std::endl;
    throw OomphLibError(err_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
  if(!CommandLineArgs::command_line_flag_has_been_set("--xmax"))
  {
    std::ostringstream err_stream;
    err_stream << "Please set --xmax." << std::endl;
    throw OomphLibError(err_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
  if(!CommandLineArgs::command_line_flag_has_been_set("--ymin"))
  {
    std::ostringstream err_stream;
    err_stream << "Please set --ymin." << std::endl;
    throw OomphLibError(err_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
  if(!CommandLineArgs::command_line_flag_has_been_set("--ymax"))
  {
    std::ostringstream err_stream;
    err_stream << "Please set --ymax." << std::endl;
    throw OomphLibError(err_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }
 
  if(CommandLineArgs::command_line_flag_has_been_set("--connectivity_mat"))
  {
    RayParam::Print_connectivity_mat = true;
  }
  else
  {
    RayParam::Print_connectivity_mat = false;
  }
 
  if(CommandLineArgs::command_line_flag_has_been_set("--natural_jacobian"))
  {
    RayParam::Print_natural_jacobian = true;
  }
  else
  {
    RayParam::Print_natural_jacobian = false;
  }
 
  if(CommandLineArgs::command_line_flag_has_been_set("--sub_blocks"))
  {
    RayParam::Print_subblocks = true;
  }
  else
  {
    RayParam::Print_subblocks = false;
  }
 
  if(CommandLineArgs::command_line_flag_has_been_set("--elemental_jacobian"))
  {
    RayParam::Print_elemental_jacobian = true;
  }
  else
  {
    RayParam::Print_elemental_jacobian = false;
  }
 
 
 
    oomph_info << "n_element: " << RayParam::Noel << std::endl;
    oomph_info << "x_min: " << RayParam::X_min << std::endl;
    oomph_info << "x_max: " << RayParam::X_max << std::endl;
    oomph_info << "y_min: " << RayParam::Y_min << std::endl;
    oomph_info << "y_max: " << RayParam::Y_max << std::endl;
 
    BiharmonicTestProblem1 problem(RayParam::X_min,
                                   RayParam::X_max,
                                   RayParam::Y_min,
                                   RayParam::Y_max,
                                   RayParam::Noel);
 
    BiharmonicPreconditioner my_prec;
    my_prec.bulk_element_mesh_pt() = problem.bulk_element_mesh_pt();
    my_prec.preconditioner_type() = 0;
 
    IterativeLinearSolver* solver_pt = new  CG<CRDoubleMatrix>;
    solver_pt->preconditioner_pt() = &my_prec;
 
    // Apply the solver
    problem.linear_solver_pt() = solver_pt;
 
    if(RayParam::Print_natural_jacobian)
    {
      std::ostringstream dump_stream;
      dump_stream << RayParam::Natural_jac_dir 
                  << "/"
                  << RayParam::Matbase_str;
      std::string dump_str = dump_stream.str();
      problem.dump_jacobian(dump_str);
    }
 
    if(RayParam::Print_connectivity_mat)
    {
      print_connectivity_matrix(&problem);
    }
 
//    if(RayParam::Print_elemental_jacobian)
//    {
//      //print_elemental_jacobian(&problem);
//    }
 
    if(RayParam::Print_subblocks)
    {
      std::ostringstream dump_stream;
      dump_stream << RayParam::Blocked_dir 
                  << "/"
                  << RayParam::Matbase_str;
      std::string dump_str = dump_stream.str();
 
//      my_prec.print_subblocks(dump_str);
 
      std::ostringstream dump_ss;
      dump_ss << RayParam::Natural_jac_dir
              << "/"
              << RayParam::Matbase_str;
      std::string another_dump_str = dump_ss.str();
 
//      my_prec.set_fullblock_dir(another_dump_str);
 
      problem.newton_solve();
 
    }
  }
  else
  {
    // Biharmonic Problem 1 (square)
    // Exact Biharmonic Preconditioner
    {
      oomph_info 
        << "/////////////////////////////////////////////////////////////////////"
        << std::endl;
      oomph_info << "TESTING: Square 2D Biharmonic Problem w/ "
        << "Exact Preconditioning"
        << std::endl;
      oomph_info 
        << "/////////////////////////////////////////////////////////////////////"
        << std::endl;
 
      // create the problem
      BiharmonicTestProblem1 problem(n_element);
 
      // setup the preconditioner
      BiharmonicPreconditioner* prec_pt = new BiharmonicPreconditioner;
      prec_pt->bulk_element_mesh_pt() = problem.bulk_element_mesh_pt();
      prec_pt->preconditioner_type() = 0;
 
      // setup the solver
      IterativeLinearSolver* solver_pt = new CG<CRDoubleMatrix>;  
      solver_pt->preconditioner_pt() = prec_pt;
 
      // apply the solver
      problem.linear_solver_pt() = solver_pt;
      problem.newton_solve();
 
      // ouput the solution
      problem.doc_solution(doc_info,BiharmonicTestFunctions1::solution);
      doc_info.number()++;
 
      // clean up
      delete solver_pt;
      delete prec_pt;
    }
 
    // Biharmonic Problem 2 (section of annulus)
    // Inexact Biharmonic Preconditioner w/ SuperLU
    {
      oomph_info 
        << "/////////////////////////////////////////////////////////////////////"
        << std::endl;
      oomph_info << "TESTING: Curved 2D Biharmonic Problem w/ "
        << "Inexact Preconditioning"
        << std::endl;
      oomph_info 
        << "/////////////////////////////////////////////////////////////////////"
        << std::endl;
 
      // create the problem
      BiharmonicTestProblem2 problem(n_element);
 
      // setup the preconditioner
      BiharmonicPreconditioner* prec_pt = new BiharmonicPreconditioner;
      prec_pt->bulk_element_mesh_pt() = problem.bulk_element_mesh_pt();
      prec_pt->preconditioner_type() = 1;
 
      // setup the solver
      IterativeLinearSolver* solver_pt = new CG<CRDoubleMatrix>;  
      solver_pt->preconditioner_pt() = prec_pt;
 
      // apply the solver
      problem.linear_solver_pt() = solver_pt;
      problem.newton_solve();
 
      // ouput the solution
      problem.doc_solution(doc_info,BiharmonicTestFunctions2::solution);
      doc_info.number()++;
 
      // clean up
      delete solver_pt;
      delete prec_pt;
    }
  }
}

References RayParam::Blocked_dir, oomph::BiharmonicPreconditioner::bulk_element_mesh_pt(), oomph::CommandLineArgs::command_line_flag_has_been_set(), RayParam::Connectivity_dir, oomph::CommandLineArgs::doc_specified_flags(), RayParam::Elemental_dir, RayParam::Matbase_str, RayParam::Natural_jac_dir, RayParam::Noel, oomph::DocInfo::number(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::CommandLineArgs::parse_and_assign(), oomph::IterativeLinearSolver::preconditioner_pt(), oomph::BiharmonicPreconditioner::preconditioner_type(), RayParam::Print_connectivity_mat, print_connectivity_matrix(), RayParam::Print_elemental_jacobian, RayParam::Print_natural_jacobian, RayParam::Print_subblocks, problem, oomph::DocInfo::set_directory(), oomph::CommandLineArgs::setup(), BiharmonicTestFunctions1::solution(), BiharmonicTestFunctions2::solution(), oomph::CommandLineArgs::specify_command_line_flag(), oomph::Global_string_for_annotation::string(), RayParam::Unclamp_right_bound, RayParam::X_max, RayParam::X_min, RayParam::Y_max, and RayParam::Y_min.

print_connectivity_matrix()

void print_connectivity_matrix

(

const Problem *const

problem_pt

)

{
  const Mesh* const mesh_pt = problem_pt->mesh_pt();
 
  const unsigned n_element = mesh_pt->nelement();
//  const unsigned n_ele_1d = sqrt(n_element);
 
  // Set up out file
  std::ostringstream filenamestream;
  filenamestream << RayParam::Connectivity_dir
    << "/"
    << RayParam::Matbase_str;
 
 
  std::ofstream outfile;
  outfile.open(filenamestream.str().c_str());
 
  // Loop through the number of elements
  for (unsigned ele_i = 0; ele_i < n_element; ele_i++) 
  {
    outfile << "Element number: " << ele_i << std::endl;
 
    // Get pointer to the element
    FiniteElement* elem_pt 
      = mesh_pt->finite_element_pt(ele_i);
 
    unsigned nnod=elem_pt->nnode();
 
    for (unsigned nod_i = 0; nod_i < nnod; nod_i++) 
    {
      outfile << "Node number: " << nod_i << std::endl;
 
      // Get the node
      Node* nod_pt = elem_pt->node_pt(nod_i);
      const unsigned nval = nod_pt->nvalue();
 
      outfile << "x: " << nod_pt->x(0) 
              << ", y: " << nod_pt->x(1) 
              << std::endl;
 
      for (unsigned val_i = 0; val_i < nval; val_i++) 
      {
        outfile << val_i << " ";
        long eqn_num = nod_pt->eqn_number(val_i);
        outfile << eqn_num << " ";
        if(!nod_pt->is_pinned(val_i))
        {
          outfile << elem_pt->local_eqn_number(eqn_num) << " ";
        }
 
        outfile << "\n";
      } // for values in node
    } // for node
  } // for elements
  outfile.close();
}

References RayParam::Connectivity_dir, oomph::Data::eqn_number(), oomph::Mesh::finite_element_pt(), oomph::Data::is_pinned(), oomph::GeneralisedElement::local_eqn_number(), RayParam::Matbase_str, oomph::Problem::mesh_pt(), oomph::Mesh::nelement(), oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), oomph::Data::nvalue(), and oomph::Node::x().

Referenced by main().

print_elemental_jacobian()

void print_elemental_jacobian	(	const unsigned &	element_number,
		const Problem *const	problem_pt
	)

{
  AssemblyHandler* const assembly_handler_pt 
    = problem_pt->assembly_handler_pt();
 
  const Mesh* const mesh_pt = problem_pt->mesh_pt();
 
  const unsigned n_element = mesh_pt->nelement();
  const unsigned n_ele_1d = sqrt(n_element);
 
  oomph_info << "Elements: 1D: " << n_ele_1d 
             << ", total: " << n_element << std::endl;
 
#ifdef PARANOID
  if(element_number >= n_element)
  {
    std::ostringstream err_stream;
    err_stream << "Supplied element number: " << element_number << ",\n"
      << "But number of elements is: " << n_element << std::endl;
    throw OomphLibError(err_stream.str(),
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
  }
#endif
 
  // Get pointer to the element
  GeneralisedElement* elem_pt 
    = mesh_pt->element_pt(element_number);
 
  // Find number of dofs in the element
  const unsigned n_element_dofs = assembly_handler_pt->ndof(elem_pt);
 
  // Set up an array
  Vector<double> element_residuals(n_element_dofs);
 
  // Set up a matrix
  DenseMatrix<double> element_jacobian(n_element_dofs);
 
  // Fill the array
  assembly_handler_pt->get_jacobian(elem_pt,
      element_residuals,
      element_jacobian);
 
  std::ostringstream ele_stream;
  ele_stream << "N"<< n_ele_1d << "_ele_" << element_number;
  element_jacobian.sparse_indexed_output(ele_stream.str(),15,true);
 
  //  // Output the coordinates just to make sure...
  //  FiniteElement* finite_element_pt 
  //    = mesh_pt->finite_element_pt(element_number);
  //
  //  const unsigned n_node = finite_element_pt->nnode();
  //
  //  // Loop over nodes and output coordinates
  //  for (unsigned n = 0; n < n_node; n++) 
  //  {
  //    Node* nod_pt = finite_element_pt->node_pt(n);
  //    const double x = nod_pt->x(0);
  //    const double y = nod_pt->x(1);
  //    oomph_info << "Node: " << n << ", " << x << ", " << y << std::endl;
  //  }
}

References oomph::Problem::assembly_handler_pt(), oomph::Mesh::element_pt(), oomph::AssemblyHandler::get_jacobian(), oomph::Problem::mesh_pt(), oomph::AssemblyHandler::ndof(), oomph::Mesh::nelement(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::Matrix< T, MATRIX_TYPE >::sparse_indexed_output(), and sqrt().