channel_with_leaflet_domain.h

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_WITH_LEAFLET_DOMAIN_HEADER
#define OOMPH_CHANNEL_WITH_LEAFLET_DOMAIN_HEADER


// Generic includes
#include "../generic/geom_objects.h"
#include "../generic/domain.h"
#include "../generic/macro_element.h"

namespace oomph
{

//===========================================================
/// Rectangular domain with a leaflet blocking the lower
/// half
//===========================================================
class ChannelWithLeafletDomain : public Domain
{

public:
 
 ///\short Constructor: Pass pointer to GeomObject that represents the leaflet,
 /// the length of the domain to left and right of the leaflet, the
 /// height of the leaflet and the overall height of the channel,
 /// the number of element columns to the left and right of the leaflet,
 /// the number of rows of elements from the bottom of the channel to 
 /// the end of the leaflet, the number of rows of elements above the
 /// end of the leaflet.
 ChannelWithLeafletDomain(GeomObject* leaflet_pt, const double& lleft,
                          const double& lright, const double& hleaflet,
                          const double& htot,
                          const unsigned& nleft, const unsigned& nright,
                          const unsigned& ny1, const unsigned&  ny2)
  {
   // Copy assignments into private storage
   Lleft = lleft; 
   Lright=lright;
   Hleaflet=hleaflet; 
   Htot=htot; 
   Nleft = nleft; 
   Nright = nright; 
   Ny1 = ny1; 
   Ny2 = ny2;
   Leaflet_pt =  leaflet_pt;
   
   // Store origin of leaflet for fast reference
   Vector<double> zeta(1);
   zeta[0]=0.0;
   Vector<double> r(2);
   Leaflet_pt->position(zeta,r);
   X_0 = r[0];

   // Number of macro elements
   unsigned nmacro = (Ny1+Ny2)*(Nleft+Nright);
   Macro_element_pt.resize(nmacro);

   // Build the 2D macro elements
   for (unsigned i=0;i<nmacro;i++)
    {
     Macro_element_pt[i]= new QMacroElement<2>(this,i);
    }
 
  }


 /// Destructor: Kill macro elements
 ~ChannelWithLeafletDomain()
  {
   unsigned nmacro = (Ny1+Ny2)*(Nleft+Nright);
   for (unsigned i=0;i<nmacro;i++){delete Macro_element_pt[i];}
  }


 /// Total height of domain (width of channel)
 double htot(){return Htot;}
 
 /// Height of leaflet
 double hleaflet(){return Hleaflet;}

 /// Length of domain to the left of leaflet
 double lleft(){return Lleft;}

 /// Length of domain to the right of leaflet
 double lright(){return Lright;}
 
 /// Pointer to the wall
 GeomObject*& leaflet_pt() { return Leaflet_pt;};

 /// Parametrisation of macro element boundaries
 void macro_element_boundary(const unsigned& t,
                             const unsigned& imacro,
                             const unsigned& idirect,
                             const Vector<double>& zeta,
                             Vector<double>& r);

protected : 
 
 /// Helper function 
 void macro_bound_I_N(const unsigned& t, const Vector<double>& zeta,
                      Vector<double>& r, const  unsigned& i,const unsigned& j); 

 /// Helper function 
 void macro_bound_I_S(const unsigned& t,const Vector<double>& zeta,
                      Vector<double>& r,const unsigned& i,const  unsigned& j);

 /// Helper function 
 void macro_bound_I_W(const unsigned& t,const Vector<double>& zeta,
                      Vector<double>& r,const  unsigned& i,const  unsigned& j);

 /// Helper function  
 void macro_bound_I_E(const unsigned& t,const Vector<double>& zeta,
                      Vector<double>& r,const unsigned& i,
                      const  unsigned& j);

 /// Helper function 
 void macro_bound_II_N(const unsigned& t, const Vector<double>& zeta,
                       Vector<double>& r,
                       const  unsigned& i,const unsigned& j); 

