FEBio/FETriphasic.h

#pragma once
#include "FEBiphasicSolute.h"
#include "FESolute.h"

//-----------------------------------------------------------------------------
// Monovalent salt (cation + anion)

class FESaltMaterialPoint : public FEMaterialPoint
{
public:
        FESaltMaterialPoint(FEMaterialPoint* ppt) : FEMaterialPoint(ppt) {}
        
        FEMaterialPoint* Copy()
        {
                FESaltMaterialPoint* pt = new FESaltMaterialPoint(*this);
                if (m_pt) pt->m_pt = m_pt->Copy();
                return pt;
        }
        
        void Serialize(DumpFile& ar)
        {
                if (ar.IsSaving())
                {
                        ar << m_c[0] << m_gradc[0] << m_j[0] << m_ca[0];
                        ar << m_c[1] << m_gradc[1] << m_j[1] << m_ca[1];
                }
                else
                {
                        ar >> m_c[0] >> m_gradc[0] >> m_j[0] >> m_ca[0];
                        ar >> m_c[1] >> m_gradc[1] >> m_j[1] >> m_ca[1];
                }
                
                if (m_pt) m_pt->Serialize(ar);
        }
        
        void Init(bool bflag)
        {
                if (bflag)
                {
                        m_c[0] = m_c[1] = 0;
                        m_ca[0] = m_ca[1] = 0;
                        m_psi = m_cF = 0;
                        m_gradc[0] = m_gradc[1] = vec3d(0,0,0);
                        m_j[0] = m_j[1] = m_Ie = vec3d(0,0,0);
                }
                
                if (m_pt) m_pt->Init(bflag);
        }
        
public:
        // salt material data
        double          m_c[2];         
        vec3d           m_gradc[2];     
        vec3d           m_j[2];         
        double          m_ca[2];        
        double          m_psi;          
        vec3d           m_Ie;           
        double          m_cF;           
};

//-----------------------------------------------------------------------------

class FETriphasic : public FEMaterial
{
public:
        FETriphasic();
        
        // returns a pointer to a new material point object
        FEMaterialPoint* CreateMaterialPointData() 
        { 
                return new FESaltMaterialPoint(new FEBiphasicMaterialPoint(m_pSolid->CreateMaterialPointData()));
        }

        // return elastic material component
        FEElasticMaterial* GetElasticMaterial() { return m_pSolid; }
        
public:
        void Init();
        
        mat3ds Stress(FEMaterialPoint& pt);
        
        tens4ds Tangent(FEMaterialPoint& pt);
        
        vec3d FluidFlux(FEMaterialPoint& pt);
        
        vec3d SoluteFlux(FEMaterialPoint& pt, const int ion);
        
        double Pressure(FEMaterialPoint& pt);
        
        double Concentration(FEMaterialPoint& pt, const int ion);
        
        double Porosity(FEMaterialPoint& pt);
        
        double FixedChargeDensity(FEMaterialPoint& pt);
        
        double ElectricPotential(FEMaterialPoint& pt, const bool eform=false);
        
        vec3d CurrentDensity(FEMaterialPoint& pt);
        
        double FluidDensity() { return m_rhoTw; }
        
        double SoluteDensity(const int ion) { return m_pSolute[ion]->Density(); }
        
        double SoluteMolarMass(const int ion) { return m_pSolute[ion]->MolarMass(); }
        
        int SoluteChargeNumber(const int ion) { return m_pSolute[ion]->ChargeNumber(); }
        
        void Serialize(DumpFile& ar);
        
public:
        double                                          m_phi0;                 
        FEElasticMaterial*                      m_pSolid;               
        double                                          m_rhoTw;                
        FEHydraulicPermeability*        m_pPerm;                
        FEOsmoticCoefficient*           m_pOsmC;                
        vector<FESolute*>                       m_pSolute;              
        double                                          m_cFr;                  
        double                                          m_Rgas;                 
        double                                          m_Tabs;                 
        double                                          m_Fc;                   
        double                                          m_penalty;              
        
        // declare as registered
        DECLARE_REGISTERED(FETriphasic);
        
        DECLARE_PARAMETER_LIST();
};