solvers/PulseWaveSolver/EquationSystems/ArterialPressureArea.cpp

 ///////////////////////////////////////////////////////////////////////////////
 //
-// File CommMpi.cpp
+// File ArterialPressureArea.cpp
 //
 // For more information, please see: http://www.nektar.info
 //
@@ -34,7 +34,7 @@
 ///////////////////////////////////////////////////////////////////////////////
 
 #include <PulseWaveSolver/EquationSystems/ArterialPressureArea.h>
-
+#define PI 3.14159265359
 namespace Nektar
 {
 
@@ -48,9 +48,14 @@ namespace Nektar
      *
      */
     ArterialPressureArea::ArterialPressureArea(Array<OneD, MultiRegions::ExpListSharedPtr> pVessel, 
-                       const LibUtilities::SessionReaderSharedPtr pSession)
-        : PulseWavePressureArea(pVessel,pSession)
+                                               const LibUtilities::SessionReaderSharedPtr pSession,
+                                               const int nDomains)
+        : PulseWavePressureArea(pVessel,pSession,nDomains)
     {
+        m_session->LoadParameter("rho", m_rho, 0.5);
+
+        m_beta = Array<OneD, Array<OneD, NekDouble> >(nDomains);
+        m_beta_trace = Array<OneD, Array<OneD, NekDouble> >(nDomains);
     }
 
     /**
@@ -61,8 +66,337 @@ namespace Nektar
 
     }
 
-    void ArterialPressureArea::v_DoPressure()
-    { 
+    void ArterialPressureArea::v_ReadParameters(int omega, int nqTrace, MultiRegions::ExpListSharedPtr &field)
+    {
+        int nq=m_vessels[2*omega]->GetNpoints();
+        m_beta[omega] = Array<OneD, NekDouble>(nq);
+        EvaluateFunction(m_vessels,m_session,"beta", m_beta[omega],"MaterialProperties",0.0,omega,nq);
+
+        m_beta_trace[omega] = Array<OneD, NekDouble>(nqTrace);
+        field->ExtractTracePhys(m_beta[omega],m_beta_trace[omega]);
+    }
+
+    void ArterialPressureArea::v_GetPacons(int omega, Array<OneD, Array<OneD, NekDouble> > &pacons_trace, int nqTrace)
+    {
+        int nparams=1;
+        pacons_trace = Array<OneD, Array<OneD, NekDouble> >(nparams);
+
+        pacons_trace[0] = Array<OneD, NekDouble>(nqTrace);
+        Vmath::Vcopy(nqTrace,m_beta_trace[omega],1,pacons_trace[0],1);
+    }
+
+
+
+
+    /**
+     *  Contains the definition of the pressure-area relationship. To implement a new pressure-area relationship
+     *  change the mathematical definitions here and also in getdpda. Also consider changing: getAfromPressure
+     *  and getCharIntegral; currently use numerical methods but it is preferable to enter analyticall expressions
+     *  if they can be derived for the chosen p-A relation. Furthermore ensure that the definition of the lower
+     *  integration boundary A_0 is updated in getAfromPressure, getAfromChars and getCharIntegral if using the
+     *  numerical methods to ensure it refers to the correct element of pacons; typically can be assumed to be
+     *  equal to the passive area.
+     *
+     *  The p-A relation cannot be a function of u and can contain any number of parameters as specified by
+     *  P-A_num_Params in the input file which should then be specified for each subdomain using:
+     *  <FUNCTION NAME="P-A_Vals">
+     *  <E VAR="Val[Domain ID][Parameter ID]" VALUE="" />
+     *  </FUNCTION>
+     */
+    
+    void ArterialPressureArea::v_getPressure(NekDouble A, Array<OneD, NekDouble> pacons, NekDouble &pressure)
+    {
+        
+        /*
+         * pacons is an array which holds, for the point for which at which we are evaluating this function,
+         * all of the constants found in the pressure-area relation. They are in the same order as the input file.
+         */
+        
+        //pacons[0] = p_ext
+        //pacons[1] = E
+        //pacons[2] = A_0
+        //pacons[3] = h0
+        //pacons[4] = nue
+        
+        NekDouble Passive = 0.0;
+        NekDouble Beta_Const = 0.0;
+        
+        Beta_Const = ((sqrt(PI)*pacons[3]*pacons[1])/((1-pow(pacons[4],2))*pacons[2]));
+        
+        Passive = (Beta_Const*(sqrt(A)-sqrt(pacons[2])));
+        
+        pressure = pacons[0] + Passive;
+        
+        //cout<<"Working?"<<pressure<<endl;
+        
+    }
+    
+    /**
+     *  Contains the definition of dp/da - needs to be updated when chaning p-A relation!
+     */
+    void ArterialPressureArea::v_getdpda(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dpda)
+    {
+        NekDouble Passive = 0.0;
+        NekDouble Beta_Const = 0.0;
+        
+        Beta_Const = ((sqrt(PI)*pacons[3]*pacons[1])/((1-pow(pacons[4],2))*pacons[2]));
+        Passive = (Beta_Const*0.5*(1/sqrt(A)));
+        
+        dpda = Passive;
+        
+    }
+    
+    /**
+     *  Calculates A in terms of pressure for for any p-A relation using a Newton iteration. If using a different
+     *  p-A relation the definition of A_0 may need modifying.
+     */
+    void ArterialPressureArea::v_getAfromPressure(NekDouble pressure, Array<OneD, NekDouble> pacons, NekDouble &A)
+    {
+        
+        //A = (pressure - pacons[2])/pacons[0] + sqrt(pacons[1]);
+        //A = A*A;
+        
+        NekDouble fa = 0.0;
+        NekDouble dpda = 0.0;
+        NekDouble p_calc = 0.0;
+        NekDouble delta_A_calc = 0.0;
+        NekDouble A_0 = 0.0;
+        
+        // Initialise to A_0
+        A_0=PI*pow(pacons[2],2)/40;
+        A=A_0;
+        
+        int proceed = 1;
+        int iter = 0;
+        int MAX_ITER = 200;
+        
+        // Tolerances for the algorithm
+        NekDouble Tol = 1.0e-10;
+        
+        while ((proceed) && (iter < MAX_ITER))
+        {
+            iter =iter+1;
+            
+            getPressure(A, pacons, p_calc);
+            fa = p_calc-pressure;
+            getdpda(A, 0, pacons, dpda);
+            
+            delta_A_calc=fa/dpda;
+            A = A - delta_A_calc;
+            
+            if (sqrt(delta_A_calc*delta_A_calc) < Tol)
+                proceed = 0;
+        }
+        
+    }
+    
+    /**
+     *  Calculates the two characteristic variables for any p-A relation, currently using numerical
+     *  integration. Unlike lambda and the derivitives we calculate both in the same function to avoid having
+     *  to do the numerical integration twice.
+     */
+    void ArterialPressureArea::v_getW(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &W1, NekDouble &W2)
+    {
+        NekDouble integralTerm=0.0;
+        getCharIntegral(A, u, pacons, integralTerm);
+        W1 = u + integralTerm;
+        W2 = u - integralTerm;
+    }
+    
+    /**
+     *  Contains the definition of A in terms of W_1 and W_2 for any p-A relation though definition of A_0 may
+     *  need updating. Solved for using a Newton iteration and numerical integration and numerical integration.
+     */
+    void ArterialPressureArea::v_getAfromChars(NekDouble W_1, NekDouble W_2, Array<OneD, NekDouble> pacons, NekDouble &A)
+    {
+        
+        //    A = (W_1 - W_2)*sqrt(m_rho/(32*pacons[0]))+sqrt(sqrt(pacons[1]));
+        //    A = A*A;
+        //    A = A*A;
+        
+        NekDouble fa = 0.0;
+        NekDouble dfa = 0.0;
+        NekDouble dpda = 0.0;
+        NekDouble integralTerm = 0.0;
+        NekDouble delta_A_calc = 0.0;
+        NekDouble A_0 = 0.0;
+        
+        // Initialise to A_0
+        A_0=PI*pow(pacons[2],2)/40;
+        A=A_0;
+        
+        int proceed = 1;
+        int iter = 0;
+        int MAX_ITER = 200;
+        
+        // Tolerances for the algorithm
+        NekDouble Tol = 1.0e-10;
+        
+        while ((proceed) && (iter < MAX_ITER))
+        {
+            iter =iter+1;
+            
+            getCharIntegral(A, 0, pacons, integralTerm);
+            fa = integralTerm - 0.5*(W_1-W_2);
+            getdpda(A, 0, pacons, dpda);
+            dfa=sqrt(dpda/(m_rho*A));
+            
+            delta_A_calc=fa/dfa;
+            A = A - delta_A_calc;
+            
+            if (sqrt(delta_A_calc*delta_A_calc) < Tol)
+                proceed = 0;
+        }
+        
+    }
+    
+    /**
+     *  Contains the definition of u in terms of W_1 and W_2 for any p-A relation
+     */
+    void ArterialPressureArea::v_getufromChars(NekDouble W_1, NekDouble W_2, Array<OneD, NekDouble> pacons, NekDouble &u)
+    {
+        
+        u = 0.5*(W_1+W_2);
+        
+    }
+    
+    /**
+     *  Contains the definition of the the first eigenvalue lambda_1 for any p-A relation
+     */
+    void ArterialPressureArea::v_getlambda1(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &lambda_1)
+    {
+        NekDouble dpda=0.0;
+        getdpda(A, u, pacons, dpda);
+        
+        lambda_1 = u + sqrt((A/m_rho)*dpda);
+        
+    }
+    
+    /**
+     *  Contains the definition of the the second eigenvalue lambda_2 for any p-A relation
+     */
+    void ArterialPressureArea::v_getlambda2(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &lambda_2)
+    {
+        NekDouble dpda=0.0;
+        getdpda(A, u, pacons, dpda);
+        
+        lambda_2 = u - sqrt((A/m_rho)*dpda);
+    }
+    
+    
+    /**
+     *  Contains the definition of dp/du for any p-A relation
+     */
+    void ArterialPressureArea::v_getdpdu(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dpdu)
+    {
+        
+        dpdu = 0;
+        
+    }
+    
+    /**
+     *  Contains the definition of dW1/da for any p-A relation
+     */
+    void ArterialPressureArea::v_getdW1da(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW1da)
+    {
+        NekDouble dpda=0.0;
+        getdpda(A, u, pacons, dpda);
+        
+        dW1da = sqrt(dpda/(A*m_rho));
+        
+    }
+    
+    /**
+     *  Contains the definition of dW1/du for any p-A relation
+     */
+    void ArterialPressureArea::v_getdW1du(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW1du)
+    {
+        
+        dW1du = 1;
+        
+    }
+    
+    /**
+     *  Contains the definition of dW2/da for any p-A relation
+     */
+    void ArterialPressureArea::v_getdW2da(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW2da)
+    {
+        
+        NekDouble dpda=0.0;
+        getdpda(A, u, pacons, dpda);
+        
+        dW2da = -1.0*sqrt(dpda/(A*m_rho));
+        
+    }
+    
+    /**
+     *  Contains the definition of dW2/du for any p-A relation
+     */
+    void ArterialPressureArea::v_getdW2du(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW2du)
+    {
+        
+        dW2du = 1;
+        
+    }
+    
+    /**
+     *  Evalautes the integral \int^A_{A_{0}}{ \sqrt{\frac{1}{\rho A} \frac{\partial p\left(A,t\right)}{\partial A}} dA}
+     *  which is found in the general definition of the characteristic variables.
+     *  Currenly does this numerically but if p-A relation allows suggest the analytical expression is entered instead.
+     *  This current implementation is suitable for any p-A relation, though A_0 may need updating as well as num of
+     *  trapeziums.
+     */
+    void ArterialPressureArea::v_getCharIntegral(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &intergralTerm)
+    {
+        
+        int numTrapeziums = 50;
+        NekDouble A_0 = 0.1*PI*pow(pacons[2],2)/4; //For active contractions can be less than static area
+        NekDouble A_n = A_0;
+        NekDouble dA = (A - A_n)/numTrapeziums;
+        NekDouble dpda_n = 0.0;
+        
+        intergralTerm = 0.0;
+        
+        for (int i = 2; i < numTrapeziums+1; i++)
+        {
+            A_n = A_n +dA;
+            getdpda(A_n, u, pacons, dpda_n);
+            intergralTerm = intergralTerm + sqrt(dpda_n/(m_rho*A_n));
+        }
+        
+        intergralTerm=intergralTerm*2.0;
+        getdpda(A_0, u, pacons, dpda_n);
+        intergralTerm=intergralTerm+sqrt(dpda_n/(m_rho*A_0));
+        getdpda(A, u, pacons, dpda_n);
+        intergralTerm=intergralTerm+sqrt(dpda_n/(m_rho*A));
+        intergralTerm = intergralTerm*0.5*(A-A_0)/numTrapeziums;
+        
+    }
+    
+    /**
+     *  Solves an Ax=b system using the linsys utility. Note that the A matrix must be passed as a 1D buffer array of
+     *  rows which is then used to form a NekMatrix
+     */
+    void ArterialPressureArea::v_solveAxb(int rows, const Array<OneD, NekDouble> &matrix_buf, const Array<OneD, NekDouble> &b_buf, Array<OneD, NekDouble> &x)
+    {
+        
+        //Cast A as NekMatrix and b as NekVector
+        boost::shared_ptr<NekMatrix<double> > N(new NekMatrix<double>(rows, rows, matrix_buf,eFULL));
+        boost::shared_ptr<NekVector<double> > b(new NekVector<double>(rows, b_buf));
+        
+        //Initialise Linear System
+        LinearSystem linsys(N);
+        
+        //Solve
+        NekVector<double> result = linsys.Solve(b);
+        result = linsys.SolveTranspose(b);
+        
+        //Convert the resulting NekVector to an array of NekDouble's for further manipulations later
+        for (int i = 0; i < rows; i++)
+        {
+            x[i]=result[i];
+        }
+        
     }
 
 }