 /// Helper function 
 void macro_bound_II_S(const unsigned& t,const Vector<double>& zeta,
                       Vector<double>& r,const unsigned& i,const  unsigned& j);

 /// Helper function 
 void macro_bound_II_W(const unsigned& t,const Vector<double>& zeta,
                       Vector<double>& r,const unsigned& i,const  unsigned& j);

 /// Helper function 
 void macro_bound_II_E(const unsigned& t,const Vector<double>& zeta,
                       Vector<double>& r,const unsigned& i,
                       const  unsigned& j);
 
 /// Helper function 
 void macro_bound_III_N(const unsigned& t, const Vector<double>& zeta,
                        Vector<double>& r,
                        const  unsigned& i,const unsigned& j); 

 /// Helper function 
 void macro_bound_III_S(const unsigned& t,const Vector<double>& zeta,
                        Vector<double>& r,const unsigned& i,const unsigned& j);

 /// Helper function 
 void macro_bound_III_W(const unsigned& t,const Vector<double>& zeta,
                        Vector<double>& r,const unsigned& i,const unsigned& j);

 /// Helper function 
 void macro_bound_III_E(const unsigned& t,const Vector<double>& zeta,
                        Vector<double>& r,const unsigned& i,
                        const unsigned& j);

 /// Helper function 
 void macro_bound_IV_N(const unsigned& t, const Vector<double>& zeta,
                       Vector<double>& r,
                       const  unsigned& i,const unsigned& j); 

 /// Helper function 
 void macro_bound_IV_S(const unsigned& t,const Vector<double>& zeta,
                       Vector<double>& r,const unsigned& i,const unsigned& j);

 /// Helper function 
 void macro_bound_IV_W(const unsigned& t,const Vector<double>& zeta,
                       Vector<double>& r,const unsigned& i,const unsigned& j);

 /// Helper function 
 void macro_bound_IV_E(const unsigned& t,const Vector<double>& zeta,
                       Vector<double>& r,const unsigned& i,
                       const unsigned& j);
 /// Helper function 
 void slanted_bound_up(const unsigned& t,
                       const Vector<double>& zeta,
                       Vector<double>& r);
                            

 
 /// Length of the domain to the right of the leaflet
 double Lright;

 /// Length of the domain to the left of the leaflet
 double Lleft;

 /// Lagrangian coordinate at end of leaflet
 double Hleaflet;

 /// Total width of the channel
 double Htot; 
 
 ///Number of macro element columnns to the right of the leaflet 
 unsigned Nright;

 ///Number of macro element columns to the left of the leaflet
 unsigned Nleft;

 ///Number of macro element rows up to the end of the leaflet
 unsigned Ny1;

 ///Number of macro element rows above the leaflet
 unsigned Ny2;

 /// \short Center of the domain : origin of the leaflet, extracted
 /// from GeomObject and stored for fast access.
 double X_0;

