equilibelem.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 MECHSYS_FEM_EQUILIB_H
#define MECHSYS_FEM_EQUILIB_H

#ifdef HAVE_CONFIG_H
  #include "config.h"
#else
  #ifndef REAL
    #define REAL double
  #endif
#endif

#include "util/string.h"
#include "models/equilibmodel.h"
#include "fem/ele/element.h"
#include "tensors/tensors.h"
#include "tensors/functions.h"
#include "util/numstreams.h"
#include "linalg/laexpr.h"

using Tensors::Stress_p_q_S_sin3th; // Mean and deviatoric Cambridge stress invariants
using Tensors::Strain_Ev_Ed;        // Volumetric and deviatoric strain invariants

namespace FEM
{

class EquilibElem : public virtual Element
{

public:
    // EquilibElem constants {{{    
    static int    NDIM;
    static int    NSTRESSCOMPS;
    static String DUX;
    static String DUY;
    static String DUZ;
    static String DFX;
    static String DFY;
    static String DFZ;
    static String DTX;
    static String DTY;
    static String DTZ;
    // }}}
    
    // Destructor
    virtual ~EquilibElem() {}

    // Virtual 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   Stiffness           (LinAlg::Matrix<REAL> & Ke, Array<size_t> & EqMap) const;
    virtual void   B_Matrix            (LinAlg::Matrix<REAL> const & derivs, LinAlg::Matrix<REAL> const & J, LinAlg::Matrix<REAL> & B) const;
    virtual void   NodalDOFs           (int iNode, Array<FEM::Node::DOFVarsStruct*> & DOFs) const;
    virtual void   SetProperties       (Array<REAL> const & EleProps) { _unit_weight=EleProps[0]; }
    virtual String OutCenter           (bool PrintCaptionOnly) const;
            void   OutNodes            (LinAlg::Matrix<REAL> & Values, Array<String> & Labels) const;
    virtual void   UpdateState         (REAL TimeInc);
    virtual void   CalcFaceNodalValues (String            const & FaceDOFName  ,         // In:  Face DOF Name, ex.: "Dtx, Dty, Dtz"
                                        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.: "Dfx, Dfy, Dfz"
                                        LinAlg::Vector<REAL>    & NodalValues  ) const;  // Out: Vector of nodal values
    virtual void Deactivate();

    // Methods
    void BackupState();
    void RestoreState();
    // Functions to assemble DAS matrices
    virtual size_t nOrder1Matrices() const { return 1; };
    virtual void   Order1Matrix(size_t index, LinAlg::Matrix<REAL> & M, Array<size_t> & RowsMap, Array<size_t> & ColsMap) const;
    // Output
    void OutTensor1 (String & Str) const; // Stress
    void OutTensor2 (String & Str) const; // Strain
    REAL OutScalar2 ()             const; // Ed

private:
    // Data
    Array<EquilibModel*> _a_model;
    REAL                 _unit_weight;

