element.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_ELEMENT_H
#define MECHSYS_FEM_ELEMENT_H

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

#include <map>

#include "util/string.h"
#include "fem/node.h"
#include "linalg/vector.h"
#include "linalg/matrix.h"
#include "linalg/lawrap.h"
#include "linalg/laexpr.h"

namespace FEM
{

class Element
{
public:
    // Destructor
    virtual ~Element() {}
    // Auxiliar structure
    struct IntegPoint
    {
        REAL r; 
        REAL s; 
        REAL t; 
        REAL w; 
    };
    int Number; //auxiliar variable
    virtual String Name() const=0;
    // Local methods
    void SetNode(int iNode, FEM::Node * ptrNode); // Copy a pointer of node iNode to the internal connects array
    // Virtual methods
    virtual bool IsEssential(String const & DOFName) const =0;
    virtual void 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    // Out: Vector of nodal values
    ) const =0;
    virtual void ReAllocateModel(String      const & ModelName, Array<REAL> const & ModelPrms, Array<Array<REAL> > const & AIniData) =0;
    void Activate()   { _is_active = true; }
    virtual void Deactivate() { _is_active = false; }
    virtual int  VTKCellType() const =0;

    bool     IsActive() const { return _is_active; }
    int        nNodes() const { return _n_nodes;   }
    size_t nIntPoints() const { return _n_int_pts; }

    // 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 =0;
    Array<Node*> const & Connects() const; // Return ana array containing pointers to nodes
            void IntegPoints(IntegPoint const * & IPs) const { IPs = _a_int_pts; }; // Return a pointer to the array of integration points
    virtual void BackupState() =0;
    virtual void UpdateState(REAL TimeInc=0) =0; // Read nodal values _IncEssential and write to _IncNatural
    virtual void RestoreState() =0;

    virtual void SetProperties(Array<REAL> const & EleProps) =0;

    virtual String OutCenter(bool PrintCaptionOnly=false) const =0;
    virtual void    OutNodes(LinAlg::Matrix<REAL> & Values, Array<String> & Labels) const { };
    virtual void       Shape(REAL r, REAL s, REAL t, LinAlg::Vector<REAL> & Shape) const =0;
    virtual void      Derivs(REAL r, REAL s, REAL t, LinAlg::Matrix<REAL> & Derivs) const =0;
    virtual void    Jacobian(REAL r, REAL s, REAL t, LinAlg::Matrix<REAL> & J) const;
    virtual void    Jacobian(LinAlg::Matrix<REAL> const & derivs, LinAlg::Matrix<REAL> & J) const;
    virtual void      Coords(LinAlg::Matrix<REAL> & coords) const;
    virtual void   InverseMap(REAL x, REAL y, REAL z, REAL & r, REAL & s, REAL & t) const;
            bool   IsInside(REAL x, REAL y, REAL z) const;
    virtual REAL   BoundDistance(REAL r, REAL s, REAL t) const { return -1; };

    // Functions to assemble DAS matrices
    virtual size_t nOrder0Matrices() const { return 0; };
    virtual size_t nOrder1Matrices() const { return 0; };
    virtual void      RHSVector(size_t index, LinAlg::Vector<REAL> & V, Array<size_t> &     Map) const { };
    virtual void   Order0Matrix(size_t index, LinAlg::Matrix<REAL> & M, Array<size_t> & RowsMap, Array<size_t> & ColsMap) const { };
    virtual void   Order1Matrix(size_t index, LinAlg::Matrix<REAL> & M, Array<size_t> & RowsMap, Array<size_t> & ColsMap) const { };
    