 /// Pointer to leaflet
 GeomObject*  Leaflet_pt;

};



//===================================================================
/// Parametrisation of macro element boundaries
//===================================================================
void ChannelWithLeafletDomain::macro_element_boundary(
 const unsigned& t,
 const unsigned& imacro,
 const unsigned& idirect,
 const Vector<double>& zeta,
 Vector<double>& r)
{

#ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
   // Warn about time argument being moved to the front
   OomphLibWarning(
    "Order of function arguments has changed between versions 0.8 and 0.85",
    "ChannelWithLeafletDomain::macro_element_boundary(...)",
    OOMPH_EXCEPTION_LOCATION);
#endif

 using namespace QuadTreeNames;

 //Number of the line and of the colum in the whole domain
 //Beware : iline starts on zero
 unsigned iline = int(floor(double(imacro)/double(Nleft+Nright)));
 unsigned icol = imacro%(Nright+Nleft);
 
 //Number of the line and of the colum in the part considered
 unsigned i,j;
 
 //Left low part of the domain : Part I
 //------------------------------------
  if( (iline<Ny1) && (icol<Nleft) )
   {
    i = iline;
    j=icol;
    
    switch(idirect)
     {
     case N:
      macro_bound_I_N(t,zeta,r,i,j);
      break;
     case S:
      macro_bound_I_S(t,zeta,r,i,j);
      break;
     case W:
      macro_bound_I_W(t,zeta,r,i,j);
      break;
     case E :
      macro_bound_I_E(t,zeta,r,i,j);  
      break;
     }
   }
  //Right low part of the domain : Part II
  //--------------------------------------
  else if( (iline<Ny1) && (icol>=Nleft) )
   {
    i = iline;
    j=icol-Nleft;
    
    switch(idirect)
     {
     case N:
      macro_bound_II_N(t,zeta,r,i,j);
      break;
     case S:
      macro_bound_II_S(t,zeta,r,i,j);
      break;
     case W:
      macro_bound_II_W(t,zeta,r,i,j);
      break;
     case E :
      macro_bound_II_E(t,zeta,r,i,j);  
      break;
     }
   }
 //Left upper part of the domain : Part III
 //----------------------------------------
  else if( (iline>=Ny1) && (icol<Nleft) )
   {
    i = iline-Ny1;
    j=icol;
    
    switch(idirect)
     {
     case N:
      macro_bound_III_N(t,zeta,r,i,j);
      break;
     case S:
      macro_bound_III_S(t,zeta,r,i,j);
      break;
     case W:
      macro_bound_III_W(t,zeta,r,i,j);
      break;
     case E :
      macro_bound_III_E(t,zeta,r,i,j);  
      break;
     }
   }
  //Right upper part of the domain : Part IV
  //-----------------------------------------
  else if( (iline>=Ny1) && (icol>=Nleft) )
   {
    i = iline-Ny1;
    j=icol-Nleft;
    
    switch(idirect)
     {
     case N:
      macro_bound_IV_N(t,zeta,r,i,j);
      break;
     case S:
      macro_bound_IV_S(t,zeta,r,i,j);
      break;
     case W:
      macro_bound_IV_W(t,zeta,r,i,j);
      break;
     case E :
      macro_bound_IV_E(t,zeta,r,i,j);  
      break;
     }
   }
}  



/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
// Helper functions for region I (lower left region)
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////



//=====================================================================
/// Helper function for eastern boundary in lower left region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_I_E(const unsigned& t, 
                                               const Vector<double>& zeta,
                                               Vector<double>& r, 
                                               const unsigned& i,
                                               const unsigned& j)
{
 
 //Find x,y on the wall corresponding to the position of the macro element
 
 //xi_wall varies from xi0 to xi1 on the wall
 double xi0, xi1;
 xi0 = double(i)*Hleaflet/double(Ny1);
 xi1 = double(i+1)* Hleaflet/double(Ny1);
 
 Vector<double> xi_wall(1);
 xi_wall[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);
 
 Vector<double> r_wall(2);
 Leaflet_pt->position(t,xi_wall,r_wall);
 
 //Find x,y on a vertical line corresponding 
 //to the position of the macro element

 //the vertical line goes from y0 to y1
 double y0,y1;
 y0 = double(i)*Hleaflet/double(Ny1);
 y1 = double(i+1)*Hleaflet/double(Ny1);
 
 Vector<double> r_vert(2);
 r_vert[0] = -Lleft+X_0;
 r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);
 
//Parameter with value 0 in -Lleft and value 1 on the wall.
 double s = double(j+1)/double(Nleft);
 
///Final expression of r
 r[0] = r_vert[0] + s*( r_wall[0]-r_vert[0]);
 r[1] = r_vert[1] + s*( r_wall[1]-r_vert[1]);
}