    // Methods
    void _set_node_vars                (int iNode);
    void _calc_initial_internal_forces ();

}; // class EquilibElem

// EquilibElem constants {{{    
int    EquilibElem::NDIM         = 3;
int    EquilibElem::NSTRESSCOMPS = 6;
String EquilibElem::DUX          = _T("Dux");
String EquilibElem::DUY          = _T("Duy");
String EquilibElem::DUZ          = _T("Duz");
String EquilibElem::DFX          = _T("Dfx");
String EquilibElem::DFY          = _T("Dfy");
String EquilibElem::DFZ          = _T("Dfz");
String EquilibElem::DTX          = _T("Dtx");
String EquilibElem::DTY          = _T("Dty");
String EquilibElem::DTZ          = _T("Dtz");
// }}}



    
inline void EquilibElem::_set_node_vars(int iNode) // {{{
{
    // Add Degree of Freedom to a node (Essential, Natural)
    _connects[iNode]->AddDOF(DUX, DFX);
    _connects[iNode]->AddDOF(DUY, DFY);
    _connects[iNode]->AddDOF(DUZ, DFZ);

    // Set vectors
    _connects[iNode]->SetEssentialVector(DUX, DUY, DUZ);
    _connects[iNode]->SetNaturalVector  (DFX, DFY, DFZ);
} // }}}

inline void EquilibElem::Deactivate() //{{{
{
    if (_is_active==false) return; // the element was already inactive
    _is_active = false; 
    //dereferencing nodes:
    for (int i_node=0; i_node<_n_nodes; i_node++)
    {
        _connects[i_node]->Refs()--;
    }
    //verifying if element is in the boundary of excavation
    bool in_boundary = false;
    for (int i_node=0; i_node<_n_nodes; i_node++)
        if (_connects[i_node]->Refs()>0) 
        {
            in_boundary=true;
            break;
        }
    //calculate internal forces for element on boudary
    if (in_boundary)
    {
        LinAlg::Vector<REAL> F(NDIM*_n_nodes);
        F.SetValues(0.0);
        
        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> B;
            B_Matrix(derivs,J, B);
                
            // Get tensor for accumulate stress
            Tensors::Tensor2 sig;
            sig = _a_model[i_ip]->Sig();

            // Calculate internal force vector;
            for (int i=0; i<F.Size(); ++i)
                for (int j=0; j<B.Rows(); ++j)
                    F(i) += B(j,i)*sig(j)*det_J*w; // dF = Sum(Bt*sig*det_J*w)
        }   

        // adding to boudary values of node if Refs()>0
        for (int i_node=0; i_node<_n_nodes; i_node++)
            if (_connects[i_node]->Refs()>0) 
            {
                _connects[i_node]->DOFVar(DFX).NaturalBry += F(i_node*NDIM  );
                _connects[i_node]->DOFVar(DFY).NaturalBry += F(i_node*NDIM+1);
                _connects[i_node]->DOFVar(DFZ).NaturalBry += F(i_node*NDIM+2);
            }
            else // eliminate dofs from nodes with no references
            {
                _connects[i_node]->ClearDOF();

            }
        //delete IP information 
    }
} //}}}

inline void EquilibElem::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(_("EquilibElem::ReAllocateModel: Array of array of initial data must have size==1 or size==%d"),_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] = AllocEquilibModel(ModelName, ModelPrms, AIniData[0]);
            else                    _a_model[i] = AllocEquilibModel(ModelName, ModelPrms, AIniData[i]);
        }

        // Calculate initial internal forces
        _calc_initial_internal_forces();
    }
    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] = AllocEquilibModel(ModelName, ModelPrms, AIniData[0]);
                else                    _a_model[i] = AllocEquilibModel(ModelName, ModelPrms, AIniData[i]);
            }
            // else: Do NOT re-allocate => Model not changed (don't do anything)
        }
    }
} // }}}

inline bool EquilibElem::IsEssential(String const & DOFName) const // {{{
{
    if (DOFName==DUX || DOFName==DUY || DOFName==DUZ) return true;
    else return false;
} // }}}

inline void EquilibElem::CalcFaceNodalValues(String            const & FaceDOFName  , // {{{
                                             REAL              const   FaceDOFValue ,
                                             Array<FEM::Node*> const & APtrFaceNodes,
                                             String                  & NodalDOFName ,
                                             LinAlg::Vector<REAL>    & NodalValues  ) const
{
    if (FaceDOFName==DTX || FaceDOFName==DTY || FaceDOFName==DTZ)
    {
        if (FaceDOFName==DTX) NodalDOFName=DFX;
        if (FaceDOFName==DTY) NodalDOFName=DFY;
        if (FaceDOFName==DTZ) NodalDOFName=DFZ;
        _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(_("EquilibElem::CalcFaceNodalValues: This method must only be called for FaceDOFName< %s > equal to Dtx, Dty or Dtz.\n %s"),
                FaceDOFName.c_str(), oss.str().c_str());
    }
} // }}}

