xda_to_poly_fsi.cc File Reference

#include <iostream>
#include <sstream>
#include <string>
#include <iterator>
#include "generic.h"

Functions
int	main ()

Function Documentation

main()

int main

(

)

{
 
 // Prompt input from user
 //-----------------------
 std::cout << "Please enter the file name without the file extension '.xda': ";
 std::string input_filename;
 std::cin >> input_filename;
 
 
 std::cout << "\n\nEnter the (uniform) wall thickness ";
 double wall_thickness;
 std::cin >> wall_thickness;
 
 
 std::cout << "\n\nDo you want to create a separate ID for each planar facet\n"
           << "on the fluid-structure interaction boundary?"
           << "\n" <<"Enter y or n: ";
 char multi_boundary_ids_marker;
 std::cin >> multi_boundary_ids_marker;
 bool do_multi_boundary_ids=true; 
 if(multi_boundary_ids_marker=='n') do_multi_boundary_ids=false;
 
 
 // Define filenames
 //-----------------
 
 // Build filename for xda input file
 std::string xda_file_name;
 xda_file_name=input_filename+".xda";
 
 std::string fluid_surface_file;
 fluid_surface_file="fluid_"+input_filename+".poly";
 
 std::string solid_surface_file;
 solid_surface_file="solid_"+input_filename+".poly";
 
 
 // Open and process xda input file
 //--------------------------------
 std::ifstream infile(xda_file_name.c_str(),std::ios_base::in);
 unsigned n_node;
 unsigned n_element;
 unsigned nbound_cond;
 
 // Dummy storage to jump lines
 char dummy[100];
 
 // Ignore file format
 infile.getline(dummy, 100);
 
 // Get number of elements
 infile>>n_element;
 
 // Ignore rest of line
 infile.getline(dummy, 100);
 
 // Get number of nodes
 infile>>n_node;
 
 // Ignore rest of line
 infile.getline(dummy, 100);
 
 // Ignore sum of element weights (whatever that is...)
 infile.getline(dummy, 100);
 
 // Get number of boundaries = number of boundary conditions
 infile>>nbound_cond;
 
 // Keep reading until "Title String"
 while (dummy[0]!= 'T') {infile.getline(dummy, 100);}
 
 
 // Read first line with node labels and count them
 string line;
 getline(infile,line);
 istringstream ostr(line);
 istream_iterator<string> it(ostr);
 istream_iterator<string> end;
 unsigned nnod_el = 0;
 Vector<unsigned> first_node;
 while (it != end)
  {
   first_node.push_back(atoi((*it).c_str()));
   it++;
   nnod_el++;
  }
 oomph_info << "Number of nodes per element: " << nnod_el << std::endl;
 
 // Storage for the global node numbers listed element-by-element
 Vector<unsigned> global_node(n_element*nnod_el);
 
 // Read in nodes
 unsigned k=0;
 
 // Copy across first nodes
 cout << "First nodes: " << std::endl;
 for(unsigned j=0;j<nnod_el;j++)
  {
   global_node[k]=first_node[k];
   cout << first_node[k] << " ";
   k++;
  }
 cout << std::endl;
 
 // Read the other ones
 for(unsigned i=1;i<n_element;i++)
  {
   for(unsigned j=0;j<nnod_el;j++)
    {
     infile >> global_node[k];
     k++;
    }
   infile.getline(dummy, 100);
  }
 
 
 // Create storage for coordinates
 Vector<double> x_node(n_node);
 Vector<double> y_node(n_node);
 Vector<double> z_node(n_node);
 
 // Get nodal coordinates
 for(unsigned i=0;i<n_node;i++)
  {
   infile>>x_node[i];
   infile>>y_node[i];
   infile>>z_node[i];
  }
 
 // Vector of bools, will tell us if we already visited a node
 std::vector<bool> done_for_fluid(n_node,false);
 std::vector<bool> done_for_solid(n_node,false);
 
 // A map, indexed by the old node number, gives the new node number 
 // in the fluid poly file
 std::map<unsigned,unsigned> fluid_node_nmbr;
 
 // Each node in the fluid surface creates two nodes in the solid 
 // surfaces, one in the inner surface and the other in the outer surface.
 // These maps, indexed by the old node number, give the new node number 
 // in the solid poly file.
 std::map<unsigned,unsigned> solid_inner_node_nmbr;
 std::map<unsigned,unsigned> solid_outer_node_nmbr;
 
 // nfluid_face[i] will store the number of facets on boundary i+1.
 // The fluid mesh has four distinct boundaries.
 Vector<unsigned> nfluid_face(4);
  
 // nsolid_linking_faces[i] will store the number of facets on the solid
 // boundaries that connect the inner and outer solid surface. 
 // With the enumeration chosen, nsolid_linking_faces[i] corresponds
 // to boundary i+2. (We don't store the number of facets in boundaries 
 // 1 and 5 because they are the same as the number of fluid facets 
 // on boundary 1.
 Vector<unsigned> nsolid_linking_faces(3); 
 
 // Initialise number of nodes in the solid surface
 unsigned n_solid_node=0;
 
 // Initalise number of nodes in the fluid surface
 unsigned n_fluid_node=0;
 
 // Create storage for outer solid nodes coordinates. Indexed by the old 
 // node number, these maps returns the nodes coordinates.
 std::map<unsigned, double> x_outer_node;
 std::map<unsigned, double> y_outer_node;
 std::map<unsigned, double> z_outer_node;
 std::map<unsigned,Vector<double> > normal_for_outer_wall;
 
 
 // Create storage for fluid faces informations :
 // fluid_faces[i][j] will store a vector containing the three node numbers 
 // of the j-th face in the (i+1)-th boundary
 Vector< Vector<Vector<unsigned> > > 
  fluid_faces(4, Vector<Vector<unsigned> >(nbound_cond,Vector<unsigned>(3) ));
 
 // Create storage for solid facets:
 // As the solid inner and outer facets come from the fluid surface, 
 // we only need to store solid facet information in boundaries 2, 3 and 4. 
 // solid_faces[i][j] will store a vector containing the three node numbers 
 // of the j-th face in the (i+2)-th boundary
 Vector< Vector<Vector<unsigned> > > 
  solid_linking_faces(3, Vector<Vector<unsigned> >(nbound_cond,
                                                   Vector<unsigned>(3) ));
 
 
 // Storage for node numbers of nodes that make up fluid fsi boundary
 // in quadratic (six node triangle) representation)
 std::set<unsigned> quadratic_fsi_boundary_nodes;
 Vector<Vector<unsigned> > quadratic_fsi_boundary_face_element_nodes;
 
 
 // Get boundary enumeration
 //--------------------------
 unsigned element_nmbr;
 unsigned side_nmbr;
 int bound_id;
 for(unsigned i=0;i<nbound_cond;i++)
  {
 
   // Number of the element
   infile>> element_nmbr ;
 
   // Which side/face on the tet are we dealing with (xda enumeratation)?
   infile>> side_nmbr ;
 
   // What's the boundary ID?
   infile>> bound_id ;
   
   // Sanity check
   if(bound_id!=1 && bound_id!=2 && bound_id!=3 && bound_id!=4 )
    {
     std::ostringstream error_stream;
     error_stream << "bound_id=" << bound_id 
                  << ". This function only takes xda type meshes with"
                  <<" one inflow (boundary 2) and two outflows (boundaries 3"
                  <<" and 4); the FSI interface must have the id 1. Your"
                  <<" input file contains a boundary with id : " 
                  << bound_id <<".\n"
                  <<"Don't panic, there are only few things to change"
                  <<"  in this well commented code, good luck ;)  \n";     
     throw OomphLibError(
      error_stream.str(),
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
        
   // Identify the "side nodes" (i.e. the nodes on the faces of
   // the bulk tet) according to the 
   // conventions in '.xda' mesh files so that orientation of the
   // faces is always the same (required for computation of 
   // outer unit normals
   Vector<unsigned> side_node(6);
 
   switch(side_nmbr)
    {
    case 0:
     side_node[0]=global_node[nnod_el*element_nmbr+1];
     side_node[1]=global_node[nnod_el*element_nmbr];
     side_node[2]=global_node[nnod_el*element_nmbr+2];
     side_node[3]=global_node[nnod_el*element_nmbr+4];
     side_node[4]=global_node[nnod_el*element_nmbr+6];
     side_node[5]=global_node[nnod_el*element_nmbr+5];
     break;
       
    case 1:
     side_node[0]=global_node[nnod_el*element_nmbr];
     side_node[1]=global_node[nnod_el*element_nmbr+1];
     side_node[2]=global_node[nnod_el*element_nmbr+3];
     side_node[3]=global_node[nnod_el*element_nmbr+4];
     side_node[4]=global_node[nnod_el*element_nmbr+8];
     side_node[5]=global_node[nnod_el*element_nmbr+7];
     break;
 
    case 2:
     side_node[0]=global_node[nnod_el*element_nmbr+1];
     side_node[1]=global_node[nnod_el*element_nmbr+2];
     side_node[2]=global_node[nnod_el*element_nmbr+3];
     side_node[3]=global_node[nnod_el*element_nmbr+5];
     side_node[4]=global_node[nnod_el*element_nmbr+9];
     side_node[5]=global_node[nnod_el*element_nmbr+8];
     break;
 
    case 3:
     side_node[0]=global_node[nnod_el*element_nmbr+2];
     side_node[1]=global_node[nnod_el*element_nmbr];
     side_node[2]=global_node[nnod_el*element_nmbr+3];
     side_node[3]=global_node[nnod_el*element_nmbr+6];
     side_node[4]=global_node[nnod_el*element_nmbr+7];
     side_node[5]=global_node[nnod_el*element_nmbr+9];
     break;
       
    default :
     throw OomphLibError(
      "Unexcpected side number in your '.xda' input file\n",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
   
   
   // Are we on the FSI interface?
   //-----------------------------
   if(bound_id==1)
    {
     // Create storage for the normal vector to the face
     Vector<double> normal(3,0.0);
     
     // Identify quadratic nodes
     quadratic_fsi_boundary_nodes.insert(side_node[0]);
     quadratic_fsi_boundary_nodes.insert(side_node[1]);
     quadratic_fsi_boundary_nodes.insert(side_node[2]);
     quadratic_fsi_boundary_nodes.insert(side_node[3]);
     quadratic_fsi_boundary_nodes.insert(side_node[4]);
     quadratic_fsi_boundary_nodes.insert(side_node[5]);
 
     // Identify quadratic facet
     Vector<unsigned> quadratic_nodes(6);
     quadratic_nodes[0]=side_node[0];
     quadratic_nodes[1]=side_node[1];
     quadratic_nodes[2]=side_node[2];
     quadratic_nodes[3]=side_node[3];
     quadratic_nodes[4]=side_node[4];
     quadratic_nodes[5]=side_node[5];
     quadratic_fsi_boundary_face_element_nodes.push_back(quadratic_nodes);
 
     // Get the vertex nodes' coordinates
     double x1=x_node[side_node[0] ]; 
     double x2=x_node[side_node[1] ];
     double x3=x_node[side_node[2] ];
 
     double y1=y_node[side_node[0] ];
     double y2=y_node[side_node[1] ];
     double y3=y_node[side_node[2] ];
 
     double z1=z_node[side_node[0] ];
     double z2=z_node[side_node[1] ];
     double z3=z_node[side_node[2] ];
 
     // Compute the outer normal vector 
     normal[0] =(y2-y1)*(z3-z1) - (z2-z1)*(y3-y1);                 
     normal[1] =(z2-z1)*(x3-x1) - (x2-x1)*(z3-z1); 
     normal[2] =(x2-x1)*(y3-y1) - (y2-y1)*(x3-x1);
     
     // Normalise
     double length= sqrt((normal[0])*(normal[0]) + 
                         (normal[1])*(normal[1]) +
                         (normal[2])*(normal[2]));
     for(unsigned idim=0; idim<3; idim++)
      {
       normal[idim]/=length;
      }
     
     // Loop over the facet's nodes vertex/simplex nodes
     for(unsigned ii=0; ii<6; ii++)
      {
       // get the node number
       unsigned inod=side_node[ii];
 
       if (ii<3)
        {
         // Have we already dealt with this node in its incarnation as a
         // fluid node?
         if(!done_for_fluid[inod])
          {
           done_for_fluid[inod]=true;
           
           // This node is a new fluid node
           n_fluid_node++;
          }
        }
 
       // Have we already dealt with this node in its incarnation as a
       // solid node?
       if(!done_for_solid[inod])
        {
         done_for_solid[inod]=true;
         normal_for_outer_wall[inod].resize(3);
         if (ii<3)
          {
           // This node is a solid node on boundary 1 and this node 
           // creates another one on boundary 5
           n_solid_node+=2;
          }
        }
       // Compute the coordinates of the new node (on boundary 5)
       x_outer_node[inod]= x_node[inod];//+ wall_thickness* normal[0];
       y_outer_node[inod]= y_node[inod];//+ wall_thickness* normal[1];
       z_outer_node[inod]= z_node[inod];//+ wall_thickness* normal[2];
       
       normal_for_outer_wall[inod][0]+=normal[0];
       normal_for_outer_wall[inod][1]+=normal[1];
       normal_for_outer_wall[inod][2]+=normal[2];
      }
      
     // We've just created a new fluid facet on boundary 1
     nfluid_face[0]++;
     
     // Store the face node numbers
     for(unsigned ii=0;ii<3; ii++)
      {
       fluid_faces[0][nfluid_face[0]-1][ii]=side_node[ii];
      }
    }
 
   // We're on a non-FSI boundary
   //----------------------------
   else
    {
     // loop over the facet's vertex/simplex nodes
     for(unsigned ii=0; ii<3; ii++)
      {
       // get the node number
       unsigned inod=side_node[ii];
       
       // Have we dealt with this fluid node?
       if(!done_for_fluid[inod])
        {
         done_for_fluid[inod]=true;
         
         // this node is a new fluid node
         n_fluid_node++;
        }
      }
     // We've just created a new facet on boundary bound_id
     nfluid_face[bound_id-1]++;
     
     // Store the face node numbers
     for(unsigned ii=0;ii<3; ii++)
      {
       fluid_faces[bound_id-1][nfluid_face[bound_id-1]-1][ii]=side_node[ii];
      }
    }
  }
 
 // Done reading xda file
 infile.close();
 
 
 // Now create storage for fluid facets on the four boundaries of the
 // fluid mesh
 for(unsigned i=0; i<4; i++)
  {
   fluid_faces[i].resize(nfluid_face[i]);
  }
 
 
 // Re-normalise unit normals to wall
 for (std::map<unsigned,Vector<double> >::iterator it=
       normal_for_outer_wall.begin();
      it!=normal_for_outer_wall.end();it++)
  {
   double sum=0.0;
   for (unsigned i=0;i<3;i++)
    {
     sum+=pow((*it).second[i],2);
    }
   for (unsigned i=0;i<3;i++)
    {
     (*it).second[i]/=sqrt(sum);
    }
  }
 
  // Output quadratic representation of fsi boundary (six node triangles)
 //---------------------------------------------------------------------
 unsigned nquadratic_nodes=quadratic_fsi_boundary_nodes.size();
 
 std::string filename;
 filename=input_filename+"_quadratic_fsi_boundary.dat";
 std::ofstream quadratic_fsi_boundary_file(filename.c_str());
 quadratic_fsi_boundary_file << nquadratic_nodes 
                             << " # number of quadratic FSI nodes\n";
 quadratic_fsi_boundary_file << quadratic_fsi_boundary_face_element_nodes.size()
                             << " # number of FSI facets\n";
 for (std::set<unsigned>::iterator it=quadratic_fsi_boundary_nodes.begin();
      it!=quadratic_fsi_boundary_nodes.end();it++)
  {
   quadratic_fsi_boundary_file <<  x_node[*it] << " "
                               <<  y_node[*it] << " "
                               <<  z_node[*it] << " "
                               << *it << std::endl;
  }
 
 
 for (Vector<Vector<unsigned> >::iterator it=
       quadratic_fsi_boundary_face_element_nodes.begin();
      it!=quadratic_fsi_boundary_face_element_nodes.end();it++)
  {
   quadratic_fsi_boundary_file << (*it)[0] << " "
                               << (*it)[1] << " "
                               << (*it)[2] << " "
                               << (*it)[3] << " "
                               << (*it)[4] << " "
                               << (*it)[5] << " " << std::endl;
  } 
 quadratic_fsi_boundary_file.close();
 
 
 
 // Output quadratic representation of outer solid boundary (six node triangles)
 //-----------------------------------------------------------------------------
 filename=input_filename+"_quadratic_outer_solid_boundary.dat";
 std::ofstream quadratic_outer_solid_boundary_file(filename.c_str());
 quadratic_outer_solid_boundary_file 
  << nquadratic_nodes 
  << " # number of quadratic nodes on outer solid boundary\n";
 quadratic_outer_solid_boundary_file 
  << quadratic_fsi_boundary_face_element_nodes.size()
  << " # number of facets on outer solid boundary\n";
 for (std::set<unsigned>::iterator it=
       quadratic_fsi_boundary_nodes.begin();
      it!=quadratic_fsi_boundary_nodes.end();it++)
  {
   quadratic_outer_solid_boundary_file 
    <<  x_outer_node[*it]+wall_thickness*normal_for_outer_wall[*it][0] << " "
    <<  y_outer_node[*it]+wall_thickness*normal_for_outer_wall[*it][1] << " "
    <<  z_outer_node[*it]+wall_thickness*normal_for_outer_wall[*it][2] << " "
    << *it << std::endl;
  }
 
 
 for (Vector<Vector<unsigned> >::iterator it=
       quadratic_fsi_boundary_face_element_nodes.begin();
      it!=quadratic_fsi_boundary_face_element_nodes.end();it++)
  {
   quadratic_outer_solid_boundary_file << (*it)[0] << " "
                                       << (*it)[1] << " "
                                       << (*it)[2] << " "
                                       << (*it)[3] << " "
                                       << (*it)[4] << " "
                                       << (*it)[5] << " " << std::endl;
  } 
 quadratic_outer_solid_boundary_file.close();
 
 
 
 
 // Start poly files, describing fluid and solid domains via
 //---------------------------------------------------------
 // planar facets
 //--------------
 std::ofstream fluid_output_stream(fluid_surface_file.c_str(),std::ios::out);
 fluid_output_stream << "# This poly file was created from a"
                     <<"  VMTK mesh." << '\n'
                     << "# VMTK assigns for the front inflow "
                     <<"face the boundary id 2, the  left " << '\n'
                     <<"# bifurcation face the id 3, the right"
                     <<" bifurcation face the id 4 and for the other"
                     <<" boundaries (FSI interface) the id 1 " << '\n'
                     << "# Oomph-lib's meshes need distinct ids  "
                     <<"for each planar facet, so we assign new boundary ids."
                     << '\n'
                     <<"# The relation between old and new ids is listed"
                     <<" at the end of this file.\n";
 
 std::ofstream solid_output_stream(solid_surface_file.c_str(),std::ios::out);
 solid_output_stream << "# This poly file was created from the fluid vmtk "
                     << " mesh assuming constant wall thickness.\n"
                     << "# Oomph-lib's meshes need distinct   "
                     <<"ids for each planar facet so we assign new \n"
                     <<"# boundary ids. The relation between old and new "
                     << " ids is listed at the end of this file.\n";
 
 
 // Write the node list for the fluid boundaries
 //---------------------------------------------
 fluid_output_stream << '\n' << "# Node list : " << '\n';
 
 // Number of nodes, dimension =3 , no attributes, with boundary markers
 fluid_output_stream <<  n_fluid_node <<  ' ' << 3 <<  ' ' << 0 
                     <<  ' ' << 1 << '\n'<< '\n';
 
 // Has a specific node been done yet?
 std::vector<bool> done_fluid_node(n_node,false);
 unsigned counter=0;
 
 // Loop over the boundaries in the fluid mesh
 for(unsigned i=0; i<4; i++)
  {
   // How many facets are on boundary 'i+1' ?
   unsigned nface=nfluid_face[i];
 
   // Loop over the facets
   for(unsigned iface=0; iface<nface; iface++)
    {
     // Get pointer to the Vector storing the three vertex nodes defining
     // the facet
     Vector<unsigned>* face_node=&(fluid_faces[i][iface]);
 
     // Loop over the facet's nodes
     for(unsigned ii=0; ii<3; ii++)
      {
       // Get the node number (in the old numbering scheme)
       unsigned inod=(*face_node)[ii];
 
       // Done yet?
       if(!done_fluid_node[inod])
        {
         done_fluid_node[inod]=true;
         counter++;
         
         // Assign the new node number on the fluid surface
         fluid_node_nmbr[inod]=counter;
 
         // Write the node coordinates (boundary marker =0 )
         fluid_output_stream <<  counter <<  ' ' 
                             << x_node[inod]  <<  ' ' 
                             << y_node[inod]  <<  ' ' 
                             << z_node[inod] << ' '
                             << 0 << '\n'; 
        }
      } 
    }
  }
 
 
 // Write the node list for the solid boundaries
 //---------------------------------------------
 solid_output_stream << '\n' << "# Node list : " << '\n';
 
 // Number of nodes, dimension =3 , no attributes, with boundary markers
 solid_output_stream <<  n_solid_node <<  ' ' << 3 <<  ' ' << 0 
                     <<  ' ' << 1 << '\n'<< '\n';
 
 // Has a specific node been done yet?
 std::vector<bool> done_solid_node(n_node,false);
 counter=0;
 
 
 // Do FSI and outer surfaces
 //--------------------------
 
 // How many facets are on boundary 0 (FSI surface)
 unsigned nface=nfluid_face[0];
 
 // Loop over the facets
 for(unsigned iface=0; iface<nface; iface++)
  {
   // Get pointer to the Vector storing the three nodes
   Vector<unsigned>* face_node=&(fluid_faces[0][iface]);
 
   // Loop over the facet's nodes
   for(unsigned ii=0; ii<3; ii++)
    {
     // get the node number (in the old numbering scheme)
     unsigned inod=(*face_node)[ii];
 
     // Done yet?
     if(!done_solid_node[inod])
      {
       done_solid_node[inod]=true;
       counter++;
       
       // Assign the new node number on the solid surface
       solid_inner_node_nmbr[inod]=counter;
       
       solid_output_stream <<  counter <<  ' ' 
                           << x_node[inod]  <<  ' ' 
                           << y_node[inod]  <<  ' ' 
                           << z_node[inod] << ' '
                           << 0 << '\n'; 
 
       // Store the new nodes on boundary 5 (outer surface)
       counter++;
       
       // Assign the new node number on the solid surface
       solid_outer_node_nmbr[inod]=counter;
       
       solid_output_stream 
        <<  counter <<  " "
        << x_outer_node[inod]+wall_thickness*normal_for_outer_wall[inod][0]  
        <<  " "
        << y_outer_node[inod]+wall_thickness*normal_for_outer_wall[inod][1]
        << " " 
        << z_outer_node[inod]+wall_thickness*normal_for_outer_wall[inod][2]  
        << 0   << '\n';
      }   
    } 
  }
 
 // Face lists for the fluid mesh 
 //------------------------------
 fluid_output_stream << '\n' << "# Face list : " << '\n';
 
 // Total number of facets and, yes, we do specify boundary IDs
 fluid_output_stream <<  nfluid_face[0]+nfluid_face[1]+nfluid_face[2]
  +nfluid_face[3] <<  ' ' << 1 << '\n'<< '\n'; 
 
 
 // All four original boundaries
 counter=0;
 for(unsigned i=0; i<4; i++)
  {
   fluid_output_stream <<'\n'
                       << "#============================="
                       <<"====================================="
                       << '\n'<<"# Faces originally on boundary " 
                       << i+1 << '\n'
                       <<"# --------------------------------" << '\n';
 
   // How many facest are on boundary 'i+1' ?
   unsigned nface=nfluid_face[i];
 
   // Loop over the all facets
   for(unsigned iface=0; iface<nface; iface++)
    {
     counter++;
     fluid_output_stream <<"# Face " << counter << '\n'; 
     
     // one polygon, zero holes
     fluid_output_stream << 1 <<  ' ' << 0 << ' ' ;
 
     // If we want a separate ID for each facet, the boundary ID is 
     // counter, otherwise, the boundary ID is i+1
     if(do_multi_boundary_ids)
      {
       fluid_output_stream << counter << '\n';
      }
     else
      {
       fluid_output_stream << i+1 << '\n';
      }
     
     // We have three vertices
     fluid_output_stream << 3 <<  ' ' ;
       
     // Get pointer to the Vector storing the three nodes
     Vector<unsigned>* face_node=&(fluid_faces[i][iface]);
     
     // This vector stores the nodes that are on boundaries 1 AND
     // another one (2,3 or 4). We only store these nodes if a
     // facet contains two or more such "double boundary" nodes.
     // These are used to identify nodes that are part of the
     // linking faces that connect the inner and outer solid surfaces.
     Vector<unsigned> double_boundary_nod(3);
 
     // Storage for the number of entries in double_boundary_nod;
     unsigned n_double_boundary_nodes=0;
     
     // loop over the three vertex/simplex nodes
     for(unsigned l=0; l<3;l++)
      {
       
       // get the old node number
       unsigned inod=(*face_node)[l];
       
       // Write the vertices indices 
       fluid_output_stream <<fluid_node_nmbr[inod] << ' ';
       
       // find out how many nodes are double boundary nodes
       if(i!=0)
        {
         if (done_for_solid[inod])
          {
           double_boundary_nod[n_double_boundary_nodes]=inod;
           n_double_boundary_nodes++;
          }
        }
      }
       
     // If we have more than one double boundary node,
     // create faces linking the solid faces on
     // boundary 1 and the solid faces on boundary 5
     if(n_double_boundary_nodes>1) 
      {
       for(unsigned idbn=0; idbn<n_double_boundary_nodes-1; idbn++)
        {
         // Each two double boundary nodes create two faces
         for(unsigned j=0; j<2; j++)
          {    
           unsigned index=0;
 
           for(unsigned l=j+idbn; l<2+idbn; l++)
            {
             // *[i-1] because *[k] stores informations for the 
             // (k+2)-th boundary
             solid_linking_faces[i-1][nsolid_linking_faces[i-1]][index]= 
              solid_inner_node_nmbr[double_boundary_nod[l]];
             index++;
            }
           for(unsigned l=idbn; l<idbn+j+1; l++)
            {
             // *[i-1] because *[k] stores informations for the 
             // (k+2)-th boundary
             solid_linking_faces[i-1][nsolid_linking_faces[i-1]][index]=
              solid_outer_node_nmbr[double_boundary_nod[l]];
             index++;
            }
           nsolid_linking_faces[i-1]++;
          }           
        }
      }
     fluid_output_stream << '\n'<< '\n';
    }
  } 
 
 
 // Write the face list of the solid surface
 //-----------------------------------------
 solid_output_stream << '\n' << "# Face list : " << '\n';
 
 // Total number of facets and, yes, we do specify boundary IDs
 solid_output_stream <<  2*nfluid_face[0]+ nsolid_linking_faces[0]
  + nsolid_linking_faces[1]+ nsolid_linking_faces[2]  
                     <<  ' ' << 1 << '\n'<< '\n';
  
 // Loop over all five original boundaries in the fluid mesh
 counter=0;
 for(unsigned i=0; i<5; i++)
  {
   solid_output_stream <<'\n'
                       << "#====================================="
                       <<"============================="
                       << '\n'<<"# Faces originally on boundary "
                       <<i+1 << '\n'
                       <<"# --------------------------------" << '\n';
 
   // Get the number of facets -- either the number of faces on the
   // fluid's FSI interface or the number of "linking" faces that
   // connect the inner and outer solid surfaces
   unsigned nface;
   if(i==0 || i==4)
    {
     nface=nfluid_face[0];
    }
   else
    {
     nface=nsolid_linking_faces[i-1];
    }
 
   // Loop over faces
   for(unsigned iface=0; iface<nface; iface++)
    {
     // Get pointer to the Vector storing the three vertex/simplex nodes
     Vector<unsigned>* face_node;
     if(i==0 || i==4)
      {
       face_node=&(fluid_faces[0][iface]);
      }
     // *[i-1] because *[k] stores informations for the (k+2)-th boundary
     else
      {
       face_node=&(solid_linking_faces[i-1][iface]);
      }
       
     counter++;
     solid_output_stream <<"# Face " <<  counter << '\n';  
      
     // One polygon, zero holes, boundary  counter
     solid_output_stream << 1 <<  ' ' << 0 << ' ' ;
     
     // If we want a separate ID for each FSI facet, the boundary ID is 
     // counter, otherwise, the boundary ID is i+1
     if(do_multi_boundary_ids)
      {
       solid_output_stream << counter << '\n';
      }
     else
      {
       solid_output_stream << i+1 << '\n'; 
      }
 
     // Three vertices and their indices in node list
     solid_output_stream << 3 <<  ' ' ;
     
     // Loop over the three vertex/simplex nodes
     for(unsigned l=0; l<3;l++)
      {
       // get the node number
       unsigned inod=(*face_node)[l];
       
       // FSI surface
       if(i==0)
        {
         // it's the old node number
         solid_output_stream << solid_inner_node_nmbr[inod] << ' ';
        }
       // Outer surface
       else if(i==4)
        {
         // it's the old node number
         solid_output_stream << solid_outer_node_nmbr[inod] << ' ';
        }
       // Linking surfaces
       else
        {
         // it's the new node number
         solid_output_stream <<inod << ' ';
        }
       
      } 
     solid_output_stream << '\n'<< '\n';
    }
  }
 
 
 // Write the hole and region lists of the fluid and solid (empty)
 //---------------------------------------------------------------
 fluid_output_stream << '\n' << "# Hole list : " << '\n';
 fluid_output_stream << 0 << '\n';
 fluid_output_stream << '\n' << "# Region list : " << '\n';
 fluid_output_stream << 0 << '\n';
 
 solid_output_stream << '\n' << "# Hole list : " << '\n';
 solid_output_stream << 0 << '\n';
 solid_output_stream << '\n' << "# Region list : " << '\n';
 solid_output_stream << 0 << '\n';
 
   
 
 
 
 // Doc enumeration of boundaries
 //------------------------------
 filename=input_filename+"_boundary_enumeration.dat";
 std::ofstream boundary_enumeration(filename.c_str());
 
 // Write the boundary informations for the fluid surface
 fluid_output_stream << '\n'<< '\n'<< '\n' 
                     <<"# Old (vmtk) enumeration of fluid boundaries:\n"
                     <<"# 1: FSI surface\n"
                     <<"# 2: Inlet\n"
                     <<"# 3: Outlet\n"
                     <<"# 4: Outlet\n"
                     <<"# The new boundary ids are as follow:"<< '\n' << '\n'; 
 
 
 // Loop over four initial fluid boundaries
 for(unsigned i=0; i<4; i++)
  {
   fluid_output_stream << "# Boundary "<< i+1
                       << " : from boundary id " ;
   unsigned id=0;
   for(unsigned j=0; j<i; j++)
    {
     id+=nfluid_face[j];
    }
   
   fluid_output_stream << id << " until boundary id " 
                       << id+nfluid_face[i]-1 << '\n';
   boundary_enumeration << id  << " " 
                        << id+nfluid_face[i]-1  << " " 
                        << "# (fluid boundary" << i+1 << ") \n";
  }
 fluid_output_stream.close();
 
 
 // Write the boundary informations for the solid surface
 solid_output_stream << '\n'<< '\n'<< '\n' 
                     <<"# The new boundary ids are as follow:"<< '\n' << '\n';
 
 solid_output_stream << "# the inner surface : from boundary id "  
                     << 0 << " until boundary id "; 
 unsigned id=nfluid_face[0]-1;
 solid_output_stream << id << "\n";
 boundary_enumeration << 0 << " " << id  << " # (solid fsi boundary)\n"; 
 
 
 solid_output_stream << "# the front inflow face : from boundary id "  
                     << id+1 << " until boundary id " ;
 boundary_enumeration << id+1 << " ";
 id+=nsolid_linking_faces[0];
 solid_output_stream <<id<< "\n";
 boundary_enumeration << id << " # (solid inflow boundary)\n"; 
 
 
 solid_output_stream << "# the left bifurcation face : from boundary id "  
                     << id+1 << " until boundary id " ;
 boundary_enumeration << id+1 << " ";
 id+=nsolid_linking_faces[1];
 solid_output_stream << id<< "\n" ;
 boundary_enumeration << id << " # (solid left outflow boundary)\n"; 
 
 solid_output_stream << "# the right bifurcation : from boundary id "  
                     << id+1 << " until boundary id " ;
 boundary_enumeration << id+1 << " ";
 id+=nsolid_linking_faces[2];
 solid_output_stream << id << "\n";
 boundary_enumeration << id << " # (solid right outflow boundary)\n"; 
 
 solid_output_stream << "# the outer surface : from boundary id "  
                     << id+1 << " until boundary id " ;
 boundary_enumeration << id+1 << " ";
 id+=nfluid_face[0];
 solid_output_stream << id << "\n";
 boundary_enumeration << id << " # (solid outer boundary)\n"; 
 solid_output_stream.close();
 boundary_enumeration.close();
 
 
} 

References calibrate::default, Eigen::placeholders::end, MergeRestartFiles::filename, i, j, k, calibrate::line, WallFunction::normal(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, out(), Eigen::bfloat16_impl::pow(), sqrt(), oomph::Global_string_for_annotation::string(), Global_parameters::x1(), and Global_parameters::x2().