//=====================================================================
/// Helper function for western boundary in lower left region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_I_W(const unsigned& t, 
                                               const Vector<double>& zeta,
                                               Vector<double>& r, 
                                               const unsigned& i,
                                               const unsigned& j)
{
 //Find x,y on the wall corresponding to the position of the macro element
 
 //xi_wall varies from xi0 to xi1 on the wall
 double xi0, xi1;
 xi0 = double(i)*Hleaflet/double(Ny1);
 xi1 = double(i+1)*Hleaflet/double(Ny1);

 Vector<double> xi_wall(1);
 xi_wall[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);

 Vector<double> r_wall(2);
 Leaflet_pt->position(t,xi_wall,r_wall);
 
 // Find x,y on a vertical line corresponding 
 //to the position of the macro element
 
 //the vertical line goes from y0 to y1
 double y0,y1;
 y0 = double(i)*Hleaflet/double(Ny1);
 y1 = double(i+1)*Hleaflet/double(Ny1);
 
 Vector<double> r_vert(2);
 r_vert[0] = -Lleft+X_0;
 r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);
 
 //Parameter with value 0 in -Lleft and value 1 on the wall.
 double s = double(j)/double(Nleft);
 
 ///Final expression of r
 r[0] = r_vert[0] + s*( r_wall[0]-r_vert[0]);
 r[1] = r_vert[1] + s*( r_wall[1]-r_vert[1]);

}



//=====================================================================
/// Helper function for northern boundary in lower left region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_I_N(const unsigned& t, 
                                               const Vector<double>& zeta,
                                               Vector<double>& r, 
                                               const unsigned& i,
                                               const unsigned& j)
{
 //Find the coordinates of the two corners of the north boundary 
 Vector<double> xi(1);
 Vector<double> r_left(2);
 Vector<double> r_right(2);
 xi[0]=1;
 macro_bound_I_W(t,xi,r_left,i,j);
 macro_bound_I_E(t,xi,r_right,i,j);
 
//Connect those two points with a straight line
 r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
 r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);
 
}


//=====================================================================
/// Helper function for southern boundary in lower left region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_I_S(const unsigned& t, 
                                               const Vector<double>& zeta,
                                               Vector<double>& r, 
                                               const unsigned& i,
                                               const unsigned& j)
{
///Find the coordinates of the two corners of the south boundary 
 Vector<double> xi(1);
 Vector<double> r_left(2);
 Vector<double> r_right(2);
 xi[0]=-1.0;
 macro_bound_I_W(t,xi,r_left,i,j);
 macro_bound_I_E(t,xi,r_right,i,j);
 
//Connect those two points with a straight line
 r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
 r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);
}






/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
// Helper functions for region II (lower right region)
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////



//=====================================================================
/// Helper function for eastern boundary in lower right region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_II_E(const unsigned& t,
                                                const Vector<double>& zeta,
                                                Vector<double>& r, 
                                                const unsigned& i,
                                                const unsigned& j)
 
{
 //Find x,y on the wall corresponding to the position of the macro element
 
 //xi_wall varies from xi0 to xi1 on the wall
 double xi0, xi1;
 xi0 = double(i)*Hleaflet/double(Ny1);
 xi1 = double(i+1)*Hleaflet/double(Ny1);
 
 Vector<double> xi_wall(1);
 xi_wall[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);
 
 Vector<double> r_wall(2);
 Leaflet_pt->position(t,xi_wall,r_wall);
 
 // Find x,y on a vertical line corresponding 
 // to the position of the macro element
 
 //the vertical line goes from y0 to y1
 double y0,y1;
 y0 = double(i)*Hleaflet/double(Ny1);
 y1 = double(i+1)*Hleaflet/double(Ny1);
 
 Vector<double> r_vert(2);
 r_vert[0] = Lright+X_0;
 r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);
 
//Parameter with value 0 in Lright and value 1 on the wall.
 double s = double(Nright-j-1)/double(Nright);  /***Change****/
 