inline void EquilibElem::Stiffness(LinAlg::Matrix<REAL> & Ke, Array<size_t> & EqMap) const // {{{
{
    // Resize Ke
    Ke.Resize(NDIM*_n_nodes, NDIM*_n_nodes); // sum(Bt*D*B*det(J)*w)
    Ke.SetValues(0.0);

    // Allocate entities used for every integration point
    LinAlg::Matrix<REAL> derivs; // size = NumLocalCoords(ex.: r,s,t) x _n_nodes
    LinAlg::Matrix<REAL> J;      // Jacobian matrix
    LinAlg::Matrix<REAL> B;      // strain-displacement matrix
    Tensors::Tensor4     tsrD;   // Fourth order constitutive tensor
    LinAlg::Matrix<REAL> D;      // Constitutive matrix

    // 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;

        Derivs(r,s,t, derivs);        // Calculate Derivatives of Shape functions w.r.t local coordinate system
        Jacobian(derivs, J);          // Calculate J (Jacobian) matrix for i_ip Integration Point
        B_Matrix(derivs,J, B);        // Calculate B matrix for i_ip Integration Point

        // Constitutive tensor 
        _a_model[i_ip]->TgStiffness(tsrD); Copy(tsrD, D);

        // Calculate Tangent Stiffness
        Ke += trn(B)*D*B*det(J)*w;
    }
    
    //Mounting a map of positions from Ke to Global
    int idx_Ke=0; // position (idx) inside Ke matrix
    EqMap.resize(Ke.Rows()); // size=Ke.Rows()=Ke.Cols()
    for (int i_node=0; i_node<_n_nodes; ++i_node)
    {
        // Fill map of Ke position to K position of DOFs components
        EqMap[idx_Ke++] = _connects[i_node]->DOFVar(DUX).EqID; 
        EqMap[idx_Ke++] = _connects[i_node]->DOFVar(DUY).EqID; 
        EqMap[idx_Ke++] = _connects[i_node]->DOFVar(DUZ).EqID; 
    }
} // }}}

inline void EquilibElem::NodalDOFs(int iNode, Array<FEM::Node::DOFVarsStruct*> & DOFs) 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

    // Resize DOFs array with 3 = 3D
    DOFs.resize(3);

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

} // }}}

inline void EquilibElem::B_Matrix(LinAlg::Matrix<REAL> const & derivs, LinAlg::Matrix<REAL> const & J, LinAlg::Matrix<REAL> & B) const // {{{
{
    /* OBS.:
     *       1) This B matrix considers Mandel representation of stress components,
     *          for this reason, some of its components are divided by sqrt(2).
     *          Ex.: Sig = {S11  S22  S33  sqrt(2)*S12  sqrt(2)*S23  sqrt(2)*S13}^(T)
     *               Eps = {E11  E22  E33  sqrt(2)*E12  sqrt(2)*E23  sqrt(2)*E13}^(T)
     *       2) This B matrix considers Soil Mechanics sign convention of stress and strains
     *          Ex.: Compressive stresses/strains are positive
     */
    
    // Resize B matrix
    B.Resize(NSTRESSCOMPS,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
    REAL sq2 = sqrt(2.0);
    REAL dNdX,dNdY,dNdZ;
    int  j=0; // j column of B
    for (int i=0; i<_n_nodes; ++i) // i row of B
    {
        // Assemble B matrix
        j=i*NDIM;
        dNdX=-cart_derivs(0,i);  dNdY=-cart_derivs(1,i);  dNdZ=-cart_derivs(2,i); // Negative values => Soil mechanics convention
        B(0,0+j) =     dNdX;     B(0,1+j) =      0.0;     B(0,2+j) =      0.0;
        B(1,0+j) =      0.0;     B(1,1+j) =     dNdY;     B(1,2+j) =      0.0;
        B(2,0+j) =      0.0;     B(2,1+j) =      0.0;     B(2,2+j) =     dNdZ;
        B(3,0+j) = dNdY/sq2;     B(3,1+j) = dNdX/sq2;     B(3,2+j) =      0.0; // sq2 => Mandel representation
        B(4,0+j) =      0.0;     B(4,1+j) = dNdZ/sq2;     B(4,2+j) = dNdY/sq2; // sq2 => Mandel representation
        B(5,0+j) = dNdZ/sq2;     B(5,1+j) =      0.0;     B(5,2+j) = dNdX/sq2; // sq2 => Mandel representation
    }
} // }}}

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

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

    // Allocate (local/element) displacements vector
    LinAlg::Vector<REAL> dU(NDIM*_n_nodes); // Delta disp. of this element
    
    // Loop along nodes
    for (int i_node=0; i_node<_n_nodes; ++i_node)
    {
        // Assemble (local/element) displacements vector. OBS.: NDIM=3
        dU(i_node*NDIM  ) = _connects[i_node]->DOFVar(DUX)._IncEssenVal;
        dU(i_node*NDIM+1) = _connects[i_node]->DOFVar(DUY)._IncEssenVal;
        dU(i_node*NDIM+2) = _connects[i_node]->DOFVar(DUZ)._IncEssenVal;
    }
    
    // Allocate (local/element) internal force vector
    LinAlg::Vector<REAL> dF(NDIM*_n_nodes); // Delta internal force of this element
    dF.SetValues(0.0);
    
    // Allocate entities used for every integration point
    LinAlg::Matrix<REAL> derivs;  // size = NumLocalCoords(ex.: r,s,t) x _n_nodes
    LinAlg::Matrix<REAL> J;       // Jacobian matrix
    LinAlg::Matrix<REAL> B;       // strain-displacement matrix
    Tensors::Tensor2     tsrDEps; // Strain tensor  tsrDEps(NSTRESSCOMPS=6) = B(NSTRESSCOMPS=6,NDIM*_n_nodes) * dU(NDIM*_n_nodes)
    Tensors::Tensor2     tsrDSig; // Stress tensor
    LinAlg::Vector<REAL> DEps;    // Strain vector in Mandel's notation 
    LinAlg::Vector<REAL> DSig;    // Stress vector in Mandel's notation 

    // 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;

        Derivs(r,s,t, derivs);        // Calculate Derivatives of Shape functions w.r.t local coordinate system
        Jacobian(derivs, J);          // Calculate J (Jacobian) matrix for i_ip Integration Point
        B_Matrix(derivs, J, B);       // Calculate B matrix for i_ip Integration Point

        // Calculate a tensor for the increments of strain
        DEps = B*dU;
        Copy(DEps, tsrDEps);  // tsrDEps = B*dU
        
        // Update model
        _a_model[i_ip]->StressUpdate(tsrDEps, tsrDSig);
        Copy(tsrDSig, DSig);

        // Calculate internal force vector;
        dF += trn(B)*DSig*det(J)*w;
    }

    // Update nodal _IncNaturVals
    for (int i_node=0; i_node<_n_nodes; ++i_node)
    {
        // Assemble (local/element) displacements vector. OBS.: NDIM=3
        _connects[i_node]->DOFVar(DFX)._IncNaturVal += dF(i_node*NDIM  );
        _connects[i_node]->DOFVar(DFY)._IncNaturVal += dF(i_node*NDIM+1);
        _connects[i_node]->DOFVar(DFZ)._IncNaturVal += dF(i_node*NDIM+2);

        // Update natural vars state
        _connects[i_node]->DOFVar(DFX).NaturalVal += dF(i_node*NDIM  );
        _connects[i_node]->DOFVar(DFY).NaturalVal += dF(i_node*NDIM+1);
        _connects[i_node]->DOFVar(DFZ).NaturalVal += dF(i_node*NDIM+2);
    }
} // }}}

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