    // Output Scalars and Tensors
    virtual REAL OutScalar1 ()             const { return 0;                               }; // to be overriden.  for example, Pore-pressure
    virtual REAL OutScalar2 ()             const { return 0;                               }; // to be overriden.  for example, Ed
    virtual void OutTensor1 (String & Str) const { Str.Printf(_(" 0 0 0  0 0 0  0 0 0 ")); }; // to be overriden.  for example, Stress
    virtual void OutTensor2 (String & Str) const { Str.Printf(_(" 0 0 0  0 0 0  0 0 0 ")); }; // to be overriden.  for example, Strain
    virtual void OutTensor3 (String & Str) const { Str.Printf(_(" 0 0 0  0 0 0  0 0 0 ")); }; // to be overriden.  for example, Plastic Strain
protected:
    // Data
    int           _n_nodes;
    int           _n_int_pts;
    int           _n_face_nodes;   // Number of nodes of a face
    int           _n_face_int_pts; // Number of integ points of a a face
    Array<Node*>  _connects;       // Connectivity (ptr to nodes of this element). size=_n_nodes
    bool          _is_active;      // Flag for active/inactive condition
    IntegPoint const * _a_int_pts; // Array of Integration Points
    // Methods
    virtual void _set_node_vars(int iNode)=0; // Setup nodal variables (dx,dy, etc.) for node [iNode]
    virtual void _face_shape(REAL r, REAL s, LinAlg::Vector<REAL> & AFaceShape) const =0;
    void         _face_jacobian(Array<FEM::Node*> const & FaceConnects, REAL const r, REAL const s, LinAlg::Matrix<REAL> & J) const;
    virtual void _distrib_val_to_face_nodal_vals(Array<Node*>         const & FaceConnects, // In:  Array of ptrs to face nodes
                                                  REAL                const   FaceValue   , // In:  A value applied on a face to be converted to nodes
                                                  LinAlg::Vector<REAL>      & NodalValues   // Out: The resultant nodal values
    ) const =0;
    virtual void _face_derivs(REAL r, REAL s,         LinAlg::Matrix<REAL> & AFaceDerivs) const =0;
    virtual void _extrapolate(LinAlg::Vector<REAL> & IPValues, LinAlg::Vector<REAL> & NodalValues) const;
}; // class Element




inline void Element::SetNode(int iNode, FEM::Node * ptrNode) // {{{
{
    // Check
    if (_connects.size()<=0)                       throw new Fatal(_("Element::SetNode: _connects.size (=%d) must be greater than zero iNode=%d"),iNode,_connects.size());
    if (iNode<0)                                   throw new Fatal(_("Element::SetNode: iNode (=%d) must be greater than or equal to zero"),iNode);
    if (iNode>=static_cast<int>(_connects.size())) throw new Fatal(_("Element::SetNode: iNode (=%d) must be smaller than _connects.size (=%d)"),iNode,_connects.size());
    if (ptrNode==NULL)                             throw new Fatal(_("Element::SetNode: ptrNode is NULL! (iNode=%d)"),iNode);

    // Connects
    _connects[iNode] = ptrNode;

    // Nodal vars
    _set_node_vars(iNode);
}

// }}}

inline Array<Node*> const & Element::Connects() const // {{{
{
    return _connects;
} // }}}

void Element::InverseMap(REAL x, REAL y, REAL z, REAL & r, REAL & s, REAL & t) const //{{{
{
    LinAlg::Vector<REAL> shape;
    LinAlg::Matrix<REAL> derivs;
    LinAlg::Matrix<REAL> J; //Jacobian matrix
    //LinAlg::Matrix<REAL> inv_jacobian;
    //LinAlg::Matrix<REAL> trn_inv_jacobian;
    LinAlg::Vector<REAL> f(3);
    LinAlg::Vector<REAL> delta(3);
    REAL tx, ty, tz; //x, y, z trial
    REAL norm_f;
    r = s = t =0; // first suposition for natural coordinates
    int k=0;
    do
    {
        k++;
        Shape (r, s, t, shape);
        Derivs(r, s, t, derivs);
        Jacobian(derivs, J);
        tx = ty = tz = 0; 
        //calculate trial of real coordinates
        for (int j=0; j<_n_nodes; j++) 
        {
            tx += shape(j)*_connects[j]->X(); //ok
            ty += shape(j)*_connects[j]->Y(); //ok
            tz += shape(j)*_connects[j]->Z(); //ok
        }
        
        // calculate the error
        f(0) = tx - x;
        f(1) = ty - y;
        f(2) = tz - z;
        
        //jacobian.Inv(inv_jacobian);
        //inv_jacobian.Trn(trn_inv_jacobian);
        delta = trn(inv(J))*f;
        //Gemv(1.0, trn_inv_jacobian, f, 0.0, delta);
        
        r -= delta(0);
        s -= delta(1);
        t -= delta(2);

        //norm_f = sqrt(f(0)*f(0)+f(1)*f(1)+f(2)*f(2));
        norm_f = sqrt(trn(f)*f);
        if (k>25) break;
    } while(norm_f>1e-4);
} 