///Final expression of r
 r[0] = r_vert[0] + s*( r_wall[0]-r_vert[0]);
 r[1] = r_vert[1] + s*( r_wall[1]-r_vert[1]);
}




//=====================================================================
/// Helper function for western boundary in lower right region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_II_W(const unsigned& t, 
                                                const Vector<double>& zeta,
                                                Vector<double>& r, 
                                                const unsigned& i,
                                                const unsigned& j)
{
 //Abscissa of the origin of the boudary east
 
 //Find x,y on the wall corresponding to the position of the macro element
 
 //xi_wall varies from xi0 to xi1 on the wall
 double xi0, xi1;
 xi0 = double(i)*Hleaflet/double(Ny1);
 xi1 = double(i+1)*Hleaflet/double(Ny1);

 Vector<double> xi_wall(1);
 xi_wall[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);

 Vector<double> r_wall(2);
 Leaflet_pt->position(t,xi_wall,r_wall);
 
 // Find x,y on a vertical line corresponding 
 // to the position of the macro element
 
 //the vertical line goes from y0 to y1
 double y0,y1;
 y0 = double(i)*Hleaflet/double(Ny1);
 y1 = double(i+1)*Hleaflet/double(Ny1);
 
 Vector<double> r_vert(2);
 r_vert[0] = Lright+X_0; 
 r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);

 //Parameter with value 0 in -Lleft and value 1 on the wall.
 double s = double(Nright-j)/double(Nright);     /***Change****/
 
 //Final expression of r
 r[0] = r_vert[0] + s*( r_wall[0]-r_vert[0]);
 r[1] = r_vert[1] + s*( r_wall[1]-r_vert[1]);
 
}



//=====================================================================
/// Helper function for northern boundary in lower right region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_II_N(const unsigned& t, 
                                                const Vector<double>& zeta,
                                                Vector<double>& r,
                                                const unsigned& i,
                                                const unsigned& j)
{
 //Find the coordinates of the two corners of the north boundary 
 Vector<double> xi(1);
 Vector<double> r_left(2);
 Vector<double> r_right(2);
 xi[0]=1;
 macro_bound_II_W(t,xi,r_left,i,j);
 macro_bound_II_E(t,xi,r_right,i,j);
 
//Connect those two points with a straight line
 r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
 r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);

}



//=====================================================================
/// Helper function for southern boundary in lower right region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_II_S(const unsigned& t, 
                                                const Vector<double>& zeta,
                                                Vector<double>& r,
                                                const unsigned& i,
                                                const unsigned& j)
{
 //Find the coordinates of the two corners of the south boundary 
 Vector<double> xi(1);
 Vector<double> r_left(2);
 Vector<double> r_right(2);
 xi[0]=-1.0;
 macro_bound_II_W(t,xi,r_left,i,j);
 macro_bound_II_E(t,xi,r_right,i,j);

//Connect those two points with a straight line
 r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
 r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);
}







/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
// Helper functions for region III (upper left region)
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////



//=====================================================================
/// Describe the line between the boundary north of the domain (at x=X_0)
/// and the top of the wall, when zeta goes from 0 to 1.
//=====================================================================
void ChannelWithLeafletDomain::slanted_bound_up(const unsigned& t,
                                                const Vector<double>& zeta,
                                                Vector<double>& r)
{
 // Coordinates of the point on the boundary beetween the upper
 // and the lower part, in the same column, at the east.
 Vector<double> xi(1);
 xi[0]=Hleaflet;

 Vector<double> r_join(2);

 Leaflet_pt->position(t,xi,r_join);

 r[0] = r_join[0] + zeta[0]*(X_0-r_join[0]);
 r[1] = r_join[1] + zeta[0]*(Htot-r_join[1]);
}