inline String EquilibElem::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()<< "p"  << _8s()<< "q"  << _8s()<< "sin3th" << _8s()<< "Ev"  << _8s()<< "Ed";
        oss << _8s()<< "Sx" << _8s()<< "Sy" << _8s()<< "Sz" << _8s()<< "Sxy" << _8s()<< "Syz" << _8s()<< "Sxz";
        oss << _8s()<< "Ex" << _8s()<< "Ey" << _8s()<< "Ez" << _8s()<< "Exy" << _8s()<< "Eyz" << _8s()<< "Exz";

        // 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
    {
        // Stress, strains and internal state values evaluated at the center of the element
        Tensors::Tensor2 sig_cen(0.0);
        Tensors::Tensor2 eps_cen(0.0);
        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)
        {
            // Stress and strains
            sig_cen += _a_model[i_ip]->Sig();
            eps_cen += _a_model[i_ip]->Eps();

            // 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 stress and strains
        sig_cen = sig_cen / _n_int_pts;
        eps_cen = eps_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;;

        // Calculate stress invariants
        REAL              p,q,sin3th;
        Tensors::Tensor2  S;
        Stress_p_q_S_sin3th(sig_cen,p,q,S,sin3th);

        // Calculate strain invariants
        REAL     Ev,Ed;
        Strain_Ev_Ed(eps_cen,Ev,Ed);

        // Output
        REAL sq2 = sqrt(2.0);
        oss << _8s()<< p          << _8s()<< q          << _8s()<< sin3th << _8s()<< Ev*100.0       << _8s()<< Ed*100.0;
        oss << _8s()<< sig_cen(0) << _8s()<< sig_cen(1) << _8s()<< sig_cen(2)         << _8s()<< sig_cen(3)/sq2 << _8s()<< sig_cen(4)/sq2 << _8s()<< sig_cen(5)/sq2;
        oss << _8s()<< eps_cen(0) << _8s()<< eps_cen(1) << _8s()<< eps_cen(2)         << _8s()<< eps_cen(3)/sq2 << _8s()<< eps_cen(4)/sq2 << _8s()<< eps_cen(5)/sq2;
        for (int j=0; j<n_int_state_vals; ++j)
            oss << _8s()<< int_state_vals_cen[j];
        oss << std::endl;
    }

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

