channel_spine_mesh.template.cc

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
#ifndef OOMPH_CHANNEL_SPINE_MESH_TEMPLATE_CC
#define OOMPH_CHANNEL_SPINE_MESH_TEMPLATE_CC

#include "channel_spine_mesh.template.h"
#include "rectangular_quadmesh.template.cc"


namespace oomph
{



//===========================================================================
/// Constructor for spine 2D mesh: Pass number of elements in x-direction, 
/// number of elements in y-direction, axial length and height of layer, 
/// and pointer to timestepper (defaults to Static timestepper).
///
/// The mesh contains a layer of spinified fluid elements (of type ELEMENT;
/// e.g  SpineElement<QCrouzeixRaviartElement<2>)
/// and a surface layer of corresponding Spine interface elements
/// of type INTERFACE_ELEMENT, e.g.
/// SpineLineFluidInterfaceElement<ELEMENT> for 2D planar
/// problems.
//===========================================================================
template<class ELEMENT>
ChannelSpineMesh<ELEMENT>::ChannelSpineMesh(
 const unsigned &nx0,
 const unsigned &nx1,
 const unsigned &nx2,
 const unsigned &ny,
 const double &lx0,
 const double &lx1,
 const double &lx2,
 const double &h,
 GeomObject* wall_pt, 
 TimeStepper* time_stepper_pt) :  
 RectangularQuadMesh<ELEMENT >(nx0+nx1+nx2,ny,0.0,lx0+lx1+lx2,
                               0.0,h,false,false,time_stepper_pt),
 Nx0(nx0), Nx1(nx1), Nx2(nx2), 
 Lx0(lx0), Lx1(lx1), Lx2(lx2), 
 Wall_pt(wall_pt)
{

 // Mesh can only be built with 2D Qelements.
 MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);

 //Mesh can only be built with spine elements
 MeshChecker::assert_geometric_element<SpineFiniteElement,ELEMENT>(2);

 // We've called the "generic" constructor for the RectangularQuadMesh
 // which doesn't do much...

 // Build the straight line object
 Straight_wall_pt = new StraightLine(h);

 // Now build the mesh: 
 build_channel_spine_mesh(time_stepper_pt);

}



//===========================================================================
/// Constuctor for spine 2D mesh: Pass number of elements in x-direction, 
/// number of elements in y-direction, axial length and height of layer, 
/// a boolean flag to make the mesh periodic in the x-direction,
/// and pointer to timestepper (defaults to Static timestepper).
///
/// The mesh contains a layer of elements (of type ELEMENT)
/// and a surface layer of corresponding Spine interface elements (of type
/// SpineLineFluidInterfaceElement<ELEMENT>).
//===========================================================================
template<class ELEMENT>
ChannelSpineMesh<ELEMENT>::ChannelSpineMesh(
 const unsigned &nx0,
 const unsigned &nx1,
 const unsigned &nx2,
 const unsigned &ny,
 const double &lx0,
 const double &lx1,
 const double &lx2,
 const double &h,
 GeomObject* wall_pt, 
 const bool& periodic_in_x,
 TimeStepper* time_stepper_pt) :  
 RectangularQuadMesh<ELEMENT >(nx0+nx1+nx2,ny,0.0,lx0+lx1+lx2,0.0,h,
                               periodic_in_x,false,time_stepper_pt),
 Nx0(nx0), Nx1(nx1), Nx2(nx2), 
 Lx0(lx0), Lx1(lx1), Lx2(lx2), 
 Wall_pt(wall_pt)
{
 // Mesh can only be built with 2D Qelements.
 MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);

 //Mesh can only be built with spine elements
 MeshChecker::assert_geometric_element<SpineFiniteElement,ELEMENT>(2);

 // We've called the "generic" constructor for the RectangularQuadMesh
 // which doesn't do much...

 // Build the straight line object
 Straight_wall_pt = new StraightLine(h);
 
 // Now build the mesh: 
 build_channel_spine_mesh(time_stepper_pt);

}