//=====================================================================
/// Helper function for eastern boundary in upper left region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_III_E(const unsigned& t,
                                                 const Vector<double>& zeta,
                                                 Vector<double>& r, 
                                                 const unsigned& i,
                                                 const unsigned& j)
{
 // Find x,y on the slanted straight line (SSL) corresponding to
 // the position of the macro element

 //xi_line varies from xi0 to xi1 on the SSL
 double xi0, xi1;
 xi0 = double(i)/double(Ny2);
 xi1 = double(i+1)/double(Ny2);

 Vector<double> xi_line(1);
 xi_line[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);

 Vector<double> r_line(2);
 slanted_bound_up(t,xi_line,r_line);

 // Find x,y on a vertical line corresponding 
 // to the position of the macro element

 //the vertical line goes from y0 to y1
 double y0,y1;
 y0 = double(i)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;
 y1 = double(i+1)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;;

 Vector<double> r_vert(2);
 r_vert[0] = -Lleft+X_0;
 r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);

//Parameter with value 0 in Lright and value 1 on the wall.
 double s = double(j+1)/double(Nleft);  /***Change****/

///Final expression of r
 r[0] = r_vert[0] + s*( r_line[0]-r_vert[0]);
 r[1] = r_vert[1] + s*( r_line[1]-r_vert[1]);
}




//=====================================================================
/// Helper function for western boundary in upper left region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_III_W(const unsigned& t,
                                                 const Vector<double>& zeta,
                                                 Vector<double>& r, 
                                                 const unsigned& i,
                                                 const unsigned& j)
{
 // Find x,y on the slanted straight line (SSL) corresponding to
 // the position of the macro element
 
 //xi_line varies from xi0 to xi1 on the SSL
 double xi0, xi1;
 xi0 = double(i)/double(Ny2);
 xi1 = double(i+1)/double(Ny2);
 
 Vector<double> xi_line(1);
 xi_line[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);
 
 Vector<double> r_line(2);
 slanted_bound_up(t,xi_line,r_line);
 
 // Find x,y on a vertical line corresponding 
 // to the position of the macro element
 
 // the vertical line goes from y0 to y1
 double y0,y1;
 y0 = double(i)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;
 y1 = double(i+1)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;;

 Vector<double> r_vert(2);
 r_vert[0] = -Lleft+X_0;
 r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);
 
//Parameter with value 0 in Lright and value 1 on the wall.
 double s = double(j)/double(Nleft);  /***Change****/
 
 //Final expression of r
 r[0] = r_vert[0] + s*( r_line[0]-r_vert[0]);
 r[1] = r_vert[1] + s*( r_line[1]-r_vert[1]);
}




//=====================================================================
/// Helper function for northern boundary in upper left region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_III_N(const unsigned& t, 
                                                 const Vector<double>& zeta,
                                                 Vector<double>& r, 
                                                 const unsigned& i,
                                                 const unsigned& j)
{
 //Find the coordinates of the two corners of the north boundary 
 Vector<double> xi(1);
 Vector<double> r_left(2);
 Vector<double> r_right(2);
 xi[0]=1;
 macro_bound_III_W(t,xi,r_left,i,j);
 macro_bound_III_E(t,xi,r_right,i,j);

 //Connect those two points with a straight line
 r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
 r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);

}


//=====================================================================
/// Helper function for southern boundary in upper left region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_III_S(const unsigned& t,
                                                 const Vector<double>& zeta,
                                                 Vector<double>& r, 
                                                 const unsigned& i,
                                                 const unsigned& j)
{
 //Find the coordinates of the two corners of the south boundary 
 Vector<double> xi(1);
 Vector<double> r_left(2);
 Vector<double> r_right(2);
 xi[0]=-1;
 macro_bound_III_W(t,xi,r_left,i,j);
 macro_bound_III_E(t,xi,r_right,i,j);

 //Connect those two points with a straight line
 r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
 r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);
}







/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
// Helper functions for region IV (upper right region)
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////