solvers/PulseWaveSolver/EquationSystems/ArterialPressureArea.h

@@ -53,20 +53,132 @@ namespace Nektar
     public:
         /// Creates an instance of this class
       static PulseWavePressureAreaSharedPtr create(Array<OneD, MultiRegions::ExpListSharedPtr>& pVessel, 
-                                               const LibUtilities::SessionReaderSharedPtr& pSession)
+                                                   const LibUtilities::SessionReaderSharedPtr& pSession,
+                                                   const int &nDomains)
         {
-            return MemoryManager<ArterialPressureArea>::AllocateSharedPtr(pVessel,pSession);
+            return MemoryManager<ArterialPressureArea>::AllocateSharedPtr(pVessel,pSession,nDomains);
         }
 
         /// Name of class
         static std::string className;
         
         ArterialPressureArea(Array<OneD, MultiRegions::ExpListSharedPtr> pVessel, 
-                 const LibUtilities::SessionReaderSharedPtr pSession); 
+                             const LibUtilities::SessionReaderSharedPtr pSession,
+                             const int nDomains); 
 
         virtual ~ArterialPressureArea();
     protected:
-        virtual void v_DoPressure();
+
+        virtual void v_ReadParameters(int omega, int nqTrace, MultiRegions::ExpListSharedPtr &field);
+
+        virtual void v_GetPacons(int omega, Array<OneD, Array<OneD, NekDouble> > &pacons_trace, int nqTrace);
+
+        /// Pressure-area relationship definition
+        virtual void v_getPressure(
+            NekDouble A, 
+            Array<OneD, NekDouble> pacons, 
+            NekDouble &pressure);
+        
+        /// Definition of dp/da
+        virtual void v_getdpda(
+            NekDouble A, 
+            NekDouble u, 
+            Array<OneD, NekDouble> pacons, 
+            NekDouble &dpda);
+        
+        /// Calculates A from pressure
+        virtual void v_getAfromPressure(
+            NekDouble pressure, 
+            Array<OneD, NekDouble> pacons, 
+            NekDouble &A);
+        
+        /// Calculates the two characteristic variables
+        virtual void v_getW(
+            NekDouble A, 
+            NekDouble u, 
+            Array<OneD, NekDouble> pacons, 
+            NekDouble &W1, 
+            NekDouble &W2);
+        
+        /// Calculates A from W_1 and W_2
+        virtual void v_getAfromChars(
+            NekDouble W_1, 
+            NekDouble W_2, 
+            Array<OneD, NekDouble> pacons, 
+            NekDouble &A);
+        
+        /// Definition of u in terms of W_1 and W_2
+        virtual void v_getufromChars(
+            NekDouble W_1, 
+            NekDouble W_2, 
+            Array<OneD, NekDouble> pacons, 
+            NekDouble &u);
+        
+        /// lambda_1 relationship definition
+        virtual void v_getlambda1(
+            NekDouble A, 
+            NekDouble u, 
+            Array<OneD, NekDouble> pacons, 
+            NekDouble &lambda_1);
+        
+        /// lambda_2 relationship definition
+        virtual void v_getlambda2(
+            NekDouble A, 
+            NekDouble u, 
+            Array<OneD, NekDouble> pacons, 
+            NekDouble &lambda_2);
+        
+        /// Definition of dp/du
+        virtual void v_getdpdu(
+            NekDouble A, 
+            NekDouble u, 
+            Array<OneD, NekDouble> pacons, 
+            NekDouble &dpdu);
+        
+        /// Definition of dW1/da
+        virtual void v_getdW1da(
+            NekDouble A, 
+            NekDouble u, 
+            Array<OneD, NekDouble> pacons, 
+            NekDouble &dW1da);
+        
+        /// Definition of dW1/du
+        virtual void v_getdW1du(
+            NekDouble A, 
+            NekDouble u, 
+            Array<OneD, NekDouble> pacons, 
+            NekDouble &dW1du);
+        
+        /// Definition of dW2/dA
+        virtual void v_getdW2da(
+            NekDouble A, 
+            NekDouble u, 
+            Array<OneD, NekDouble> pacons, 
+            NekDouble &dW2da);
+        
+        /// Definition of dW2/du
+        virtual void v_getdW2du(
+            NekDouble A, 
+            NekDouble u, 
+            Array<OneD, NekDouble> pacons, 
+            NekDouble &dW2du);
+        
+        /// Evaluates the integral \int^A_{A_{0}}{ \sqrt{\frac{1}{\rho A} \frac{\partial p\left(A,t\right)}{\partial A}} dA}
+        virtual void v_getCharIntegral(
+            NekDouble A, 
+            NekDouble u, 
+            Array<OneD, NekDouble> pacons, 
+            NekDouble &intergralTerm);
+        
+        //void solveAxb(int rows, Array<OneD, double> A_buf, Array<OneD, double> b_buf);
+        virtual void v_solveAxb(
+            int rows, 
+            const Array<OneD, NekDouble> &matrix_buf, 
+            const Array<OneD, NekDouble> &b_buf, 
+            Array<OneD, NekDouble> &x);
+
+        Array<OneD, Array<OneD, NekDouble> > m_beta;
+        Array<OneD, Array<OneD, NekDouble> > m_beta_trace;
 
     private:
 