bool Element::IsInside(REAL x, REAL y, REAL z) const //{{{
{
    REAL tiny = 1e-4;
    REAL huge = 1e+20;
    REAL max, min;
    
    //fast search in X -----------------------------------------------------------------------
    max = -huge;
    for (int i=0; i < nNodes(); i++) if (_connects[i]->X() > max) max=_connects[i]->X();
    if ( x > max ) return false;
    min = +huge;
    for (int i=0; i < nNodes(); i++) if (_connects[i]->X() < min) min=_connects[i]->X();
    if ( x < min ) return false;

    //fast search in Y -----------------------------------------------------------------------
    max = -huge;
    for (int i=0; i < nNodes(); i++) if (_connects[i]->Y() > max) max=_connects[i]->Y();
    if ( y > max ) return false;
    min = +huge;
    for (int i=0; i < nNodes(); i++) if (_connects[i]->Y() < min) min=_connects[i]->Y();
    if ( y < min ) return false;
    
    //fast search in Z -----------------------------------------------------------------------
    max = -huge;
    for (int i=0; i < nNodes(); i++) if (_connects[i]->Z() > max) max=_connects[i]->Z();
    if ( z > max ) return false;
    min = +huge;
    for (int i=0; i < nNodes(); i++) if (_connects[i]->Z() < min) min=_connects[i]->Z();
    if ( z < min ) return false;
    
    REAL r, s, t;
    InverseMap(x,y,z,r,s,t);
    if (BoundDistance(r,s,t)>-tiny) return true;
    else return false;
} // }}}

inline void Element::_face_jacobian(Array<FEM::Node*> const & FaceConnects, REAL const r, REAL const s, LinAlg::Matrix<REAL> & J) const // {{{
{
    // Calculate the shape function derivatives for a face
    LinAlg::Matrix<REAL> m_face_derivs;
    _face_derivs(r, s, m_face_derivs);

    // Get the face coordinates
    LinAlg::Matrix<REAL> m_face_coords(_n_face_nodes,3);
    for (int i=0; i<_n_face_nodes; i++)
    {
        m_face_coords(i,0) = FaceConnects[i]->X();
        m_face_coords(i,1) = FaceConnects[i]->Y();
        m_face_coords(i,2) = FaceConnects[i]->Z();
    }

    // Determine face jacobian
    J.Resize(2,3);
    //LinAlg::Gemm(1, m_face_derivs, m_face_coords, 0, J);
    J = m_face_derivs*m_face_coords;
} // }}}

inline void Element::Coords(LinAlg::Matrix<REAL> & coords) const // {{{
{
    // Calculate a matrix with nodal coordinates
    coords.Resize(_n_nodes,3);   // 3 => X,Y,Z (all nodes have X,Y and Z)

    // Loop along all nodes of this element
    for (int i=0; i<_n_nodes; i++)
    {
        coords(i,0) = _connects[i]->X();
        coords(i,1) = _connects[i]->Y();
        coords(i,2) = _connects[i]->Z();
    }
} // }}}