//=====================================================================
/// Helper function for eastern boundary in upper right region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_IV_E(const unsigned& t,
                                                const Vector<double>& zeta,
                                                Vector<double>& r, 
                                                const unsigned& i,
                                                const unsigned& j)
{
 // Find x,y on the slanted straight line (SSL) corresponding to
 // the position of the macro element

 //xi_line varies from xi0 to xi1 on the SSL
 double xi0, xi1;
 xi0 = double(i)/double(Ny2);
 xi1 = double(i+1)/double(Ny2);

 Vector<double> xi_line(1);
 xi_line[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);

 Vector<double> r_line(2);
 slanted_bound_up(t,xi_line,r_line);

 // Find x,y on a vertical line corresponding 
 // to the position of the macro element

 //the vertical line goes from y0 to y1
 double y0,y1;
 y0 = double(i)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;
 y1 = double(i+1)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;;

 Vector<double> r_vert(2);
 r_vert[0] = Lright+X_0;
 r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);

 //Parameter with value 0 in Lright and value 1 on the wall.
 double s = double(Nright -j-1)/double(Nright);  /***Change****/

 //Final expression of r
 r[0] = r_vert[0] + s*( r_line[0]-r_vert[0]);
 r[1] = r_vert[1] + s*( r_line[1]-r_vert[1]);
}



//=====================================================================
/// Helper function for western boundary in upper right region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_IV_W(const unsigned& t,
                                                const Vector<double>& zeta,
                                                Vector<double>& r,
                                                const unsigned& i,
                                                const unsigned& j)
{
 // Find x,y on the slanted straight line (SSL) corresponding to
 // the position of the macro element

 //xi_line varies from xi0 to xi1 on the SSL
 double xi0, xi1;
 xi0 = double(i)/double(Ny2);
 xi1 = double(i+1)/double(Ny2);

 Vector<double> xi_line(1);
 xi_line[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);

 Vector<double> r_line(2);
 slanted_bound_up(t,xi_line,r_line);

 // Find x,y on a vertical line corresponding 
 // to the position of the macro element

 // The vertical line goes from y0 to y1
 double y0,y1;
 y0 = double(i)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;
 y1 = double(i+1)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;;
 
 Vector<double> r_vert(2);
 r_vert[0] = Lright+X_0;
 r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);

 //Parameter with value 0 in Lright and value 1 on the wall.
 double s = double(Nright-j)/double(Nright);  /***Change****/

 //Final expression of r
 r[0] = r_vert[0] + s*( r_line[0]-r_vert[0]);
 r[1] = r_vert[1] + s*( r_line[1]-r_vert[1]);
}



//=====================================================================
/// Helper function for northern boundary in upper right region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_IV_N(const unsigned& t, 
                                                const Vector<double>& zeta,
                                                Vector<double>& r, 
                                                const unsigned& i,
                                                const unsigned& j)
{
 // Find the coordinates of the two corners of the north boundary 
 Vector<double> xi(1);
 Vector<double> r_left(2);
 Vector<double> r_right(2);
 xi[0]=1;
 macro_bound_IV_W(t,xi,r_left,i,j);
 macro_bound_IV_E(t,xi,r_right,i,j);

 //Connect those two points with a straight line
 r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
 r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);

}


//=====================================================================
/// Helper function for southern boundary in upper right region
//=====================================================================
void ChannelWithLeafletDomain::macro_bound_IV_S(const unsigned& t, 
                                                const Vector<double>& zeta,
                                                Vector<double>& r, 
                                                const unsigned& i,
                                                const unsigned& j)
{
 //Find the coordinates of the two corners of the south boundary 
 Vector<double> xi(1);
 Vector<double> r_left(2);
 Vector<double> r_right(2);
 xi[0]=-1;
 macro_bound_IV_W(t,xi,r_left,i,j);
 macro_bound_IV_E(t,xi,r_right,i,j);

 //Connect those two points with a straight line
 r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
 r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);
}


}


#endif