inline void EquilibElem::OutNodes(LinAlg::Matrix<REAL> & Values, Array<String> & Labels) const // {{{
{
    int const DATA_COMPS=18;
    Values.Resize(_n_nodes,DATA_COMPS);
    Labels.resize(DATA_COMPS);
    Labels[ 0] = DUX ; Labels[ 1] = DUY ; Labels[ 2] = DUZ; 
    Labels[ 3] = DFX ; Labels[ 4] = DFY ; Labels[ 5] = DFZ;
    Labels[ 6] = "Ex"; Labels[ 7] = "Ey"; Labels[ 8] = "Ez"; Labels[ 9] = "Exy"; Labels[10] = "Eyz"; Labels[11] = "Exz";
    Labels[12] = "Sx"; Labels[13] = "Sy"; Labels[14] = "Sz"; Labels[15] = "Sxy"; Labels[16] = "Syz"; Labels[17] = "Sxz";
    for (int i_node=0; i_node<_n_nodes; i_node++)
    {
        Values(i_node,0) = _connects[i_node]->DOFVar(DUX).EssentialVal;
        Values(i_node,1) = _connects[i_node]->DOFVar(DUY).EssentialVal;
        Values(i_node,2) = _connects[i_node]->DOFVar(DUZ).EssentialVal;
        Values(i_node,3) = _connects[i_node]->DOFVar(DFX).NaturalVal;
        Values(i_node,4) = _connects[i_node]->DOFVar(DFY).NaturalVal;
        Values(i_node,5) = _connects[i_node]->DOFVar(DFZ).NaturalVal;
    }

    //Extrapolation
    LinAlg::Vector<REAL> ip_values(_n_int_pts);
    LinAlg::Vector<REAL> nodal_values(_n_nodes);
    
    // Strains
    for (int i_comp=0; i_comp<NSTRESSCOMPS; i_comp++)
    {
        for (int j_ip=0; j_ip<_n_int_pts; j_ip++)
        {
            Tensors::Tensor2 const & eps = _a_model[j_ip]->Eps();
            ip_values(j_ip) = eps(i_comp); //getting IP values
        }
        _extrapolate(ip_values, nodal_values);
        for (int j_node=0; j_node<_n_nodes; j_node++)
            Values(j_node,i_comp+6) = nodal_values(j_node);
    }
    
    // Stresses
    for (int i_comp=0; i_comp<NSTRESSCOMPS; i_comp++)
    {
        for (int j_ip=0; j_ip<_n_int_pts; j_ip++)
        {
            Tensors::Tensor2 const & sig = _a_model[j_ip]->Sig();
            ip_values(j_ip) = sig(i_comp); //getting IP values
        }
        _extrapolate(ip_values, nodal_values);
        for (int j_node=0; j_node<_n_nodes; j_node++)
            Values(j_node,i_comp+12) = nodal_values(j_node);
    }
} // }}}

inline void EquilibElem::_calc_initial_internal_forces() // {{{
{
    // Allocate (local/element) internal force vector
    LinAlg::Vector<REAL> F(NDIM*_n_nodes);
    F.SetValues(0.0);

    // Allocate entities used for every integration point
    LinAlg::Matrix<REAL> derivs;  // size = NumLocalCoords(ex.: r,s,t) x _n_nodes
    LinAlg::Matrix<REAL> J;       // Jacobian matrix
    LinAlg::Matrix<REAL> B;       // strain-displacement matrix
    LinAlg::Vector<REAL> Sig;     // Stress vector in Mandel's notation 

    // 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;

        Derivs(r,s,t, derivs);        // Calculate Derivatives of Shape functions w.r.t local coordinate system
        Jacobian(derivs, J);          // Calculate J (Jacobian) matrix for i_ip Integration Point
        B_Matrix(derivs, J, B);       // Calculate B matrix for i_ip Integration Point

        Tensors::Tensor2 const & tsrSig = _a_model[i_ip]->Sig(); // Declare a reference to the actual (initial) stress tensor inside EquilibModel
        Copy(tsrSig, Sig);

        // Calculate internal force vector;
        F += trn(B)*Sig*det(J)*w;
    }

    // Update nodal NaturVals
    for (int i_node=0; i_node<_n_nodes; ++i_node)
    {
        // Assemble (local/element) displacements vector. OBS.: NDIM=3
        _connects[i_node]->DOFVar(DFX).NaturalVal += F(i_node*NDIM  ); // NaturalVal must be set to zero during AddDOF routine
        _connects[i_node]->DOFVar(DFY).NaturalVal += F(i_node*NDIM+1);
        _connects[i_node]->DOFVar(DFZ).NaturalVal += F(i_node*NDIM+2);
    }
} // }}}