inline void Element::Jacobian(LinAlg::Matrix<REAL> const & derivs, LinAlg::Matrix<REAL> & J) const // {{{
{
    // Calculate a matrix with nodal coordinates
    LinAlg::Matrix<REAL> cmatrix; // size = _n_nodes x 3
    cmatrix.Resize(_n_nodes,3);   // 3 => X,Y,Z (all nodes have X,Y and Z)

    // Loop along all nodes of this element
    for (int i=0; i<_n_nodes; i++)
    {
        cmatrix(i,0) = _connects[i]->X();
        cmatrix(i,1) = _connects[i]->Y();
        cmatrix(i,2) = _connects[i]->Z();
    }

    // Calculate the Jacobian; Size = NumLocalCoords x 3
    J = derivs * cmatrix;
} // }}}

inline void Element::Jacobian(REAL r, REAL s, REAL t, LinAlg::Matrix<REAL> & J) const // {{{
{
    LinAlg::Matrix<REAL> derivs; Derivs(r, s, t, derivs);
    LinAlg::Matrix<REAL> coords; Coords(coords);
    J = derivs*coords;
} // }}}

void Element::_extrapolate(LinAlg::Vector<REAL> & IPValues, LinAlg::Vector<REAL> & NodalValues) const  // {{{
{
    //Extrapolation:
    //                  IPValues = E * NodalValues;
    //  where:
    //                              t           t
    //                         E = N * inv(N * N )
    //
    //  and            N = [shape functions matrix]
    //                        1     2       ...     nNodes
    //                 1    [[N_11 N_12
    //                 2     [N_21
    //                 :     [
    //                nIP    [N_ ...                    ]]

    LinAlg::Matrix<REAL> N(_n_int_pts, _n_nodes);  // matrix of all IP shape functions
    LinAlg::Matrix<REAL> E(_n_nodes, _n_int_pts);  // Extrapolator matrix
    LinAlg::Vector<REAL> shape(_n_nodes);
    
    //filling N matrix
    for (int i_ip=0; i_ip<_n_int_pts; i_ip++)
    {
        REAL r = _a_int_pts[i_ip].r;
        REAL s = _a_int_pts[i_ip].s;
        REAL t = _a_int_pts[i_ip].t;
        Shape(r, s, t, shape);
        for (int j_node=0; j_node<_n_nodes; j_node++)
            N(i_ip, j_node) = shape(j_node);
    }

    //calculate extrapolator matrix
    assert(_n_nodes>=_n_int_pts); // TODO: replace with trow
    E = trn(N)*inv(N*trn(N));

    //perform extrapolation
    NodalValues = E*IPValues;
    
}; 



// Define a pointer to a function that makes (allocate) a new element
typedef Element * (*ElementMakerPtr)();

// Typdef of the array map that contains all the pointers to the functions that makes elements
typedef std::map<String, ElementMakerPtr, std::less<String> > ElementFactory_t;

// Instantiate the array map that contains all the pointers to the functions that makes elements
ElementFactory_t ElementFactory;

// Allocate a new element according to a string giving the name of the element
Element * AllocElement(String const & ElementName) // {{{
{
    ElementMakerPtr ptr=NULL;
    ptr = ElementFactory[ElementName];
    if (ptr==NULL)
        throw new Fatal(_("FEM::AllocElement: There is no < %s > implemented in this library"), ElementName.c_str());
    return (*ptr)();
} // }}}

}; // namespace FEM




// Define a structure to contain dof information from derivate elements
struct DOFInfo
{
    Array<String> NodeEssential;
    Array<String> NodeNatural;
    Array<String> FaceEssential;
    Array<String> FaceNatural;
};

// Typdef of the map that contains all dof info from problem type elements
typedef std::map<String, DOFInfo, std::less<String> > DOFInfoMap_t;

// Instantiate the map that contains all dof info from problem type elements
DOFInfoMap_t DOFInfoMap;

#endif // MECHSYS_FEM_ELEMENT

// vim:fdm=marker

---