solvers/PulseWaveSolver/EquationSystems/LymphaticPressureArea.cpp

 ///////////////////////////////////////////////////////////////////////////////
 //
-// File CommMpi.cpp
+// File LymphaticPressureArea.cpp
 //
 // For more information, please see: http://www.nektar.info
 //
@@ -34,7 +34,7 @@
 ///////////////////////////////////////////////////////////////////////////////
 
 #include <PulseWaveSolver/EquationSystems/LymphaticPressureArea.h>
-
+#define PI 3.14159265359
 namespace Nektar
 {
 
@@ -48,9 +48,11 @@ namespace Nektar
      *
      */
     LymphaticPressureArea::LymphaticPressureArea(Array<OneD, MultiRegions::ExpListSharedPtr> pVessel, 
-                       const LibUtilities::SessionReaderSharedPtr pSession)
-        : PulseWavePressureArea(pVessel,pSession)
+                                                 const LibUtilities::SessionReaderSharedPtr pSession,
+                                                 const int nDomains)
+        : PulseWavePressureArea(pVessel,pSession,nDomains)
     {
+        m_session->LoadParameter("rho", m_rho, 0.5);
     }
 
     /**
@@ -61,9 +63,325 @@ namespace Nektar
 
     }
 
-    void LymphaticPressureArea::v_DoPressure()
+
+    void LymphaticPressureArea::v_ReadParameters(int omega, int nqTrace, MultiRegions::ExpListSharedPtr &field)
+    {
+    }
+
+    void LymphaticPressureArea::v_GetPacons(int omega, Array<OneD, Array<OneD, NekDouble> > &pacons_trace, int nqTrace)
+    {
+    }
+
+
+    /**
+     *  Contains the definition of the pressure-area relationship. To implement a new pressure-area relationship
+     *  change the mathematical definitions here and also in getdpda. Also consider changing: getAfromPressure
+     *  and getCharIntegral; currently use numerical methods but it is preferable to enter analyticall expressions
+     *  if they can be derived for the chosen p-A relation. Furthermore ensure that the definition of the lower
+     *  integration boundary A_0 is updated in getAfromPressure, getAfromChars and getCharIntegral if using the
+     *  numerical methods to ensure it refers to the correct element of pacons; typically can be assumed to be
+     *  equal to the passive area.
+     *
+     *  The p-A relation cannot be a function of u and can contain any number of parameters as specified by
+     *  P-A_num_Params in the input file which should then be specified for each subdomain using:
+     *  <FUNCTION NAME="P-A_Vals">
+     *  <E VAR="Val[Domain ID][Parameter ID]" VALUE="" />
+     *  </FUNCTION>
+     */
+    
+    void LymphaticPressureArea::v_getPressure(NekDouble A, Array<OneD, NekDouble> pacons, NekDouble &pressure)
+    {
+        
+        /*
+         * pacons is an array which holds, for the point for which at which we are evaluating this function,
+         * all of the constants found in the pressure-area relation. They are in the same order as the input file.
+         */
+        
+        //pacons[0] = p_{ext}
+        //pacons[1] = P_d
+        //pacons[2] = D_d
+        //pacons[3] = M
+        //pacons[4] = f
+        //pacons[5] = t_0
+        
+        NekDouble Passive = 0.0;
+        NekDouble Active = 0.0;
+        
+        Passive = pacons[1]*(exp(2*sqrt(A)/(pacons[2]*sqrt(PI))) -pow(pacons[2],3)*pow(PI,1.5)*pow(A,-1.5)/8);
+        Active = -(pacons[3]*sqrt(PI)/(2*sqrt(A)))*(cos(2*PI*pacons[4]*(m_time-pacons[5]))-1);
+        
+        pressure = pacons[0] + Passive + Active;
+        
+        //cout<<"Working?"<<pressure<<endl;
+        
+    }
+    
+    /**
+     *  Contains the definition of dp/da - needs to be updated when chaning p-A relation!
+     */
+    void LymphaticPressureArea::v_getdpda(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dpda)
+    {
+        NekDouble Passive = 0.0;
+        NekDouble Active = 0.0;
+        
+        Passive = pacons[1]*(exp(2*sqrt(A)/(pacons[2]*sqrt(PI)))/(sqrt(A*PI)*pacons[2]) + 3.0*pow(pacons[2],3)*pow(PI,1.5)*pow(A,-2.5)/16);
+        Active = (pacons[3]*sqrt(PI)/(4*pow(A,1.5)))*(cos(2*PI*pacons[4]*(m_time-pacons[5]))-1);
+        
+        dpda = Passive + Active;
+        
+    }
+    
+    /**
+     *  Calculates A in terms of pressure for for any p-A relation using a Newton iteration. If using a different
+     *  p-A relation the definition of A_0 may need modifying.
+     */
+    void LymphaticPressureArea::v_getAfromPressure(NekDouble pressure, Array<OneD, NekDouble> pacons, NekDouble &A)
+    {
+        
+        //A = (pressure - pacons[2])/pacons[0] + sqrt(pacons[1]);
+        //A = A*A;
+        
+        NekDouble fa = 0.0;
+        NekDouble dpda = 0.0;
+        NekDouble p_calc = 0.0;
+        NekDouble delta_A_calc = 0.0;
+        NekDouble A_0 = 0.0;
+        
+        // Initialise to A_0
+        A_0=PI*pow(pacons[2],2)/40;
+        A=A_0;
+        
+        int proceed = 1;
+        int iter = 0;
+        int MAX_ITER = 200;
+        
+        // Tolerances for the algorithm
+        NekDouble Tol = 1.0e-10;
+        
+        while ((proceed) && (iter < MAX_ITER))
+        {
+            iter =iter+1;
+            
+            getPressure(A, pacons, p_calc);
+            fa = p_calc-pressure;
+            getdpda(A, 0, pacons, dpda);
+            
+            delta_A_calc=fa/dpda;
+            A = A - delta_A_calc;
+            
+            if (sqrt(delta_A_calc*delta_A_calc) < Tol)
+                proceed = 0;
+        }
+        
+    }
+    
+    /**
+     *  Calculates the two characteristic variables for any p-A relation, currently using numerical
+     *  integration. Unlike lambda and the derivitives we calculate both in the same function to avoid having
+     *  to do the numerical integration twice.
+     */
+    void LymphaticPressureArea::v_getW(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &W1, NekDouble &W2)
+    {
+        NekDouble integralTerm=0.0;
+        getCharIntegral(A, u, pacons, integralTerm);
+        W1 = u + integralTerm;
+        W2 = u - integralTerm;
+    }
+    
+    /**
+     *  Contains the definition of A in terms of W_1 and W_2 for any p-A relation though definition of A_0 may
+     *  need updating. Solved for using a Newton iteration and numerical integration and numerical integration.
+     */
+    void LymphaticPressureArea::v_getAfromChars(NekDouble W_1, NekDouble W_2, Array<OneD, NekDouble> pacons, NekDouble &A)
+    {
+        
+        //    A = (W_1 - W_2)*sqrt(m_rho/(32*pacons[0]))+sqrt(sqrt(pacons[1]));
+        //    A = A*A;
+        //    A = A*A;
+        
+        NekDouble fa = 0.0;
+        NekDouble dfa = 0.0;
+        NekDouble dpda = 0.0;
+        NekDouble integralTerm = 0.0;
+        NekDouble delta_A_calc = 0.0;
+        NekDouble A_0 = 0.0;
+        
+        // Initialise to A_0
+        A_0=PI*pow(pacons[2],2)/40;
+        A=A_0;
+        
+        int proceed = 1;
+        int iter = 0;
+        int MAX_ITER = 200;
+        
+        // Tolerances for the algorithm
+        NekDouble Tol = 1.0e-10;
+        
+        while ((proceed) && (iter < MAX_ITER))
+        {
+            iter =iter+1;
+            
+            getCharIntegral(A, 0, pacons, integralTerm);
+            fa = integralTerm - 0.5*(W_1-W_2);
+            getdpda(A, 0, pacons, dpda);
+            dfa=sqrt(dpda/(m_rho*A));
+            
+            delta_A_calc=fa/dfa;
+            A = A - delta_A_calc;
+            
+            if (sqrt(delta_A_calc*delta_A_calc) < Tol)
+                proceed = 0;
+        }
+        
+    }
+    
+    /**
+     *  Contains the definition of u in terms of W_1 and W_2 for any p-A relation
+     */
+    void LymphaticPressureArea::v_getufromChars(NekDouble W_1, NekDouble W_2, Array<OneD, NekDouble> pacons, NekDouble &u)
+    {
+        
+        u = 0.5*(W_1+W_2);
+        
+    }
+    
+    /**
+     *  Contains the definition of the the first eigenvalue lambda_1 for any p-A relation
+     */
+    void LymphaticPressureArea::v_getlambda1(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &lambda_1)
+    {
+        NekDouble dpda=0.0;
+        getdpda(A, u, pacons, dpda);
+        
+        lambda_1 = u + sqrt((A/m_rho)*dpda);
+        
+    }
+    
+    /**
+     *  Contains the definition of the the second eigenvalue lambda_2 for any p-A relation
+     */
+    void LymphaticPressureArea::v_getlambda2(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &lambda_2)
+    {
+        NekDouble dpda=0.0;
+        getdpda(A, u, pacons, dpda);
+        
+        lambda_2 = u - sqrt((A/m_rho)*dpda);
+    }
+    
+    
+    /**
+     *  Contains the definition of dp/du for any p-A relation
+     */
+    void LymphaticPressureArea::v_getdpdu(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dpdu)
+    {
+        
+        dpdu = 0;
+        
+    }
+    
+    /**
+     *  Contains the definition of dW1/da for any p-A relation
+     */
+    void LymphaticPressureArea::v_getdW1da(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW1da)
+    {
+        NekDouble dpda=0.0;
+        getdpda(A, u, pacons, dpda);
+        
+        dW1da = sqrt(dpda/(A*m_rho));
+        
+    }
+    
+    /**
+     *  Contains the definition of dW1/du for any p-A relation
+     */
+    void LymphaticPressureArea::v_getdW1du(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW1du)
+    {
+        
+        dW1du = 1;
+        
+    }
+    
+    /**
+     *  Contains the definition of dW2/da for any p-A relation
+     */
+    void LymphaticPressureArea::v_getdW2da(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW2da)
+    {
+        
+        NekDouble dpda=0.0;
+        getdpda(A, u, pacons, dpda);
+        
+        dW2da = -1.0*sqrt(dpda/(A*m_rho));
+        
+    }
+    
+    /**
+     *  Contains the definition of dW2/du for any p-A relation
+     */
+    void LymphaticPressureArea::v_getdW2du(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW2du)
+    {
+        
+        dW2du = 1;
+        
+    }
+    
+    /**
+     *  Evalautes the integral \int^A_{A_{0}}{ \sqrt{\frac{1}{\rho A} \frac{\partial p\left(A,t\right)}{\partial A}} dA}
+     *  which is found in the general definition of the characteristic variables.
+     *  Currenly does this numerically but if p-A relation allows suggest the analytical expression is entered instead.
+     *  This current implementation is suitable for any p-A relation, though A_0 may need updating as well as num of
+     *  trapeziums.
+     */
+    void LymphaticPressureArea::v_getCharIntegral(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &intergralTerm)
+    {
+        
+        int numTrapeziums = 50;
+        NekDouble A_0 = 0.1*PI*pow(pacons[2],2)/4; //For active contractions can be less than static area
+        NekDouble A_n = A_0;
+        NekDouble dA = (A - A_n)/numTrapeziums;
+        NekDouble dpda_n = 0.0;
+        
+        intergralTerm = 0.0;
+        
+        for (int i = 2; i < numTrapeziums+1; i++)
+        {
+            A_n = A_n +dA;
+            getdpda(A_n, u, pacons, dpda_n);
+            intergralTerm = intergralTerm + sqrt(dpda_n/(m_rho*A_n));
+        }
+        
+        intergralTerm=intergralTerm*2.0;
+        getdpda(A_0, u, pacons, dpda_n);
+        intergralTerm=intergralTerm+sqrt(dpda_n/(m_rho*A_0));
+        getdpda(A, u, pacons, dpda_n);
+        intergralTerm=intergralTerm+sqrt(dpda_n/(m_rho*A));
+        intergralTerm = intergralTerm*0.5*(A-A_0)/numTrapeziums;
+        
+    }
+    
+    /**
+     *  Solves an Ax=b system using the linsys utility. Note that the A matrix must be passed as a 1D buffer array of
+     *  rows which is then used to form a NekMatrix
+     */
+    void LymphaticPressureArea::v_solveAxb(int rows, const Array<OneD, NekDouble> &matrix_buf, const Array<OneD, NekDouble> &b_buf, Array<OneD, NekDouble> &x)
     {
-        cout<<"lymphatic pressure area"<<endl; 
+        
+        //Cast A as NekMatrix and b as NekVector
+        boost::shared_ptr<NekMatrix<double> > N(new NekMatrix<double>(rows, rows, matrix_buf,eFULL));
+        boost::shared_ptr<NekVector<double> > b(new NekVector<double>(rows, b_buf));   
+        
+        //Initialise Linear System
+        LinearSystem linsys(N);
+        
+        //Solve
+        NekVector<double> result = linsys.Solve(b);
+        result = linsys.SolveTranspose(b);
+        
+        //Convert the resulting NekVector to an array of NekDouble's for further manipulations later
+        for (int i = 0; i < rows; i++)
+        {
+            x[i]=result[i];
+        }
+        
     }
 
 }

