snap_mesh.cc File Reference

#include "generic.h"
#include "solid.h"
#include "constitutive.h"
#include "navier_stokes.h"
#include "poisson.h"
#include "linear_elasticity.h"
#include "mesh_smoothing.h"
#include "meshes/tetgen_mesh.h"

Functions
void	demo_smoothing_with_linear_elasticity ()
void	demo_smoothing_with_poisson ()
void	demo_smoothing_with_nonlinear_elasticity ()
int	main (int argc, char **argv)

Function Documentation

demo_smoothing_with_linear_elasticity()

void demo_smoothing_with_linear_elasticity

(

)

=============================================================

Demonstrate mesh smoothing with linear elasticity

Storage for mesh

 {
 
 // Label for output
 DocInfo doc_info;
 
 
 // Read boundary enumeration, generated by xda_to_poly_fsi
 //--------------------------------------------------------
 string input_string;
 
 // Open input file
 ifstream* input_file_pt=new 
  ifstream("short_coarse_mesh_files/boundary_enumeration.dat"); 
 
 // Check if it's been opened succesfully
 if (input_file_pt==0)
  {
   oomph_info << "ERROR while trying to open boundary enumeration file " 
              << std::endl;
   exit(1); 
  }
   
 getline(*input_file_pt,input_string,' ');
 unsigned fluid_fsi_lo=atoi(input_string.c_str());
 getline(*input_file_pt,input_string,'#');
 unsigned fluid_fsi_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 //unsigned fluid_in_lo=atoi(input_string.c_str());
 getline(*input_file_pt,input_string,'#');
 //unsigned fluid_in_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 //unsigned fluid_out_lo=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 getline(*input_file_pt,input_string,' ');
 getline(*input_file_pt,input_string,'#');
 //unsigned fluid_out_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 unsigned solid_fsi_lo=atoi(input_string.c_str());
 getline(*input_file_pt,input_string,'#');
 unsigned solid_fsi_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 //unsigned solid_pin_lo=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 input_file_pt->ignore(80,'\n');
 getline(*input_file_pt,input_string,' ');
 getline(*input_file_pt,input_string,'#');
 //unsigned solid_pin_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 unsigned solid_outer_lo=atoi(input_string.c_str());
 getline(*input_file_pt,input_string,'#');
 unsigned solid_outer_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 delete input_file_pt;
 
 // Identify boundaries that make up the FSI interfaces in the
 // fluid and solid meshes
 Vector<unsigned> fluid_fsi_boundary_id; 
 Vector<unsigned> solid_fsi_boundary_id;
 Vector<unsigned> solid_outer_boundary_id;
 
 // The FSI boundaries of the fluid mesh
 unsigned nfluid_fsi=fluid_fsi_hi-fluid_fsi_lo+1;
 fluid_fsi_boundary_id.resize(nfluid_fsi);
 for(unsigned i=0; i<nfluid_fsi; i++)
  {
   fluid_fsi_boundary_id[i]=fluid_fsi_lo+i;
  }
 
  // The solid and fluid fsi boundaries are numbered in the same way.
 unsigned nsolid_fsi=solid_fsi_hi-solid_fsi_lo+1;
 solid_fsi_boundary_id.resize(nsolid_fsi);
 for(unsigned i=0; i<nsolid_fsi; i++)
  {
   solid_fsi_boundary_id[i]=solid_fsi_lo+i;
  }
 
 // Outer solid boundary
 unsigned nsolid_outer=solid_outer_hi-solid_outer_lo+1;
 solid_outer_boundary_id.resize(nsolid_outer);
 for(unsigned i=0; i<nsolid_outer; i++)
  {
   solid_outer_boundary_id[i]=solid_outer_lo+i;
  }
 
 
 // Variables to create meshes
 string node_file_name;
 string element_file_name;
 string face_file_name;
 bool split_corner_elements;
 bool switch_normal;
 
 SolidTetgenMesh<TPVDElement<3,3> >* Orig_mesh_pt;
 
 // Loop over fluid and solid meshes
 for (unsigned do_solid=0;do_solid<2;do_solid++)
  {
 
   // Build original meshes
   if (do_solid==1)
    {
 
     oomph_info << "Doing solid\n";
 
     // Output directory
     doc_info.set_directory("RESLT_solid_with_linear_elasticity");
 
     //Create bulk mesh, sub-dividing "corner" elements
     node_file_name="short_coarse_mesh_files/solid.1.node";
     element_file_name="short_coarse_mesh_files/solid.1.ele";
     face_file_name="short_coarse_mesh_files/solid.1.face";
     split_corner_elements=true;
     Orig_mesh_pt = new SolidTetgenMesh<TPVDElement<3,3> >(node_file_name,
                                                        element_file_name,
                                                        face_file_name,
                                                        split_corner_elements);
    }
   else
    {
 
     oomph_info << "Doing fluid\n";
 
     // Output directory
     doc_info.set_directory("RESLT_fluid_with_linear_elasticity");
     
     //Create bulk mesh, sub-dividing "corner" elements
     node_file_name="short_coarse_mesh_files/fluid.1.node";
     element_file_name="short_coarse_mesh_files/fluid.1.ele";
     face_file_name="short_coarse_mesh_files/fluid.1.face";
     split_corner_elements=true;
     switch_normal=true;
     Orig_mesh_pt =  new SolidTetgenMesh<TPVDElement<3,3> >(node_file_name,
                                                         element_file_name,
                                                         face_file_name,
                                                         split_corner_elements,
                                                         switch_normal);
    }
   
   
   // Check for inverted elements before quadratic trafo
   std::string filename;
   filename=doc_info.directory();
   filename+="/mesh_before_snap.dat";
   Orig_mesh_pt->output(filename.c_str());
   bool mesh_has_inverted_elements;
   Orig_mesh_pt->check_inverted_elements(mesh_has_inverted_elements);  
   cout << "Before quadratic snapping mesh does ";
   if (!mesh_has_inverted_elements) cout << "not ";
   cout << "have inverted elements. \n";
   
 
   if (do_solid==1) 
    {
   // Snap FSI interface to quadratic surface
     switch_normal=false; 
     Orig_mesh_pt->template snap_to_quadratic_surface<TPVDElement<3,3> >(
      solid_fsi_boundary_id,
      "short_coarse_mesh_files/quadratic_fsi_boundary.dat",
      switch_normal);
 
     // Snap outer solid interface to quadratic surface
     switch_normal=true;
     Orig_mesh_pt->template snap_to_quadratic_surface<TPVDElement<3,3> >(
      solid_outer_boundary_id,
      "short_coarse_mesh_files/quadratic_outer_solid_boundary.dat",
      switch_normal);
    }
   else
    {
     // Snap FSI interface to quadratic surface
     switch_normal=false; 
     Orig_mesh_pt->template snap_to_quadratic_surface<TPVDElement<3,3> >(
      fluid_fsi_boundary_id,
      "short_coarse_mesh_files/quadratic_fsi_boundary.dat",
      switch_normal);
    }
   
   
   // Check
   filename=doc_info.directory();
   filename+="/mesh_after_snap.dat";
   Orig_mesh_pt->output(filename.c_str());
   std::ofstream inverted_mesh_elements;
   filename=doc_info.directory();
   filename+="/inverted_elements_after_snap.dat";
   inverted_mesh_elements.open(filename.c_str());
   Orig_mesh_pt->check_inverted_elements(mesh_has_inverted_elements,
                                         inverted_mesh_elements);  
   inverted_mesh_elements.close(); 
   cout << "After quadratic snapping mesh does ";
   if (!mesh_has_inverted_elements) cout << "not ";
   cout << "have inverted elements. \n";
   
 
   
   // Smooth it: 
   //-----------
   
   
   // Create set containing the nodes whose displacements are contrained
   // (nodes on quadratic boundaries)
   std::set<Node*> pinned_nodes;
   
   if (do_solid==1) 
    {
     
     // Loop over nodes on the FSI boundary in the solid mesh
     unsigned nbound_pinned=solid_fsi_boundary_id.size();
     for(unsigned i=0;i<nbound_pinned;i++)
      {
       //Get the mesh boundary
       unsigned b = solid_fsi_boundary_id[i];
       unsigned num_nod=Orig_mesh_pt->nboundary_node(b);
       for (unsigned inod=0;inod<num_nod;inod++)
        {
         // Get node
         pinned_nodes.insert(Orig_mesh_pt->boundary_node_pt(b,inod));
        }
      }
     
     // Loop over nodes on outer boundary of solid mesh
     nbound_pinned=solid_outer_boundary_id.size();
     for(unsigned i=0;i<nbound_pinned;i++)
      {
       //Get the mesh boundary
       unsigned b = solid_outer_boundary_id[i];
       unsigned num_nod=Orig_mesh_pt->nboundary_node(b);
       for (unsigned inod=0;inod<num_nod;inod++)
        {
         // Get node     
         pinned_nodes.insert(Orig_mesh_pt->boundary_node_pt(b,inod));
        }
      }
     
    }
   else
    {
     // Loop over nodes on the FSI boundary in the solid mesh
     unsigned nbound_pinned=fluid_fsi_boundary_id.size();
     for(unsigned i=0;i<nbound_pinned;i++)
      {
       //Get the mesh boundary
       unsigned b = fluid_fsi_boundary_id[i];
       unsigned num_nod=Orig_mesh_pt->nboundary_node(b);
       for (unsigned inod=0;inod<num_nod;inod++)
        {
         // Get node
         pinned_nodes.insert(Orig_mesh_pt->boundary_node_pt(b,inod));
        }
      }
    }
     
   // Smooth
   LinearElasticitySmoothMesh<TLinearElasticityElement<3,3> >()(Orig_mesh_pt,
                                                                pinned_nodes);
 
   // Check
   filename=doc_info.directory();
   filename+="/mesh_after_smooth.dat";
   Orig_mesh_pt->output(filename.c_str());
   filename=doc_info.directory();
   filename+="/inverted_elements_after_smooth.dat";
   inverted_mesh_elements.open(filename.c_str());
   Orig_mesh_pt->check_inverted_elements(mesh_has_inverted_elements,
                                         inverted_mesh_elements);  
   inverted_mesh_elements.close(); 
   cout << "After smoothing mesh does ";
   if (!mesh_has_inverted_elements) cout << "not ";
   cout << "have inverted elements. \n";
 
   // Cleanup
   delete Orig_mesh_pt;
   Orig_mesh_pt=0;
 
  }
 
 }