inline void EquilibElem::Order1Matrix(size_t index, LinAlg::Matrix<REAL> & M, Array<size_t> & RowsMap, Array<size_t> & ColsMap) const // {{{
{
    assert(index == 0);
    Stiffness(M, RowsMap);
    ColsMap.resize(RowsMap.size());
    for (size_t i=0; i<RowsMap.size(); i++) ColsMap[i] = RowsMap[i];
} // }}}
    
inline void EquilibElem::OutTensor1(String & Str) const // {{{
{
    // Stress evaluated at the center of the element
    Tensors::Tensor2 s(0.0);

    // Loop over integration points
    for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
        s += _a_model[i_ip]->Sig();
    
    // Average stress
    s = s / _n_int_pts;

    // Output
    REAL sq2 = sqrt(2.0);
    Str.Printf(_(" %e %e %e  %e %e %e  %e %e %e "), s(0),s(3)/sq2,s(5)/sq2,  s(3)/sq2,s(1),s(4)/sq2,  s(5)/sq2,s(4)/sq2,s(2));

} // }}}

inline void EquilibElem::OutTensor2(String & Str) const // {{{
{
    // Strains evaluated at the center of the element
    Tensors::Tensor2 e(0.0);

    // Loop over integration points
    for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
        e += 100.0*_a_model[i_ip]->Eps();
    
    // Average strains
    e = e / _n_int_pts;

    // Output
    REAL sq2 = sqrt(2.0);
    Str.Printf(_(" %e %e %e  %e %e %e  %e %e %e "), e(0),e(3)/sq2,e(5)/sq2,  e(3)/sq2,e(1),e(4)/sq2,  e(5)/sq2,e(4)/sq2,e(2));

} // }}}

inline REAL EquilibElem::OutScalar2() const // {{{
{
    // Strains evaluated at the center of the element
    Tensors::Tensor2 e(0.0);

    // Loop over integration points
    for (int i_ip=0; i_ip<_n_int_pts; ++i_ip)
        e += 100.0*_a_model[i_ip]->Eps();
    
    // Average strains
    e = e / _n_int_pts;

    // Calculate strain invariants
    REAL   Ev,Ed;
    Strain_Ev_Ed(e,Ev,Ed);

    // Output
    return Ed;

} // }}}




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

    // Nodal
    D.NodeEssential.push_back(EquilibElem::DUX + _("@Nodal displacement increment in x direction"));
    D.NodeEssential.push_back(EquilibElem::DUY + _("@Nodal displacement increment in y direction"));
    D.NodeEssential.push_back(EquilibElem::DUZ + _("@Nodal displacement increment in z direction"));
    D.NodeNatural  .push_back(EquilibElem::DFX + _("@Nodal force increment in x direction"));
    D.NodeNatural  .push_back(EquilibElem::DFY + _("@Nodal force increment in y direction"));
    D.NodeNatural  .push_back(EquilibElem::DFZ + _("@Nodal force increment in z direction"));

    // Face
    D.FaceEssential.push_back(EquilibElem::DUX + _("@Displacement increment in x direction on face"));
    D.FaceEssential.push_back(EquilibElem::DUY + _("@Displacement increment in y direction on face"));
    D.FaceEssential.push_back(EquilibElem::DUZ + _("@Displacement increment in z direction on face"));
    D.FaceNatural  .push_back(EquilibElem::DTX + _("@Traction increment in x direction on face"));
    D.FaceNatural  .push_back(EquilibElem::DTY + _("@Traction increment in y direction on face"));
    D.FaceNatural  .push_back(EquilibElem::DTZ + _("@Traction increment in z direction on face"));

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

    return 0;
} // }}}

// Execute the autoregistration
int __EquilibElemDOFInfo_dummy_int  = EquilibDOFInfoRegister();

}; // namespace FEM

#endif // MECHSYS_FEM_EQUILIB_H

// vim:fdm=marker

---