//===========================================================================
/// Helper function that actually builds the channel-spine mesh
/// based on the parameters set in the various constructors
//===========================================================================
template<class ELEMENT>
void ChannelSpineMesh<ELEMENT>::build_channel_spine_mesh(
 TimeStepper* time_stepper_pt) 
{
 // Build the underlying quad mesh: 
 RectangularQuadMesh<ELEMENT >::build_mesh(time_stepper_pt);

 // Read out the number of elements in the x-direction and y-direction
 // and in each of the left, centre and right regions
 unsigned n_x = this->Nx;
 unsigned n_y = this->Ny;
 unsigned n_x0 = this->Nx0;
 unsigned n_x1 = this->Nx1;
 unsigned n_x2 = this->Nx2;

 // Set up the pointers to elements in the left region
 unsigned nleft=n_x0*n_y;;
 Left_element_pt.reserve(nleft);
 unsigned ncentre=n_x1*n_y;;
 Centre_element_pt.reserve(ncentre);
 unsigned nright=n_x2*n_y;;
 Right_element_pt.reserve(nright); 
 for(unsigned irow=0;irow<n_y;irow++)
  {
   for (unsigned e=0;e<n_x0;e++)
    {
     Left_element_pt.push_back(
      this->finite_element_pt(irow*n_x+e));
    }
   for (unsigned e=0;e<n_x1;e++)
    {
     Centre_element_pt.push_back(
      this->finite_element_pt(irow*n_x+(n_x0+e)));
    }
   for (unsigned e=0;e<n_x2;e++)
    {
     Right_element_pt.push_back(
      this->finite_element_pt(irow*n_x+(n_x0+n_x1+e)));
    }
  }
 
#ifdef PARANOID
 // Check that we have the correct number of elements
 if (nelement()!=nleft+ncentre+nright)
  {
   throw OomphLibError("Incorrect number of element pointers!",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif

 // Allocate memory for the spines and fractions along spines
 //---------------------------------------------------------

 // Read out number of linear points in the element
 unsigned n_p = dynamic_cast<ELEMENT*>(finite_element_pt(0))->nnode_1d();
 
 unsigned nspine;
 // Allocate store for the spines:
 if (this->Xperiodic)
  {
   nspine = (n_p-1)*n_x;
   Spine_pt.reserve(nspine);
   // Number of spines in each region
   // NOTE that boundary spines are in both regions
   Nleft_spine = (n_p-1)*n_x0+1;
   Ncentre_spine = (n_p-1)*n_x1+1;
   Nright_spine = (n_p-1)*n_x2;
  }
 else
  {
   nspine = (n_p-1)*n_x+1;
   Spine_pt.reserve(nspine);
   // Number of spines in each region
   // NOTE that boundary spines are in both regions
   Nleft_spine = (n_p-1)*n_x0+1;
   Ncentre_spine = (n_p-1)*n_x1+1;
   Nright_spine = (n_p-1)*n_x2+1;
  }

 // end Allocating memory


 // set up the vectors of geometric data & objects for building spines
 Vector<double> r_wall(2), zeta(1), s_wall(1);
 GeomObject* geometric_object_pt=0;

 // LEFT REGION
 // ===========

 // SPINES IN LEFT REGION
 // ---------------------

 // Set up zeta increments
 double zeta_lo=0.0;
 double dzeta=Lx0/n_x0;

 // Initialise number of elements in previous regions:
 unsigned n_prev_elements=0;


 //FIRST SPINE
 //-----------

 //Element 0
 //Node 0
 //Assign the new spine with unit length
 Spine* new_spine_pt=new Spine(1.0);
 new_spine_pt->spine_height_pt()->pin(0);
 Spine_pt.push_back(new_spine_pt);

 // Get pointer to node
 SpineNode* nod_pt=element_node_pt(0,0);
 //Set the pointer to the spine
 nod_pt->spine_pt() = new_spine_pt;
 //Set the fraction
 nod_pt->fraction() = 0.0;
 // Pointer to the mesh that implements the update fct
 nod_pt->spine_mesh_pt() = this; 
 // Set update fct id
 nod_pt->node_update_fct_id()=0;

 // Provide spine with additional storage for wall coordinate 
 // and wall geom object:
 
 {
  // Get the Lagrangian coordinate in the Lower Wall
  zeta[0] =  0.0;
  //Get the geometric object and local coordinate
  Straight_wall_pt->locate_zeta(zeta,geometric_object_pt,s_wall);
  
  //The local coordinate is a geometric parameter
  //This needs to be set (rather than added) because the
  //same spine may be visited more than once
  Vector<double> parameters(1,s_wall[0]);
  nod_pt->spine_pt()->set_geom_parameter(parameters);
  
  // Get position of wall
  Straight_wall_pt->position(s_wall,r_wall);
  
  // Adjust spine height
  nod_pt->spine_pt()->height()=r_wall[1];
  
  // The sub geom object is one (and only) geom object
  // for spine:
  Vector<GeomObject*> geom_object_pt(1);
  geom_object_pt[0] = geometric_object_pt;
  
  // Pass geom object(s) to spine
  nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
 }

 //Loop vertically along the spine
 //Loop over the elements 
 for(unsigned long i=0;i<n_y;i++)
  {
   //Loop over the vertical nodes, apart from the first
   for(unsigned l1=1;l1<n_p;l1++)
    {
     // Get pointer to node
     SpineNode* nod_pt=element_node_pt(i*n_x,l1*n_p);
     //Set the pointer to the spine
     nod_pt->spine_pt() = new_spine_pt;
     //Set the fraction
     nod_pt->fraction()=(double(i)+double(l1)/double(n_p-1))/double(n_y);
     // Pointer to the mesh that implements the update fct
     nod_pt->spine_mesh_pt() = this; 
     // Set update fct id
     nod_pt->node_update_fct_id()=0;
    }
  } // end loop over elements


 //LOOP OVER OTHER SPINES
 //----------------------

 //Now loop over the elements horizontally
 for(unsigned long j=0;j<n_x0;j++)
  {
   //Loop over the nodes in the elements horizontally, ignoring 
   //the first column
   unsigned n_pmax=n_p;
   for(unsigned l2=1;l2<n_pmax;l2++)
    {
     //Assign the new spine with unit height
     new_spine_pt=new Spine(1.0);
     new_spine_pt->spine_height_pt()->pin(0);
     Spine_pt.push_back(new_spine_pt);
     
     // Get the node
     SpineNode* nod_pt=element_node_pt(j,l2);
     //Set the pointer to spine
     nod_pt->spine_pt() = new_spine_pt;
     //Set the fraction
     nod_pt->fraction() = 0.0;
     // Pointer to the mesh that implements the update fct
     nod_pt->spine_mesh_pt() = this; 
     // Set update fct id
     nod_pt->node_update_fct_id()=0;

     {
      // Provide spine with additional storage for wall coordinate 
      // and wall geom object:
      
      // Get the Lagrangian coordinate in the Lower Wall
      zeta[0] =  zeta_lo + double(j)*dzeta + double(l2)*dzeta/(n_p-1.0);
      //Get the geometric object and local coordinate
      Straight_wall_pt->locate_zeta(zeta,geometric_object_pt,s_wall);
      
      //The local coordinate is a geometric parameter
      //This needs to be set (rather than added) because the
      //same spine may be visited more than once
      Vector<double> parameters(1,s_wall[0]);
      nod_pt->spine_pt()->set_geom_parameter(parameters);
      
      // Get position of wall
      Straight_wall_pt->position(s_wall,r_wall);
     
      // Adjust spine height
      nod_pt->spine_pt()->height()=r_wall[1];
      
      // The sub geom object is one (and only) geom object
      // for spine:
      Vector<GeomObject*> geom_object_pt(1);
      geom_object_pt[0] = geometric_object_pt;
      
      // Pass geom object(s) to spine
      nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
     }

     //Loop vertically along the spine
     //Loop over the elements 
     for(unsigned long i=0;i<n_y;i++)
      {
       //Loop over the vertical nodes, apart from the first
       for(unsigned l1=1;l1<n_p;l1++)
        {
         // Get the node
         SpineNode* nod_pt=element_node_pt(i*n_x+j,l1*n_p+l2);
         //Set the pointer to the spine
         nod_pt->spine_pt() = new_spine_pt;
         //Set the fraction
         nod_pt->fraction()=(double(i)+double(l1)/double(n_p-1))/double(n_y);
         // Pointer to the mesh that implements the update fct
         nod_pt->spine_mesh_pt() = this; 
         // Set update fct id
         nod_pt->node_update_fct_id()=0;
        }  
      }
    }
  }

 // Increment number of previous elements
 n_prev_elements+=n_x0*n_y;


 // CENTRE REGION
 // ===========

 zeta_lo=Lx0;
 dzeta=Lx1/n_x1;

 // SPINES IN LEFT REGION
 // ---------------------

 //LOOP OVER OTHER SPINES
 //----------------------

 //Now loop over the elements horizontally
 for(unsigned long j=n_x0;j<n_x0+n_x1;j++)
  {
   //Loop over the nodes in the elements horizontally, ignoring 
   //the first column
   unsigned n_pmax=n_p;
   for(unsigned l2=1;l2<n_pmax;l2++)
    {
     //Assign the new spine with unit height
     new_spine_pt=new Spine(1.0);
     new_spine_pt->spine_height_pt()->pin(0);
     Spine_pt.push_back(new_spine_pt);
     
     // Get the node
     SpineNode* nod_pt=element_node_pt(j,l2);
     //Set the pointer to spine
     nod_pt->spine_pt() = new_spine_pt;
     //Set the fraction
     nod_pt->fraction() = 0.0;
     // Pointer to the mesh that implements the update fct
     nod_pt->spine_mesh_pt() = this; 
     // Set update fct id
     nod_pt->node_update_fct_id()=1;

     {
      // Provide spine with additional storage for wall coordinate 
      // and wall geom object:
      
      // Get the Lagrangian coordinate in the Lower Wall
      zeta[0] =  zeta_lo + double(j-n_x0)*dzeta + double(l2)*dzeta/(n_p-1.0);
      //Get the geometric object and local coordinate
      Wall_pt->locate_zeta(zeta,geometric_object_pt,s_wall);
      
      //The local coordinate is a geometric parameter
      //This needs to be set (rather than added) because the
      //same spine may be visited more than once
      Vector<double> parameters(1,s_wall[0]);
      nod_pt->spine_pt()->set_geom_parameter(parameters);
      
      // Get position of wall
      Wall_pt->position(s_wall,r_wall);
      
      // Adjust spine height
      nod_pt->spine_pt()->height()=r_wall[1];
      
      // The sub geom object is one (and only) geom object
      // for spine:
      Vector<GeomObject*> geom_object_pt(1);
      geom_object_pt[0] = geometric_object_pt;
      
      // Pass geom object(s) to spine
      nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
     }

     //Loop vertically along the spine
     //Loop over the elements 
     for(unsigned long i=0;i<n_y;i++)
      {
       //Loop over the vertical nodes, apart from the first
       for(unsigned l1=1;l1<n_p;l1++)
        {
         // Get the node
         SpineNode* nod_pt=element_node_pt(i*n_x+j,l1*n_p+l2);
         //Set the pointer to the spine
         nod_pt->spine_pt() = new_spine_pt;
         //Set the fraction
         nod_pt->fraction()=(double(i)+double(l1)/double(n_p-1))/double(n_y);
         // Pointer to the mesh that implements the update fct
         nod_pt->spine_mesh_pt() = this; 
         // Set update fct id
         nod_pt->node_update_fct_id()=1;
        }  
      }
    }
  }

 // Increment number of previous elements
 n_prev_elements+=n_x1*n_y;


 // RIGHT REGION
 // ===========

 // SPINES IN RIGHT REGION
 // ----------------------

 // Set up zeta increments
 zeta_lo=Lx0+Lx1;
 dzeta=Lx2/n_x2;

 //LOOP OVER OTHER SPINES
 //----------------------

 //Now loop over the elements horizontally
 for(unsigned long j=n_x0+n_x1;j<n_x0+n_x1+n_x2;j++)
  {
   //Need to copy last spine in previous element over to first spine 
   //in next elements, for all elements other than the first
   if(j>0)
    {
     // For first spine, re-assign the geometric objects, since 
     // we treat it as part of the right region.
     if (j==n_x0+n_x1)
      {
       SpineNode* nod_pt=element_node_pt(j,0);
       // Set update fct id
       nod_pt->node_update_fct_id()=2;
       {
        // Provide spine with additional storage for wall coordinate 
        // and wall geom object:
      
        // Get the Lagrangian coordinate in the Lower Wall
        zeta[0] =  zeta_lo + double(j-n_x0-n_x1)*dzeta;
        //Get the geometric object and local coordinate
        Straight_wall_pt->locate_zeta(zeta,geometric_object_pt,s_wall);
        
        //The local coordinate is a geometric parameter
        //This needs to be set (rather than added) because the
        //same spine may be visited more than once
        Vector<double> parameters(1,s_wall[0]);
        nod_pt->spine_pt()->set_geom_parameter(parameters);
        
        // Get position of wall
        Straight_wall_pt->position(s_wall,r_wall);
        
        // Adjust spine height
        nod_pt->spine_pt()->height()=r_wall[1];
        
        // The sub geom object is one (and only) geom object
        // for spine:
        Vector<GeomObject*> geom_object_pt(1);
        geom_object_pt[0] = geometric_object_pt;
        
        // Pass geom object(s) to spine
        nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
       }
      }
    }
   //Loop over the nodes in the elements horizontally, ignoring 
   //the first column
   
   // Last spine needs special treatment in x-periodic meshes:
   unsigned n_pmax=n_p;
   if ((this->Xperiodic)&&(j==n_x-1)) n_pmax=n_p-1;
    
   for(unsigned l2=1;l2<n_pmax;l2++)
    {
     //Assign the new spine with unit height
     new_spine_pt=new Spine(1.0);
     new_spine_pt->spine_height_pt()->pin(0);
     Spine_pt.push_back(new_spine_pt);
     
     // Get the node
     SpineNode* nod_pt=element_node_pt(j,l2);
     //Set the pointer to spine
     nod_pt->spine_pt() = new_spine_pt;
     //Set the fraction
     nod_pt->fraction() = 0.0;
     // Pointer to the mesh that implements the update fct
     nod_pt->spine_mesh_pt() = this; 
     // Set update fct id
     nod_pt->node_update_fct_id()=2;

     {
      // Provide spine with additional storage for wall coordinate 
      // and wall geom object:
      
      // Get the Lagrangian coordinate in the Lower Wall
      zeta[0] =  zeta_lo + double(j-n_x0-n_x1)*dzeta + double(l2)*dzeta/(n_p-1.0);
      //Get the geometric object and local coordinate
      Straight_wall_pt->locate_zeta(zeta,geometric_object_pt,s_wall);
      
      //The local coordinate is a geometric parameter
      //This needs to be set (rather than added) because the
      //same spine may be visited more than once
      Vector<double> parameters(1,s_wall[0]);
      nod_pt->spine_pt()->set_geom_parameter(parameters);
      
      // Get position of wall
      Straight_wall_pt->position(s_wall,r_wall);
      
      // Adjust spine height
      nod_pt->spine_pt()->height()=r_wall[1];
      
      // The sub geom object is one (and only) geom object
      // for spine:
      Vector<GeomObject*> geom_object_pt(1);
      geom_object_pt[0] = geometric_object_pt;
      
      // Pass geom object(s) to spine
      nod_pt->spine_pt()->set_geom_object_pt(geom_object_pt);
     }
     
     //Loop vertically along the spine
     //Loop over the elements 
     for(unsigned long i=0;i<n_y;i++)
      {
       //Loop over the vertical nodes, apart from the first
       for(unsigned l1=1;l1<n_p;l1++)
        {
         // Get the node
         SpineNode* nod_pt=element_node_pt(i*n_x+j,l1*n_p+l2);
         //Set the pointer to the spine
         nod_pt->spine_pt() = new_spine_pt;
         //Set the fraction
         nod_pt->fraction()=(double(i)+double(l1)/double(n_p-1))/double(n_y);
         // Pointer to the mesh that implements the update fct
         nod_pt->spine_mesh_pt() = this; 
         // Set update fct id
         nod_pt->node_update_fct_id()=2;
        }  
      }
    }
  }

 // Increment number of previous elements
 n_prev_elements+=n_x2*n_y;

 // Last spine needs special treatment for periodic meshes
 // because it's the same as the first one...
 if (this->Xperiodic)
  {
   // Last spine is the same as first one...
   Spine* final_spine_pt=Spine_pt[0];  

   // Get the node
   SpineNode* nod_pt=element_node_pt((n_x-1),(n_p-1));

   //Set the pointer for the first node
   nod_pt->spine_pt() = final_spine_pt;
   //Set the fraction to be the same as for the nodes on the first row
   nod_pt->fraction() = element_node_pt(0,0)->fraction();
   // Pointer to the mesh that implements the update fct
   nod_pt->spine_mesh_pt() = element_node_pt(0,0)->spine_mesh_pt();
   
   //Now loop vertically along the spine
   for(unsigned i=0;i<n_y;i++)
    {
     //Loop over the vertical nodes, apart from the first
     for(unsigned l1=1;l1<n_p;l1++)
      {
       // Get the node
       SpineNode* nod_pt=element_node_pt(i*n_x+(n_x-1),l1*n_p+(n_p-1));

       //Set the pointer to the spine
       nod_pt->spine_pt() = final_spine_pt;
       //Set the fraction to be the same as in first row
       nod_pt->fraction() = element_node_pt(i*n_x,l1*n_p)->fraction();
       // Pointer to the mesh that implements the update fct
       nod_pt->spine_mesh_pt() = 
        element_node_pt(i*n_x,l1*n_p)->spine_mesh_pt();
      }
    }
   
  }



} // end of build_channel_spine_mesh


//======================================================================
/// Reorder the elements, so we loop over them vertically first
/// (advantageous in "wide" domains if a frontal solver is used).
//======================================================================
template <class ELEMENT>
void ChannelSpineMesh<ELEMENT>::element_reorder()
{

 unsigned n_x = this->Nx;
 unsigned n_y = this->Ny;
 //Find out how many elements there are
 unsigned long Nelement = nelement();
 //Find out how many fluid elements there are
 unsigned long Nfluid = n_x*n_y;
 //Create a dummy array of elements
 Vector<FiniteElement *> dummy;

 //Loop over the elements in horizontal order
 for(unsigned long j=0;j<n_x;j++)
  {
   //Loop over the elements in lower layer vertically
   for(unsigned long i=0;i<n_y;i++)
    {
     //Push back onto the new stack
     dummy.push_back(finite_element_pt(n_x*i + j));
    }

   //Push back the line element onto the stack
   dummy.push_back(finite_element_pt(Nfluid+j));
  }

 //Now copy the reordered elements into the element_pt
 for(unsigned long e=0;e<Nelement;e++)
  {
   Element_pt[e] = dummy[e];
  }

} // end of element_reorder


}

#endif