solvers/PulseWaveSolver/EquationSystems/LymphaticPressureArea.h

@@ -52,22 +52,70 @@ namespace Nektar
     {
     public:
         /// Creates an instance of this class
-      static PulseWavePressureAreaSharedPtr create(Array<OneD, MultiRegions::ExpListSharedPtr>& pVessel, 
-                                               const LibUtilities::SessionReaderSharedPtr& pSession)
+        static PulseWavePressureAreaSharedPtr create(Array<OneD, MultiRegions::ExpListSharedPtr>& pVessel, 
+                                                     const LibUtilities::SessionReaderSharedPtr& pSession,
+                                                     const int &nDomains)
         {
-            return MemoryManager<LymphaticPressureArea>::AllocateSharedPtr(pVessel,pSession);
+            return MemoryManager<LymphaticPressureArea>::AllocateSharedPtr(pVessel,pSession,nDomains);
         }
 
         /// Name of class
         static std::string className;
         
         LymphaticPressureArea(Array<OneD, MultiRegions::ExpListSharedPtr> pVessel, 
-                 const LibUtilities::SessionReaderSharedPtr pSession); 
+                              const LibUtilities::SessionReaderSharedPtr pSession,
+                              const int nDomains); 
 
         virtual ~LymphaticPressureArea();
     protected:
-        virtual void v_DoPressure();
+        virtual void v_ReadParameters(int omega, int nqTrace, MultiRegions::ExpListSharedPtr &field);
 
+        virtual void v_GetPacons(int omega, Array<OneD, Array<OneD, NekDouble> > &pacons_trace, int nqTrace);
+
+        /// Pressure-area relationship definition
+        virtual void v_getPressure(NekDouble A, Array<OneD, NekDouble> pacons, NekDouble &pressure);
+        
+        /// Definition of dp/da
+        virtual void v_getdpda(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dpda);
+        
+        /// Calculates A from pressure
+        virtual void v_getAfromPressure(NekDouble pressure, Array<OneD, NekDouble> pacons, NekDouble &A);
+        
+        /// Calculates the two characteristic variables
+        virtual void v_getW(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &W1, NekDouble &W2);
+        
+        /// Calculates A from W_1 and W_2
+        virtual void v_getAfromChars(NekDouble W_1, NekDouble W_2, Array<OneD, NekDouble> pacons, NekDouble &A);
+        
+        /// Definition of u in terms of W_1 and W_2
+        virtual void v_getufromChars(NekDouble W_1, NekDouble W_2, Array<OneD, NekDouble> pacons, NekDouble &u);
+        
+        /// lambda_1 relationship definition
+        virtual void v_getlambda1(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &lambda_1);
+        
+        /// lambda_2 relationship definition
+        virtual void v_getlambda2(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &lambda_2);
+        
+        /// Definition of dp/du
+        virtual void v_getdpdu(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dpdu);
+        
+        /// Definition of dW1/da
+        virtual void v_getdW1da(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW1da);
+        
+        /// Definition of dW1/du
+        virtual void v_getdW1du(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW1du);
+        
+        /// Definition of dW2/dA
+        virtual void v_getdW2da(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW2da);
+        
+        /// Definition of dW2/du
+        virtual void v_getdW2du(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW2du);
+        
+        /// Evaluates the integral \int^A_{A_{0}}{ \sqrt{\frac{1}{\rho A} \frac{\partial p\left(A,t\right)}{\partial A}} dA}
+        virtual void v_getCharIntegral(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &intergralTerm);
+        
+        //void solveAxb(int rows, Array<OneD, double> A_buf, Array<OneD, double> b_buf);
+        virtual void v_solveAxb(int rows, const Array<OneD, NekDouble> &matrix_buf, const Array<OneD, NekDouble> &b_buf, Array<OneD, NekDouble> &x);
     private:
 
     };

