flowelem.h

/*************************************************************************************
 * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
 * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
 * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
 *                                                                                   *
 * This file is part of MechSys.                                                     *
 *                                                                                   *
 * MechSys is free software; you can redistribute it and/or modify it under the      *
 * terms of the GNU General Public License as published by the Free Software         *
 * Foundation; either version 2 of the License, or (at your option) any later        *
 * version.                                                                          *
 *                                                                                   *
 * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY        *
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   *
 * PARTICULAR PURPOSE. See the GNU General Public License for more details.          *
 *                                                                                   *
 * You should have received a copy of the GNU General Public License along with      *
 * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
 * Fifth Floor, Boston, MA 02110-1301, USA                                           *
 *************************************************************************************/

#ifndef FEM_FLOWELEM_H
#define FEM_FLOWELEM_H

#ifndef REAL
    #define REAL double
#endif

#include <sstream>

#include "models/flow/flowmodel.h"
#include "fem/ele/element.h"
#include "tensors/tensors.h"
#include "tensors/functions.h" 
#include "util/numstreams.h"

namespace FEM
{

class FlowElem: public virtual Element
{
public:
    // FlowElem constants {{{   
    static int    NDIM;
    static String DWP;
    static String DWD;
    static String DFWD;
    // }}}

    // Destructor
    virtual ~FlowElem() {}
    // Methods
    virtual void ReAllocateModel(String const & ModelName, Array<REAL> const & ModelPrms, Array<Array<REAL> > const & AIniData);
    virtual bool IsEssential(String const & DOFName) const;
    virtual void CalcFaceNodalValues(String            const & FaceDOFName  , // In:  Face DOF Name, ex.: "Dfwd"
                             REAL              const   FaceDOFValue ,         // In:  Face DOF Value, ex.: 10 m^3
                             Array<FEM::Node*> const & APtrFaceNodes,         // In:  Array of ptrs to face nodes
                             String                  & NodalDOFName ,         // Out: Nodal DOF Name, ex.: "Dwd"
                             LinAlg::Vector<REAL>    & NodalValues            // Out: Vector of nodal values
    ) const;
    virtual void Bp_Matrix(LinAlg::Matrix<REAL> const & derivs, LinAlg::Matrix<REAL> const & J, LinAlg::Matrix<REAL> & Bp) const;
    virtual void Permeability(LinAlg::Matrix<REAL> & H, Array<size_t> & EqMap) const; 
    // Create a structure to hold the DOFs information of Node [j_node] inside Element [i_ele]
    // These DOFs are the only ones that Element [i_ele] uses => NOT the global ones
    virtual void NodalDOFs(int iNode, Array<FEM::Node::DOFVarsStruct*> & DOFs) const;
    virtual void BackupState();
    virtual void UpdateState(REAL TimeInc); // Read nodal values _IncEssential and write to _IncNatural
    virtual void RestoreState();
    virtual void SetProperties(Array<REAL> const & EleProps) { }
    virtual String OutCenter(bool PrintCaptionOnly) const;
    virtual void OutNodes(LinAlg::Matrix<REAL> & Values, Array<String> & Labels) const;
    // Functions to assemble DAS matrices
    size_t nOrder0Matrices() const { return 1; };
    void      RHSVector(size_t index, LinAlg::Vector<REAL> & V, Array<size_t> & Map) const;
    void   Order0Matrix(size_t index, LinAlg::Matrix<REAL> & V, Array<size_t> & RowsMap, Array<size_t> & ColsMap) const;
private:
    // Data
    Array<FlowModel*> _a_model;
    // Methods
    virtual void _set_node_vars(int iNode);
    virtual void _calc_Nu_matrix(LinAlg::Vector<REAL> & Shape, LinAlg::Matrix<REAL> & Nu) const;
    virtual void _calc_initial_internal_pore_pressures();
}; // class FlowElem

// FlowElem constants {{{   
int    FlowElem::NDIM = 3;
String FlowElem::DWP  = _T("Dwp");
String FlowElem::DWD  = _T("Dwd");
String FlowElem::DFWD  = _T("Dfwd");
// }}}

// ========================================================================== Implementation

inline void FlowElem::_set_node_vars(int iNode) // {{{
{
    // Add Degree of Freedom to a node (Essential, Natural) 
    _connects[iNode]->AddDOF(DWP,DWD);
} // }}}

void FlowElem::ReAllocateModel(String const & ModelName, Array<REAL> const & ModelPrms, Array<Array<REAL> > const & AIniData) // {{{
{
    // Check size of AIniData
    if (!(AIniData.size()==1 || static_cast<int>(AIniData.size())==_n_int_pts))
        throw new Fatal(_("FlowElem::ReAllocateModel: Array of array of initial data must have size==1 or size== %d \n"),_n_int_pts);

    // If pointers to model was not already defined => No model was allocated
    if (_a_model.size()==0)
    {
        // Resize the array of model pointers
        _a_model.resize(_n_int_pts);

        // Loop along integration points
        for (int i=0; i<_n_int_pts; ++i)
        {
            // Allocate a new model and set parameters
            if (AIniData.size()==1) _a_model[i] = AllocFlowModel(ModelName, ModelPrms, AIniData[0]);
            else                    _a_model[i] = AllocFlowModel(ModelName, ModelPrms, AIniData[i]);
        }

        // Calculate initial internal pore-pressures
        _calc_initial_internal_pore_pressures();
    }
    else // Model objects already defined (allocated)
    {
        // Loop along integration points
        for (int i=0; i<_n_int_pts; ++i)
        {
            if (_a_model[i]->Name()!=ModelName) // Do re-allocate
            {
                // Delete old model
                delete _a_model[i];

                // Allocate a new model and set parameters
                if (AIniData.size()==1) _a_model[i] = AllocFlowModel(ModelName, ModelPrms, AIniData[0]);
                else                    _a_model[i] = AllocFlowModel(ModelName, ModelPrms, AIniData[i]);
            }
            // else: Do NOT re-allocate => Model not changed (don't do anything)
        }
    }
} // }}}

inline bool FlowElem::IsEssential(String const & DOFName) const // {{{
{
    if (DOFName==DWP) return true;
    else return false;
} // }}}

inline void FlowElem::CalcFaceNodalValues(String            const & FaceDOFName ,         // In:  Face DOF Name, ex.: "tx, ty, tz" {{{
                                         REAL              const   FaceDOFValue ,         // In:  Face DOF Value, ex.: 10 kPa, 20 kPa
                                         Array<FEM::Node*> const & APtrFaceNodes,         // In:  Array of ptrs to face nodes
                                         String                  & NodalDOFName ,         // Out: Nodal DOF Name, ex.: "fx, fy, fz"
                                         LinAlg::Vector<REAL>    & NodalValues  ) const   // Out: Vector of nodal values
{
    if (FaceDOFName==DFWD)
    {
        if (FaceDOFName==DFWD) NodalDOFName=DWD;
        _distrib_val_to_face_nodal_vals(APtrFaceNodes, FaceDOFValue, NodalValues);
    }
    else
    {
        std::ostringstream oss; oss << "Face nodes coordinates:\n";
        for (size_t i_node=0; i_node<APtrFaceNodes.size(); ++i_node)
            oss << "X=" << APtrFaceNodes[i_node]->X() << ", Y=" << APtrFaceNodes[i_node]->Y() << ", Z=" << APtrFaceNodes[i_node]->Z() << std::endl;
        throw new Fatal(_("FlowElem::CalcFaceNodalValues: This method must only be called for FaceDOFName< %d > equal to tq\n %s \n"), FaceDOFName.c_str(), oss.str().c_str());
    }
} // }}}

inline void FlowElem::Permeability(LinAlg::Matrix<REAL> & He, Array<size_t> & EqMap) const // {{{
{
    //  
    //   Permeability Matrix K:
    //   ============================
    //       
    //                   /    T                   
    //         [P]   =   | [B]  * [K] * {B}  * dV
    //                   /    p            p  
    //           
    //   OBS.: [K] = [k]/gamaW
    
    // Resize He
    He.Resize(_n_nodes, _n_nodes); // sum(Bpt*k*Bp*det(J)*w)
    He.SetValues(0.0);

    // Loop along integration points
    for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
    {
        // Temporary Integration Points
        REAL r = _a_int_pts[i_ip].r;
        REAL s = _a_int_pts[i_ip].s;
        REAL t = _a_int_pts[i_ip].t;
        REAL w = _a_int_pts[i_ip].w;

        // Calculate Derivatives of Shape functions w.r.t local coordinate system
        LinAlg::Matrix<REAL> derivs; // size = NumLocalCoords(ex.: r,s,t) x _n_nodes
        Derivs(r,s,t, derivs);

        // Calculate J (Jacobian) matrix for i_ip Integration Point
        LinAlg::Matrix<REAL> J;
        Jacobian(derivs, J);
        REAL det_J = J.Det();

        // Calculate B matrix for i_ip Integration Point
        LinAlg::Matrix<REAL> Bp;
        Bp_Matrix(derivs,J, Bp);

        // Calculate Tangent Permeability
        blitz::TinyVector<REAL,6> K;
        _a_model[i_ip]->TgPermeability(K);

        REAL gammaW = _a_model[i_ip]->GammaW();

        // Calculate K*Bp
        LinAlg::Matrix<REAL> KBp(Bp.Rows(),Bp.Cols());
        for (int i=0; i<Bp.Rows(); ++i)
        for (int j=0; j<Bp.Cols(); ++j)
            KBp(i,j) = K(i)*Bp(i,j); //this multiplication consider that permeability tensor is diagonal!

        // Calculate Bpt*K*Bp and find He // [BptKBp]_{Bp_cols x Bp_cols}
        for (int i=0; i<Bp.Cols(); ++i)
        for (int j=0; j<Bp.Cols(); ++j)
        {
            for (int k=0; k<Bp.Rows(); ++k)
                He(i,j) += -Bp(k,i)*KBp(k,j)*1.0/gammaW*det_J*w; //negative signal value because Darcy's law
        }
    }
    
    //Mounting a map of positions from He to Global
    
    int idx_He=0; // position (idx) inside He matrix
    EqMap.resize(He.Rows()); // size=He.Rows()=He.Cols()
    for (int i_node=0; i_node<_n_nodes; ++i_node)
    {
        // Fill map of He position to Global position of DOFs components
        EqMap[idx_He++] = _connects[i_node]->DOFVar(DWP).EqID; 
    }
} // }}}

inline void FlowElem::_calc_Nu_matrix(LinAlg::Vector<REAL> & Shape, //{{{
                                    LinAlg::Matrix<REAL> & Nu//Nu must be (3 x 3*_n_nodes);
                                    ) const
{
    Nu.Resize(NDIM, NDIM*_n_nodes);
    for(int i=0; i<_n_nodes; i++)
    {
        Nu(0,0 + 3*i) = Shape(i);  Nu(0,1 + 3*i) = 0;         Nu(0,2 + 3*i) = 0;    
        Nu(1,0 + 3*i) = 0;         Nu(1,1 + 3*i) = Shape(i);  Nu(1,2 + 3*i) = 0;    
        Nu(2,0 + 3*i) = 0;         Nu(2,1 + 3*i) = 0;         Nu(2,2 + 3*i) = Shape(i); 
    }
} //}}}

inline void FlowElem::Bp_Matrix(LinAlg::Matrix<REAL> const & derivs,  //{{{
                              LinAlg::Matrix<REAL> const & J, 
                              LinAlg::Matrix<REAL> & Bp) const
{
    // Resize Bp matrix
    Bp.Resize(NDIM, _n_nodes);

    // Inverse of Jacobian
    LinAlg::Matrix<REAL> inv_J(J.Rows(),J.Cols());
    J.Inv(inv_J);

    // Cartesian derivatives
    LinAlg::Matrix<REAL> cart_derivs;
    cart_derivs.Resize(J.Rows(), derivs.Cols());
    LinAlg::Gemm(1.0, inv_J, derivs, 0.0, cart_derivs); // cart_derivs <- inv_J*derivs

    // Loop along all nodes of the element
    for (int i=0; i<_n_nodes; ++i) //
    {
        Bp(0, i) = cart_derivs(0,i);     
        Bp(1, i) = cart_derivs(1,i);     
        Bp(2, i) = cart_derivs(2,i);     
    }
} //}}}

inline void FlowElem::NodalDOFs(int iNode, Array<FEM::Node::DOFVarsStruct*> & DOFs) const // {{{
{
    // Resize DOFs array:  1 pore-pressure
    DOFs.resize(1);

    // Set DOFs array with pointers to DOFVarsStruct inside this element Node [iNode]
    DOFs[0] = &_connects[iNode]->DOFVar(DWP);

} // }}}

inline void FlowElem::BackupState() // {{{
{
    for (int i=0; i<_n_int_pts; ++i)
        _a_model[i]->BackupState();
} // }}}

inline void FlowElem::UpdateState(REAL TimeInc) // Read nodal values _IncEssential and write to _IncNatural  {{{
{

    //   Internal volume integration:
    //   ============================
    //       
    //            int    /    T
    //         {dQ}   =  | [B]  * h * {v} * dV               
    //                   /    p                   
    //           
    //         h = TimeInc
    
    // Allocate (local/element) pore-pressures vectors
    LinAlg::Vector<REAL>  P(_n_nodes); 
    LinAlg::Vector<REAL> dP(_n_nodes); 
    
    // Loop along nodes
    for (int i_node=0; i_node<_n_nodes; ++i_node)
    {
        // Assemble (local/element) pore-pressures vectors. 
        dP(i_node) = _connects[i_node]->DOFVar(DWP)._IncEssenVal;
        P(i_node)  = _connects[i_node]->DOFVar(DWP).EssentialVal;
    }
    
    // Allocate (local/element) internal volume vector
    LinAlg::Vector<REAL> dQ(_n_nodes); // Delta internal volume of this element
    dQ.SetValues(0.0);
    
    // Loop along integration points
    for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
    {
        // Temporary Integration Points
        REAL r = _a_int_pts[i_ip].r;
        REAL s = _a_int_pts[i_ip].s;
        REAL t = _a_int_pts[i_ip].t;
        REAL w = _a_int_pts[i_ip].w;

        // Calculate Shape functions w.r.t local coordinate system
        LinAlg::Vector<REAL> shape; // size = NumLocalCoords(ex.: r,s,t) x _n_nodes
        Shape(r,s,t, shape);

        // Calculate Derivatives of Shape functions w.r.t local coordinate system
        LinAlg::Matrix<REAL> derivs; // size = NumLocalCoords(ex.: r,s,t) x _n_nodes
        Derivs(r,s,t, derivs);

        // Calculate J (Jacobian) matrix for i_ip Integration Point
        LinAlg::Matrix<REAL> J;
        Jacobian(derivs, J);
        REAL det_J = J.Det();

        //Calulate Bp matrix
        LinAlg::Matrix<REAL> Bp;
        Bp_Matrix(derivs, J, Bp);
        
        REAL gammaW = _a_model[i_ip]->GammaW();

        // Calc gradient vector
        Tensors::Tensor1 Grad;
        Grad=0,0,0;
        for (int i=0; i<Bp.Rows(); ++i)
            for (int j=0; j<Bp.Cols(); ++j)
                Grad(i) += Bp(i,j)*P(j)/gammaW; 

        // Update flow model
        REAL dp=0;                        // increment of pore-pressure
        for (int i=0; i<_n_nodes; ++i)    // Calculate the pore-pressure increment at current integration point
            dp += shape(i)*dP(i);
        _a_model[i_ip]->FlowUpdate(dp);   //Update the model
        
        //Calculate Velocity
        Tensors::Tensor1 Vel;
        _a_model[i_ip]->FlowVelocity(Grad,Vel);
        
        //Calculate internal volume vector
        for (int i=0; i<dQ.Size(); ++i)
            for (int j=0; j<NDIM; ++j)
                dQ(i) += Bp(j,i)*TimeInc*Vel(j)*det_J*w; 
    }

    // Update nodal _IncNaturVals
    for (int i_node=0; i_node<_n_nodes; ++i_node)
    {
        // Assemble (local/element) volumes vector. 
        _connects[i_node]->DOFVar(DWD)._IncNaturVal += dQ(i_node);
        // Update natural vars state
        _connects[i_node]->DOFVar(DWD).NaturalVal   += dQ(i_node);
    }

} // }}}

inline void FlowElem::RestoreState() // {{{
{
    for (int i=0; i<_n_int_pts; ++i)
        _a_model[i]->RestoreState();
} // }}}

inline String FlowElem::OutCenter(bool PrintCaptionOnly=false) const // {{{
{
    // Auxiliar variables
    std::ostringstream oss;

    // Number of state values
    int n_int_state_vals = _a_model[0]->nInternalStateValues();
    
    // Print caption
    if (PrintCaptionOnly)
    {
        // Stress and strains
        oss << _8s()<< DWP;

        // Internal state values
        Array<String> str_state_names;   _a_model[0]->InternalStateNames(str_state_names);
        for (int i=0; i<n_int_state_vals; ++i)
            oss << _8s()<< str_state_names[i];
        oss << std::endl;
    }
    else
    {
        // pore-pressures, volumes and internal state values evaluated at the center of the element
        REAL         p_cen = 0; //pore-pressure
        Array<REAL>  int_state_vals_cen;   int_state_vals_cen.assign(n_int_state_vals,0.0);
        
        // Loop over integration points
        for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
        {
            // pore-pressures
            p_cen += _a_model[i_ip]->Pp();

            // Internal state values
            Array<REAL> int_state_vals;    _a_model[i_ip]->InternalStateValues(int_state_vals);
            for (int j=0; j<n_int_state_vals; ++j)
                int_state_vals_cen[j] += int_state_vals[j];
        }
        
        // Average pore-pressures
        p_cen = p_cen/_n_int_pts; 
        
        // Average internal state values
        for (int j=0; j<n_int_state_vals; ++j)
            int_state_vals_cen[j] = int_state_vals_cen[j] / _n_int_pts;;

        // Output
        oss << _8s()<< p_cen;
        for (int j=0; j<n_int_state_vals; ++j)
            oss << _8s()<< int_state_vals_cen[j];
        oss << std::endl;
    }

    return oss.str(); 
} // }}}

void FlowElem::OutNodes(LinAlg::Matrix<REAL> & Values, Array<String> & Labels) const // {{{
{
    int const DATA_COMPS=2;
    Values.Resize(_n_nodes,DATA_COMPS);
    Labels.resize(DATA_COMPS);
    Labels[ 0] = DWP; Labels[ 1] = DWD; 
    for (int i_node=0; i_node<_n_nodes; i_node++)
    {
        Values(i_node,0) = _connects[i_node]->DOFVar(DWP).EssentialVal;
        Values(i_node,1) = _connects[i_node]->DOFVar(DWD).EssentialVal;
    }
} // }}}

inline void FlowElem::_calc_initial_internal_pore_pressures() // {{{
{
    // Here must be done an extrapolation form pore-pressures values into integration
    // points to nodal points
    // For simplicity, here is considered the mean value of integration points
    //
    // Here is considered an initial stationary state, it is the reason ...
    // ... to leave the calculus of initital internal volumes.
    
    // Allocate (local/element) internal force vector
    LinAlg::Vector<REAL> P(_n_nodes);
    P.SetValues(0.0);

    REAL mean_pp_value=0;
    // Loop along integration points
    for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
    {
        mean_pp_value+= _a_model[i_ip]->Pp()/_n_int_pts;
    }

    for (int i_node=0; i_node<_n_int_pts; ++i_node)
    {
        P(i_node) = mean_pp_value;
    }

    // Update nodal NaturVals
    for (int i_node=0; i_node<_n_nodes; ++i_node)
    {
        // Assemble (local/element) pore-pressures vector
        int node_refs = _connects[i_node]->Refs();
        _connects[i_node]->DOFVar(DWP).EssentialVal += P(i_node)/node_refs; // NaturalVal must be set to zero during AddDOF routine
    }
} // }}}

inline void FlowElem::RHSVector(size_t index, LinAlg::Vector<REAL> & V, Array<size_t> & Map) const // {{{
{
    assert(index == 0);
    V.Resize(_n_nodes);
    
    // Loop along nodes
    for (int i_node=0; i_node<_n_nodes; ++i_node)
    {
        // Assemble (local/element) total pore-presure vector. 
        V(i_node) = _connects[i_node]->DOFVar("p").EssentialVal;
    }
    //Mounting a map of positions f
    int idx=0; // position (idx) inside V vector
    Map.resize(V.Size()); 
    for (int i_node=0; i_node<_n_nodes; ++i_node)
    {
        Map[idx++] = _connects[i_node]->DOFVar("p").EqID; 
    }
} // }}}

inline void FlowElem::Order0Matrix(size_t index, LinAlg::Matrix<REAL> & M, Array<size_t> & RowsMap, Array<size_t> & ColsMap) const // {{{
{
    assert(index == 0);
    Permeability(M, RowsMap);
    ColsMap.resize(RowsMap.size());
    for (size_t i=0; i<RowsMap.size(); i++) ColsMap[i] = RowsMap[i];
} // }}}




// Register the DOF information of FlowElem into DOFInfoMap
int FlowDOFInfoRegister() // {{{
{
    // Temporary 
    DOFInfo D; 

    // Nodal
    D.NodeEssential.push_back(FlowElem::DWP  + _("@Nodal pore-presure increment"));
    D.NodeNatural  .push_back(FlowElem::DWD  + _("@Nodal water discharge increment"));

    // Face
    D.FaceEssential.push_back(FlowElem::DWP  + _("@Pore-presure increment per area"));
    D.FaceNatural  .push_back(FlowElem::DFWD + _("@Water discharge increment per area"));

    // Insert into DOFInfoMap
    DOFInfoMap["WaterFlow"] = D;

    return 0;
} // }}}

// Execute the autoregistration
int __FlowElemDOFInfo_dummy_int  = FlowDOFInfoRegister();

}; // namespace FEM

#endif // FEM_FLOWELEM_H

// vim:fdm=marker

---