References b, oomph::SolidMesh::boundary_node_pt(), oomph::Mesh::check_inverted_elements(), oomph::DocInfo::directory(), MergeRestartFiles::filename, i, oomph::Mesh::nboundary_node(), oomph::oomph_info, oomph::Mesh::output(), oomph::DocInfo::set_directory(), and oomph::Global_string_for_annotation::string().

Referenced by main().

demo_smoothing_with_nonlinear_elasticity()

void demo_smoothing_with_nonlinear_elasticity

(

)

/////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////// =============================================================

Demonstrate smoothing with nonlinear elasticity

Storage for mesh and copy

...and another copy

...and another copy

 {
 
 // Label for output
 DocInfo doc_info;
 
 
 // Read boundary enumeration, generated by xda_to_poly_fsi
 //--------------------------------------------------------
 string input_string;
 
 // Open input file
 ifstream* input_file_pt=new 
  ifstream("short_coarse_mesh_files/boundary_enumeration.dat"); 
 
 // Check if it's been opened succesfully
 if (input_file_pt==0)
  {
   oomph_info << "ERROR while trying to open boundary enumeration file " 
              << std::endl;
   exit(1); 
  }
   
 
 getline(*input_file_pt,input_string,' ');
 unsigned fluid_fsi_lo=atoi(input_string.c_str());
 getline(*input_file_pt,input_string,'#');
 unsigned fluid_fsi_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 //unsigned fluid_in_lo=atoi(input_string.c_str());
 getline(*input_file_pt,input_string,'#');
 //unsigned fluid_in_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 //unsigned fluid_out_lo=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 getline(*input_file_pt,input_string,' ');
 getline(*input_file_pt,input_string,'#');
 //unsigned fluid_out_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 unsigned solid_fsi_lo=atoi(input_string.c_str());
 getline(*input_file_pt,input_string,'#');
 unsigned solid_fsi_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 //unsigned solid_pin_lo=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 input_file_pt->ignore(80,'\n');
 getline(*input_file_pt,input_string,' ');
 getline(*input_file_pt,input_string,'#');
 //unsigned solid_pin_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 //unsigned solid_outer_lo=atoi(input_string.c_str());
 getline(*input_file_pt,input_string,'#');
 //unsigned solid_outer_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 delete input_file_pt;
 
 // Identify boundaries that make up the FSI interfaces in the
 // fluid and solid meshes
 Vector<unsigned> fluid_fsi_boundary_id; 
 Vector<unsigned> solid_fsi_boundary_id;
 
 // The FSI boundaries of the fluid mesh
 unsigned nfluid_fsi=fluid_fsi_hi-fluid_fsi_lo+1;
 fluid_fsi_boundary_id.resize(nfluid_fsi);
 for(unsigned i=0; i<nfluid_fsi; i++)
  {
   fluid_fsi_boundary_id[i]=fluid_fsi_lo+i;
  }
 
  // The solid and fluid fsi boundaries are numbered in the same way.
 unsigned nsolid_fsi=solid_fsi_hi-solid_fsi_lo+1;
 solid_fsi_boundary_id.resize(nsolid_fsi);
 for(unsigned i=0; i<nsolid_fsi; i++)
  {
   solid_fsi_boundary_id[i]=solid_fsi_lo+i;
  }
 
 // Variables to create meshes
 string node_file_name;
 string element_file_name;
 string face_file_name;
 bool split_corner_elements;
 bool switch_normal;
 
 SolidTetgenMesh<TPVDElement<3,3> >* Orig_mesh_pt;
 SolidTetgenMesh<TPVDElement<3,3> >* Copy_of_mesh_pt;
 
 // Loop over fluid and solid meshes
 for (unsigned do_solid=0;do_solid<2;do_solid++)
  {
 
   // Build original meshes
   if (do_solid==1)
    {
 
     oomph_info << "Doing solid\n";
 
     // Output directory
     doc_info.set_directory("RESLT_solid");
 
     //Create bulk mesh, sub-dividing "corner" elements
     node_file_name="short_coarse_mesh_files/solid.1.node";
     element_file_name="short_coarse_mesh_files/solid.1.ele";
     face_file_name="short_coarse_mesh_files/solid.1.face";
     split_corner_elements=true;
     Orig_mesh_pt = new SolidTetgenMesh<TPVDElement<3,3> >(node_file_name,
                                                        element_file_name,
                                                        face_file_name,
                                                        split_corner_elements);
     Copy_of_mesh_pt=new SolidTetgenMesh<TPVDElement<3,3> >(node_file_name,
                                                         element_file_name,
                                                         face_file_name,
                                                         split_corner_elements);
    }
   else
    {
 
     oomph_info << "Doing fluid\n";
 
     // Output directory
     doc_info.set_directory("RESLT_fluid");
     
     //Create bulk mesh, sub-dividing "corner" elements
     node_file_name="short_coarse_mesh_files/fluid.1.node";
     element_file_name="short_coarse_mesh_files/fluid.1.ele";
     face_file_name="short_coarse_mesh_files/fluid.1.face";
     split_corner_elements=true;
     switch_normal=true;
     Orig_mesh_pt =  new SolidTetgenMesh<TPVDElement<3,3> >(node_file_name,
                                                         element_file_name,
                                                         face_file_name,
                                                         split_corner_elements,
                                                         switch_normal);
     Copy_of_mesh_pt=new SolidTetgenMesh<TPVDElement<3,3> >(node_file_name,
                                                         element_file_name,
                                                         face_file_name,
                                                         split_corner_elements,
                                                         switch_normal);
    }
   
   
   // Check for inverted elements before quadratic trafo
   std::string filename;
   filename=doc_info.directory();
   filename+="/mesh_before_snap.dat";
   Orig_mesh_pt->output(filename.c_str());
   bool mesh_has_inverted_elements;
   Orig_mesh_pt->check_inverted_elements(mesh_has_inverted_elements);  
   cout << "Before quadratic snapping mesh does ";
   if (!mesh_has_inverted_elements) cout << "not ";
   cout << "have inverted elements. \n";
   
   // Snap FSI interface to quadratic surface
   if (do_solid==1) 
    {
     switch_normal=false; 
     Orig_mesh_pt->template snap_to_quadratic_surface<TPVDElement<3,3> >(
      solid_fsi_boundary_id,
      "short_coarse_mesh_files/quadratic_fsi_boundary.dat",
      switch_normal);
    }
   else
    {
     switch_normal=false; 
     Orig_mesh_pt->template snap_to_quadratic_surface<TPVDElement<3,3> >(
      fluid_fsi_boundary_id,
      "short_coarse_mesh_files/quadratic_fsi_boundary.dat",
      switch_normal);
    }
   
   
   // Check
   filename=doc_info.directory();
   filename+="/mesh_after_snap.dat";
   Orig_mesh_pt->output(filename.c_str());
   std::ofstream inverted_mesh_elements;
   filename=doc_info.directory();
   filename+="/inverted_elements_after_snap.dat";
   inverted_mesh_elements.open(filename.c_str());
   Orig_mesh_pt->check_inverted_elements(mesh_has_inverted_elements,
                                         inverted_mesh_elements);  
   inverted_mesh_elements.close(); 
   cout << "After quadratic snapping mesh does ";
   if (!mesh_has_inverted_elements) cout << "not ";
   cout << "have inverted elements. \n";
   
   // Smooth it: Specify the two meshes (the original and a working
   // copy that can get deleted afterwards and the boundary IDs
   // of the fsi boundary (same in fluid and solid!) 
   unsigned max_steps=100;
   if (CommandLineArgs::Argc!=1)
    {
     cout << "Validation, just doing one step of mesh smoothing\n";
     max_steps=1;
    }
   NonLinearElasticitySmoothMesh<TPVDElement<3,3> >()(Orig_mesh_pt,
                                                      Copy_of_mesh_pt,
                                                      solid_fsi_boundary_id,
                                                      doc_info,
                                                      max_steps);
   
 
   // Check
   filename=doc_info.directory();
   filename+="/mesh_after_smooth.dat";
   Orig_mesh_pt->output(filename.c_str());
 
   filename=doc_info.directory();
   filename+="/inverted_elements_after_smooth.dat";
   inverted_mesh_elements.open(filename.c_str());
   Orig_mesh_pt->check_inverted_elements(mesh_has_inverted_elements,
                                         inverted_mesh_elements);  
   inverted_mesh_elements.close(); 
   cout << "After smoothing mesh does ";
   if (!mesh_has_inverted_elements) cout << "not ";
   cout << "have inverted elements. \n";
 
 
   filename=doc_info.directory();
   filename+="/new_nodal_positions.dat";
   ofstream node_file;
   node_file.open(filename.c_str());
   unsigned nnod=Orig_mesh_pt->nnode();
   for (unsigned j=0;j<nnod;j++)    
    {
     Node* nod_pt=Orig_mesh_pt->node_pt(j);
     node_file << nod_pt->x(0) << " " 
               << nod_pt->x(1) << " " 
               << nod_pt->x(2) << " " 
               << std::endl;
    }
   node_file.close();
    
    
   // Cleanup
   delete Orig_mesh_pt;
   Orig_mesh_pt=0;
   delete Copy_of_mesh_pt;
   Copy_of_mesh_pt=0;
   
  }
 
}