solvers/PulseWaveSolver/EquationSystems/PulseWavePressureArea.cpp

@@ -45,9 +45,11 @@ namespace Nektar
      */
   PulseWavePressureArea::PulseWavePressureArea(
       Array<OneD, MultiRegions::ExpListSharedPtr> &pVessel, 
-      const LibUtilities::SessionReaderSharedPtr &pSession)
+      const LibUtilities::SessionReaderSharedPtr &pSession,
+      const int &nDomains)
     : m_vessels(pVessel),
-      m_session(pSession)
+      m_session(pSession),
+      m_domains(nDomains)
     {
     }
 
@@ -66,4 +68,98 @@ namespace Nektar
                                       Loki::NoDestroy > Type;
         return Type::Instance();
     }
+
+    void PulseWavePressureArea::EvaluateFunction(
+        Array<OneD, MultiRegions::ExpListSharedPtr>       pFields,
+        LibUtilities::SessionReaderSharedPtr              pSession,
+        std::string                                       pFieldName,
+        Array<OneD, NekDouble>&                           pArray,
+        const std::string&                                pFunctionName,
+        NekDouble                                         pTime,
+        const int domain,
+        const int nq)
+    {
+        ASSERTL0(pSession->DefinesFunction(pFunctionName),
+                 "Function '" + pFunctionName + "' does not exist.");
+
+        //unsigned int nq = pFields[0]->GetNpoints();
+        //if (pArray.num_elements() != nq)
+        //{
+        //pArray = Array<OneD, NekDouble> (nq);
+        //}
+
+        LibUtilities::FunctionType vType;
+        vType = pSession->GetFunctionType(pFunctionName, pFieldName);
+        if (vType == LibUtilities::eFunctionTypeExpression)
+        {
+            Array<OneD, NekDouble> x0(nq);
+            Array<OneD, NekDouble> x1(nq);
+            Array<OneD, NekDouble> x2(nq);
+
+            pFields[0]->GetCoords(x0, x1, x2);
+            LibUtilities::EquationSharedPtr ffunc =
+                pSession->GetFunction(pFunctionName, pFieldName,domain);
+
+            ffunc->Evaluate(x0, x1, x2, pTime, pArray);
+        }
+        else if (vType == LibUtilities::eFunctionTypeFile)
+        {
+            std::string filename = pSession->GetFunctionFilename(
+                pFunctionName,
+                pFieldName,
+                domain);
+
+            std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
+            std::vector<std::vector<NekDouble> > FieldData;
+            Array<OneD, NekDouble> vCoeffs(pFields[0]->GetNcoeffs());
+            Vmath::Zero(vCoeffs.num_elements(), vCoeffs, 1);
+
+            int numexp = pFields[0]->GetExpSize();
+            Array<OneD,int> ElementGIDs(numexp);
+
+            // Define list of global element ids
+            for(int i = 0; i < numexp; ++i)
+            {
+                ElementGIDs[i] = pFields[0]->GetExp(i)->GetGeom()->GetGlobalID();
+            }
+
+            LibUtilities::FieldIOSharedPtr fld =
+                MemoryManager<LibUtilities::FieldIO>::AllocateSharedPtr(m_session->GetComm());
+
+            fld->Import(filename,FieldDef,FieldData,
+                          LibUtilities::NullFieldMetaDataMap,
+                          ElementGIDs);
+            int idx = -1;
+
+            // Loop over all the expansions
+            for (int i = 0; i < FieldDef.size(); ++i)
+            {
+                // Find the index of the required field in the
+                // expansion segment
+                for(int j = 0; j < FieldDef[i]->m_fields.size(); ++j)
+                {
+                    if (FieldDef[i]->m_fields[j] == pFieldName)
+                    {
+                        idx = j;
+                    }
+                }
+
+                if (idx >= 0)
+                {
+                    pFields[0]->ExtractDataToCoeffs(
+                        FieldDef[i], FieldData[i],
+                        FieldDef[i]->m_fields[idx], vCoeffs);
+                }
+                else
+                {
+                    cout << "Field " + pFieldName + " not found." << endl;
+                }
+            }
+
+
+            pFields[0]->BwdTrans_IterPerExp(vCoeffs, pArray);
+
+        }
+    }
+
 }

solvers/PulseWaveSolver/EquationSystems/PulseWavePressureArea.h

@@ -51,37 +51,145 @@ namespace Nektar
     typedef LibUtilities::NekFactory< std::string, 
         PulseWavePressureArea, 
         Array<OneD, MultiRegions::ExpListSharedPtr>&, 
-        const LibUtilities::SessionReaderSharedPtr& > PressureAreaFactory;
+        const LibUtilities::SessionReaderSharedPtr&,
+        const int& > PressureAreaFactory;
     PressureAreaFactory& GetPressureAreaFactory();
-    
+   
     class PulseWavePressureArea
     {
     public:
         PulseWavePressureArea(Array<OneD, MultiRegions::ExpListSharedPtr> &pVessel,
-                          const LibUtilities::SessionReaderSharedPtr &pSession);
+                              const LibUtilities::SessionReaderSharedPtr &pSession,
+                              const int &nDomains);
 
         virtual ~PulseWavePressureArea();
 
-        inline void DoPressure();
+        inline void ReadParameters(int omega, int nqTrace, MultiRegions::ExpListSharedPtr &field);
+        inline void GetPacons(int omega, Array<OneD, Array<OneD, NekDouble> > &pacons_trace, int nqTrace);
+        inline void getPressure(NekDouble A, Array<OneD, NekDouble> pacons, NekDouble &pressure);
+        inline void getdpda(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dpda);
+        inline void getAfromPressure(NekDouble pressure, Array<OneD, NekDouble> pacons, NekDouble &A);
+        inline void getW(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &W1, NekDouble &W2);
+        inline void getAfromChars(NekDouble W_1, NekDouble W_2, Array<OneD, NekDouble> pacons, NekDouble &A);
+        inline void getufromChars(NekDouble W_1, NekDouble W_2, Array<OneD, NekDouble> pacons, NekDouble &u);
+        inline void getlambda1(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &lambda_1);
+        inline void getlambda2(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &lambda_2);
+        inline void getdpdu(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dpdu);
+        inline void getdW1da(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW1da);
+        inline void getdW1du(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW1du);
+        inline void getdW2da(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW2da);
+        inline void getdW2du(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW2du);
+        inline void getCharIntegral(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &intergralTerm);
+        inline void solveAxb(int rows, const Array<OneD, NekDouble> &matrix_buf, const Array<OneD, NekDouble> &b_buf, Array<OneD, NekDouble> &x);
 
     protected:
-        virtual void v_DoPressure()=0;
+        virtual void v_ReadParameters(int omega, int nqTrace, MultiRegions::ExpListSharedPtr &field)=0;
+        virtual void v_GetPacons(int omega, Array<OneD, Array<OneD, NekDouble> > &pacons_trace, int nqTrace)=0;
+        virtual void v_getPressure(NekDouble A, Array<OneD, NekDouble> pacons, NekDouble &pressure)=0;
+        virtual void v_getdpda(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dpda)=0;
+        virtual void v_getAfromPressure(NekDouble pressure, Array<OneD, NekDouble> pacons, NekDouble &A)=0;
+        virtual void v_getW(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &W1, NekDouble &W2)=0;
+        virtual void v_getAfromChars(NekDouble W_1, NekDouble W_2, Array<OneD, NekDouble> pacons, NekDouble &A)=0;
+        virtual void v_getufromChars(NekDouble W_1, NekDouble W_2, Array<OneD, NekDouble> pacons, NekDouble &u)=0;
+        virtual void v_getlambda1(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &lambda_1)=0;
+        virtual void v_getlambda2(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &lambda_2)=0;
+        virtual void v_getdpdu(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dpdu)=0;
+        virtual void v_getdW1da(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW1da)=0;
+        virtual void v_getdW1du(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW1du)=0;
+        virtual void v_getdW2da(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW2da)=0;
+        virtual void v_getdW2du(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW2du)=0;
+        virtual void v_getCharIntegral(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &intergralTerm)=0;
+        virtual void v_solveAxb(int rows, const Array<OneD, NekDouble> &matrix_buf, const Array<OneD, NekDouble> &b_buf, Array<OneD, NekDouble> &x)=0;
 
         Array<OneD, MultiRegions::ExpListSharedPtr> m_vessels;
 	LibUtilities::SessionReaderSharedPtr m_session;
+        const int m_domains;
+
+        NekDouble m_rho;
+        NekDouble m_time;
+
+        void EvaluateFunction(
+            Array<OneD, MultiRegions::ExpListSharedPtr> pFields,
+            LibUtilities::SessionReaderSharedPtr pSession,
+            std::string pFieldName, 
+            Array<OneD, NekDouble>& pArray,
+            const std::string& pFunctionName,
+            NekDouble pTime = NekDouble(0),
+            const int domain = 0,
+            const int nq=0);
 
     private:
     };
 
-    /**
-     *
-     */
-    inline void PulseWavePressureArea::DoPressure()
+    inline void PulseWavePressureArea::ReadParameters(int omega, int nqTrace, MultiRegions::ExpListSharedPtr &field)
     {
-        v_DoPressure();
+        v_ReadParameters(omega, nqTrace, field);
     }
 
+    inline void PulseWavePressureArea::GetPacons(int omega, Array<OneD, Array<OneD, NekDouble> > &pacons_trace, int nqTrace)
+    {
+        v_GetPacons(omega, pacons_trace, nqTrace);
+    }
 