References oomph::CommandLineArgs::Argc, oomph::Mesh::check_inverted_elements(), oomph::DocInfo::directory(), MergeRestartFiles::filename, i, j, oomph::Mesh::nnode(), oomph::SolidMesh::node_pt(), oomph::oomph_info, oomph::Mesh::output(), oomph::DocInfo::set_directory(), oomph::Global_string_for_annotation::string(), and oomph::Node::x().

Referenced by main().

demo_smoothing_with_poisson()

void demo_smoothing_with_poisson

(

)

/////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////// =============================================================

Demonstrate smoothing with Poisson

Storage for mesh

 {
 
 // Label for output
 DocInfo doc_info;
 
 
 // Read boundary enumeration, generated by xda_to_poly_fsi
 //--------------------------------------------------------
 string input_string;
 
 // Open input file
 ifstream* input_file_pt=new 
  ifstream("short_coarse_mesh_files/boundary_enumeration.dat"); 
 
 // Check if it's been opened succesfully
 if (input_file_pt==0)
  {
   oomph_info << "ERROR while trying to open boundary enumeration file " 
              << std::endl;
   exit(1); 
  }
   
 
 getline(*input_file_pt,input_string,' ');
 unsigned fluid_fsi_lo=atoi(input_string.c_str());
 getline(*input_file_pt,input_string,'#');
 unsigned fluid_fsi_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 //unsigned fluid_in_lo=atoi(input_string.c_str());
 getline(*input_file_pt,input_string,'#');
 //unsigned fluid_in_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 //unsigned fluid_out_lo=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 getline(*input_file_pt,input_string,' ');
 getline(*input_file_pt,input_string,'#');
 //unsigned fluid_out_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 unsigned solid_fsi_lo=atoi(input_string.c_str());
 getline(*input_file_pt,input_string,'#');
 unsigned solid_fsi_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 //unsigned solid_pin_lo=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 input_file_pt->ignore(80,'\n');
 getline(*input_file_pt,input_string,' ');
 getline(*input_file_pt,input_string,'#');
 //unsigned solid_pin_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 getline(*input_file_pt,input_string,' ');
 unsigned solid_outer_lo=atoi(input_string.c_str());
 getline(*input_file_pt,input_string,'#');
 unsigned solid_outer_hi=atoi(input_string.c_str());
 input_file_pt->ignore(80,'\n');
 
 delete input_file_pt;
 
 // Identify boundaries that make up the FSI interfaces in the
 // fluid and solid meshes
 Vector<unsigned> fluid_fsi_boundary_id; 
 Vector<unsigned> solid_fsi_boundary_id;
 Vector<unsigned> solid_outer_boundary_id;
 
 // The FSI boundaries of the fluid mesh
 unsigned nfluid_fsi=fluid_fsi_hi-fluid_fsi_lo+1;
 fluid_fsi_boundary_id.resize(nfluid_fsi);
 for(unsigned i=0; i<nfluid_fsi; i++)
  {
   fluid_fsi_boundary_id[i]=fluid_fsi_lo+i;
  }
 
  // The solid and fluid fsi boundaries are numbered in the same way.
 unsigned nsolid_fsi=solid_fsi_hi-solid_fsi_lo+1;
 solid_fsi_boundary_id.resize(nsolid_fsi);
 for(unsigned i=0; i<nsolid_fsi; i++)
  {
   solid_fsi_boundary_id[i]=solid_fsi_lo+i;
  }
 
 // Outer solid boundary
 unsigned nsolid_outer=solid_outer_hi-solid_outer_lo+1;
 solid_outer_boundary_id.resize(nsolid_outer);
 for(unsigned i=0; i<nsolid_outer; i++)
  {
   solid_outer_boundary_id[i]=solid_outer_lo+i;
  }
 
 
 // Variables to create meshes
 string node_file_name;
 string element_file_name;
 string face_file_name;
 bool split_corner_elements;
 bool switch_normal;
 
 SolidTetgenMesh<TPVDElement<3,3> >* Orig_mesh_pt;
 
 // Loop over fluid and solid meshes
 for (unsigned do_solid=0;do_solid<2;do_solid++)
  {
 
   // Build original meshes
   if (do_solid==1)
    {
 
     oomph_info << "Doing solid\n";
 
     // Output directory
     doc_info.set_directory("RESLT_solid_with_poisson");
 
     //Create bulk mesh, sub-dividing "corner" elements
     node_file_name="short_coarse_mesh_files/solid.1.node";
     element_file_name="short_coarse_mesh_files/solid.1.ele";
     face_file_name="short_coarse_mesh_files/solid.1.face";
     split_corner_elements=true;
     Orig_mesh_pt = new SolidTetgenMesh<TPVDElement<3,3> >(node_file_name,
                                                        element_file_name,
                                                        face_file_name,
                                                        split_corner_elements);
    }
   else
    {
 
     oomph_info << "Doing fluid\n";
 
     // Output directory
     doc_info.set_directory("RESLT_fluid_with_poisson");
     
     //Create bulk mesh, sub-dividing "corner" elements
     node_file_name="short_coarse_mesh_files/fluid.1.node";
     element_file_name="short_coarse_mesh_files/fluid.1.ele";
     face_file_name="short_coarse_mesh_files/fluid.1.face";
     split_corner_elements=true;
     switch_normal=true;
     Orig_mesh_pt =  new SolidTetgenMesh<TPVDElement<3,3> >(node_file_name,
                                                         element_file_name,
                                                         face_file_name,
                                                         split_corner_elements,
                                                         switch_normal);
    }
   
   
   // Check for inverted elements before quadratic trafo
   std::string filename;
   filename=doc_info.directory();
   filename+="/mesh_before_snap.dat";
   Orig_mesh_pt->output(filename.c_str());
   bool mesh_has_inverted_elements;
   Orig_mesh_pt->check_inverted_elements(mesh_has_inverted_elements);  
   cout << "Before quadratic snapping mesh does ";
   if (!mesh_has_inverted_elements) cout << "not ";
   cout << "have inverted elements. \n";
   
 
   if (do_solid==1) 
    {
   // Snap FSI interface to quadratic surface
     switch_normal=false; 
     Orig_mesh_pt->template snap_to_quadratic_surface<TPVDElement<3,3> >(
      solid_fsi_boundary_id,
      "short_coarse_mesh_files/quadratic_fsi_boundary.dat",
      switch_normal);
 
     // Snap outer solid interface to quadratic surface
     switch_normal=true;
     Orig_mesh_pt->template snap_to_quadratic_surface<TPVDElement<3,3> >(
      solid_outer_boundary_id,
      "short_coarse_mesh_files/quadratic_outer_solid_boundary.dat",
      switch_normal);
    }
   else
    {
     // Snap FSI interface to quadratic surface
     switch_normal=false; 
     Orig_mesh_pt->template snap_to_quadratic_surface<TPVDElement<3,3> >(
      fluid_fsi_boundary_id,
      "short_coarse_mesh_files/quadratic_fsi_boundary.dat",
      switch_normal);
    }
   
   
   // Check
   filename=doc_info.directory();
   filename+="/mesh_after_snap.dat";
   Orig_mesh_pt->output(filename.c_str());
   std::ofstream inverted_mesh_elements;
   filename=doc_info.directory();
   filename+="/inverted_elements_after_snap.dat";
   inverted_mesh_elements.open(filename.c_str());
   Orig_mesh_pt->check_inverted_elements(mesh_has_inverted_elements,
                                         inverted_mesh_elements);  
   inverted_mesh_elements.close(); 
   cout << "After quadratic snapping mesh does ";
   if (!mesh_has_inverted_elements) cout << "not ";
   cout << "have inverted elements. \n";
   
 
   
   // Smooth it: 
   //-----------
      
   // Create set containing the nodes whose displacements are contrained
   // (nodes on quadratic boundaries)
   std::set<Node*> pinned_nodes;
   
   if (do_solid==1) 
    {
     
     // Loop over nodes on the FSI boundary in the solid mesh
     unsigned nbound_pinned=solid_fsi_boundary_id.size();
     for(unsigned i=0;i<nbound_pinned;i++)
      {
       //Get the mesh boundary
       unsigned b = solid_fsi_boundary_id[i];
       unsigned num_nod=Orig_mesh_pt->nboundary_node(b);
       for (unsigned inod=0;inod<num_nod;inod++)
        {
         // Get node
         pinned_nodes.insert(Orig_mesh_pt->boundary_node_pt(b,inod));
        }
      }
     
     // Loop over nodes on outer boundary of solid mesh
     nbound_pinned=solid_outer_boundary_id.size();
     for(unsigned i=0;i<nbound_pinned;i++)
      {
       //Get the mesh boundary
       unsigned b = solid_outer_boundary_id[i];
       unsigned num_nod=Orig_mesh_pt->nboundary_node(b);
       for (unsigned inod=0;inod<num_nod;inod++)
        {
         // Get node     
         pinned_nodes.insert(Orig_mesh_pt->boundary_node_pt(b,inod));
        }
      }
     
    }
   else
    {
     // Loop over nodes on the FSI boundary in the solid mesh
     unsigned nbound_pinned=fluid_fsi_boundary_id.size();
     for(unsigned i=0;i<nbound_pinned;i++)
      {
       //Get the mesh boundary
       unsigned b = fluid_fsi_boundary_id[i];
       unsigned num_nod=Orig_mesh_pt->nboundary_node(b);
       for (unsigned inod=0;inod<num_nod;inod++)
        {
         // Get node
         pinned_nodes.insert(Orig_mesh_pt->boundary_node_pt(b,inod));
        }
      }
    }
     
   // Smooth
   PoissonSmoothMesh<TPoissonElement<3,3> >()(Orig_mesh_pt,pinned_nodes);
 
 
   // Check
   filename=doc_info.directory();
   filename+="/mesh_after_smooth.dat";
   Orig_mesh_pt->output(filename.c_str());
   filename=doc_info.directory();
   filename+="/inverted_elements_after_smooth.dat";
   inverted_mesh_elements.open(filename.c_str());
   Orig_mesh_pt->check_inverted_elements(mesh_has_inverted_elements,
                                         inverted_mesh_elements);  
   inverted_mesh_elements.close(); 
   cout << "After smoothing mesh does ";
   if (!mesh_has_inverted_elements) cout << "not ";
   cout << "have inverted elements. \n";
 
   // Cleanup
   delete Orig_mesh_pt;
   Orig_mesh_pt=0;
   
  }
 
 }

References b, oomph::SolidMesh::boundary_node_pt(), oomph::Mesh::check_inverted_elements(), oomph::DocInfo::directory(), MergeRestartFiles::filename, i, oomph::Mesh::nboundary_node(), oomph::oomph_info, oomph::Mesh::output(), oomph::DocInfo::set_directory(), and oomph::Global_string_for_annotation::string().

Referenced by main().

main()

int main	(	int argc	,
		char **	argv
	)

/////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////// Demonstrate how to update FSI fluid and solid meshes to quadratic representation of FSI interface.

{
 // Store command line arguments
 CommandLineArgs::setup(argc,argv);
 
 // Demonstrate smoothing with nonlinear elastictiy
 demo_smoothing_with_nonlinear_elasticity();
 
 // Demonstrate smoothing with nonlinear elastictiy
 demo_smoothing_with_poisson();
 
 // Demonstrate smoothing with nonlinear elastictiy
 demo_smoothing_with_linear_elasticity();
 
}

References demo_smoothing_with_linear_elasticity(), demo_smoothing_with_nonlinear_elasticity(), demo_smoothing_with_poisson(), and Flag_definition::setup().