-
+    inline void PulseWavePressureArea::getPressure(NekDouble A, Array<OneD, NekDouble> pacons, NekDouble &pressure)
+    {
+        v_getPressure(A, pacons, pressure);
+    };
+    inline void PulseWavePressureArea::getdpda(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dpda)
+    {
+        v_getdpda(A, u, pacons, dpda);
+    };
+    inline void PulseWavePressureArea::getAfromPressure(NekDouble pressure, Array<OneD, NekDouble> pacons, NekDouble &A)
+    {
+        v_getAfromPressure(pressure, pacons, A);
+    };
+    inline void PulseWavePressureArea::getW(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &W1, NekDouble &W2)
+    {
+        v_getW(A, u, pacons, W1, W2);
+    };
+    inline void PulseWavePressureArea::getAfromChars(NekDouble W_1, NekDouble W_2, Array<OneD, NekDouble> pacons, NekDouble &A)
+    {
+        v_getAfromChars(W_1, W_2, pacons, A);
+    };
+    inline void PulseWavePressureArea::getufromChars(NekDouble W_1, NekDouble W_2, Array<OneD, NekDouble> pacons, NekDouble &u)
+    {
+        v_getufromChars(W_1, W_2, pacons, u);
+    };
+    inline void PulseWavePressureArea::getlambda1(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &lambda_1)
+    {
+        v_getlambda1(A, u, pacons, lambda_1);
+    };
+    inline void PulseWavePressureArea::getlambda2(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &lambda_2)
+    {
+        v_getlambda2(A, u, pacons, lambda_2);
+    };
+    inline void PulseWavePressureArea::getdpdu(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dpdu)
+    {
+        v_getdpdu(A, u, pacons, dpdu);
+    };
+    inline void PulseWavePressureArea::getdW1da(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW1da)
+    {
+        v_getdW1da(A, u, pacons, dW1da);
+    };
+    inline void PulseWavePressureArea::getdW1du(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW1du)
+    {
+        v_getdW1du(A, u, pacons, dW1du);
+    };
+    inline void PulseWavePressureArea::getdW2da(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW2da)
+    {
+        v_getdW2da(A, u, pacons, dW2da);
+    };
+    inline void PulseWavePressureArea::getdW2du(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &dW2du)
+    {
+        v_getdW2du(A, u, pacons, dW2du);
+    };
+    inline void PulseWavePressureArea::getCharIntegral(NekDouble A, NekDouble u, Array<OneD, NekDouble> pacons, NekDouble &intergralTerm)
+    {
+        v_getCharIntegral(A, u, pacons, intergralTerm);
+    };
+    inline void PulseWavePressureArea::solveAxb(int rows, const Array<OneD, NekDouble> &matrix_buf, const Array<OneD, NekDouble> &b_buf, Array<OneD, NekDouble> &x)
+    {
+        v_solveAxb(rows, matrix_buf, b_buf, x);
+    };
 }
 #endif

solvers/PulseWaveSolver/EquationSystems/PulseWavePropagation.cpp

@@ -61,9 +61,6 @@ namespace Nektar
     {
         PulseWaveSystem::v_InitObject();
 	
-        m_pressureArea=GetPressureAreaFactory().CreateInstance("Lymphatic",m_vessels,m_session);
-        m_pressureArea->DoPressure();
-	
         if (m_explicitAdvection)
         {
             m_ode.DefineOdeRhs       (&PulseWavePropagation::DoOdeRhs, this);
@@ -112,6 +109,13 @@ namespace Nektar
         // do advection evauation in all domains
         for(int omega=0; omega < m_nDomains; ++omega)
         {
+            m_pressureArea->GetPacons(omega, m_pacons_trace, GetTraceTotPoints());
+
+            for(int i=0; i<m_pacons_trace[0].num_elements(); ++i)
+            {
+                cout<<"m_pacons_trace[0][i]: "<<m_pacons_trace[0][i]<<endl;
+            }
+
             m_currentDomain = omega;
             int nq = m_vessels[omega*m_nVariables]->GetTotPoints();
             int ncoeffs = m_vessels[omega*m_nVariables]->GetNcoeffs();

solvers/PulseWaveSolver/EquationSystems/PulseWavePropagation.h

@@ -38,7 +38,7 @@
 
 #include <PulseWaveSolver/EquationSystems/PulseWaveSystem.h>
 #include <PulseWaveSolver/EquationSystems/PulseWaveBoundary.h>
-#include <PulseWaveSolver/EquationSystems/PulseWavePressureArea.h>
+
 
 using namespace Nektar::SolverUtils;
 
@@ -93,8 +93,8 @@ namespace Nektar
                                  NekDouble n);
         
         Array<OneD, PulseWaveBoundarySharedPtr> m_Boundary;
+        Array<OneD, Array<OneD, NekDouble> > m_pacons_trace;
 
-        PulseWavePressureAreaSharedPtr m_pressureArea;
         virtual void v_GenerateSummary(SolverUtils::SummaryList& s);
     };
 }

solvers/PulseWaveSolver/EquationSystems/PulseWaveSystem.cpp

@@ -196,6 +196,15 @@ namespace Nektar
         // Load external pressure
         m_session->LoadParameter("pext", m_pext, 0.0);
         
+        std::string vPressureArea = "Arterial";
+
+        if (m_session->DefinesSolverInfo("PressureArea"))
+        {
+            vPressureArea = m_session->GetSolverInfo("PressureArea");
+        }
+
+        m_pressureArea=GetPressureAreaFactory().CreateInstance(vPressureArea,m_vessels,m_session,m_nDomains);
+
         int nq = 0;
         /**
          *  Gets the Material Properties of each arterial segment
@@ -226,6 +235,8 @@ namespace Nektar
 
             int nqTrace = GetTraceTotPoints();
 
+            m_pressureArea->ReadParameters(omega,nqTrace,m_fields[0]);
+
             m_beta_trace[omega] = Array<OneD, NekDouble>(nqTrace);
             m_fields[0]->ExtractTracePhys(m_beta[omega],m_beta_trace[omega]);
             

solvers/PulseWaveSolver/EquationSystems/PulseWaveSystem.h

@@ -37,6 +37,7 @@
 #define NEKTAR_SOLVERS_PULSEWAVESOLVER_EQUATIONSYSTEMS_PULSEWAVESYSTEM_H
 
 #include <SolverUtils/UnsteadySystem.h>
+#include <PulseWaveSolver/EquationSystems/PulseWavePressureArea.h>
 
 using namespace Nektar::SolverUtils;
 
@@ -172,6 +173,8 @@ namespace Nektar
         void WriteVessels(const std::string &outname);
         
         void EnforceInterfaceConditions(const Array<OneD, const Array<OneD, NekDouble> > &fields);
+
+        PulseWavePressureAreaSharedPtr m_pressureArea;
         
     private:
         void SetUpDomainInterfaces(void);