stagesmanager.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_STAGESMANAGER_H
#define MECHSYS_FEM_STAGESMANAGER_H

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

#include <algorithm>

#include "util/string.h"
#include "fem/inputdata.h"
#include "fem/filesdata.h"
#include "fem/data.h"

namespace FEM
{


class StagesManager
{
public:
    // Friend prototype for operator <<
    friend std::ostream & operator<< (std::ostream & os, FEM::StagesManager const & SM);
    // Methods
    void Initialize      (InputData * ID, FilesData * FDat, Data * data); // FilesData could be altered in analyis using embedded elements
    int  nStages         () const       { return _idat->nSTAGES;      }
    int  CurrentStageNum () const       { return _current_stage_num;  }
    void ActivateStage   (int iStage);
    void DumpStage       (int nStage) {}
    void RestoreStage    (int nStage) {}
private:
    // Data
    int               _current_stage_num;
    InputData       * _idat;
    FilesData       * _fdat;
    Data            * _data;
    // Functions
    void _alloc_and_fill_nodal_coords();
    void _alloc_and_init_elements(); // Allocate elements, initialize connectivities and setup nodal variables (dx,dy,fx,fz, ...)
    void _alloc_and_init_embededs();
    void _clear_bry_conditions(); 
    void _update_nodes_given_face_bry(int iStage); // 1st) Read .face bry information and set Nodes bry values (convert acording to element type)
    void _update_nodes_given_node_bry(int iStage); // 2nd) Read .node bry information and set Nodes bry values
    void _set_element_attributes(int iStage);
}; // class StagesManager




inline void StagesManager::Initialize(InputData * ID, FilesData * FDat, Data * data) // {{{
{
    _current_stage_num = -1;
                 _idat = ID;
                 _fdat = FDat;
                 _data = data;
    _alloc_and_fill_nodal_coords(); // Allocate and fill Nodes inside _data;
    _alloc_and_init_elements();     // Allocate elements array with they respective type, ex.: Hex8, Bar2, etc.
} // }}}

inline void StagesManager::_alloc_and_fill_nodal_coords() // {{{
{
    // Allocate _data->Nodes array with information from FilesData _fdat
    _data->Nodes.resize(_fdat->Nodes.size());

    // Loop along the number of nodes
    for (size_t i=0; i<_data->Nodes.size(); ++i)
    {
        // Copy coordinates from FilesData _fdat to _data structure
        _data->Nodes[i].Initialize(               i,
                                  _fdat->Nodes[i].X,
                                  _fdat->Nodes[i].Y,
                                  _fdat->Nodes[i].Z);
    }
} // }}}

inline void StagesManager::_alloc_and_init_elements() // {{{ 
{
    // Allocate _data->Elements array with information from FilesData _fdat
    _data->Elements.resize(_fdat->Elements.size());

    // Declare a reference to the first stage structure from FilesData _fdat
    FilesData::Stage const & S = _fdat->Stages[0];

    // Loop along the number of elements
    for (size_t i=0; i<_data->Elements.size(); ++i)
    {
        // Declare a reference to the [i] element of _fdat structure
        FilesData::Element const & E = _fdat->Elements[i];

        // For [i] element, allocate a new element of type equal to Attribute.EleName, ex.: Hex8, Bar2, etc.
        // After this allocation, the array of connectivity will be already set
        _data->Elements[i] = AllocElement(S.Attributes[E.IdxAtt].EleName);
        _data->Elements[i]->Number = i+1;

        // Initialize connectivities of element [i]
        for (size_t j=0; j<E.Connects.size(); ++j)
        {
            // Copy a pointer of node iNode to the internal connects array
            // Remember that inside E.Connects the node numbers start at one (1)
            // This method also setup nodal variables (dx,dy, etc.) for node [j] inside element [i] 
            _data->Elements[i]->SetNode(static_cast<int>(j), &_data->Nodes[E.Connects[j]-1]);
            _data->Nodes[E.Connects[j]-1].Refs()++;
        }
    }

    // Loop along number of embeddeds elements
    if (_fdat->EmbElements.size()>0) _alloc_and_init_embededs();
    
} // }}}
    
inline void StagesManager::_alloc_and_init_embededs() // {{{
{
    // determination of ebedded portions
    //
    // int iElem, fElem, Elem, nElem, lElem; //initial final actual and next and last element  
    FilesData::Stage const & S = _fdat->Stages[0];
    int  n_elems = _fdat->Elements.size();
    int  n_embs  = _fdat->EmbElements.size();
    REAL tiny    = 1E-4;

    for (int i_emb=0; i_emb<n_embs; i_emb++)
    {
        FilesData::Element & great_embed = _fdat->EmbElements[i_emb];
        great_embed.IniVals.resize(1);
        great_embed.IniVals[0].resize(1);
        great_embed.IniVals[0][0]=0.0;
        FilesData::Node & init_node  = _fdat->Nodes[great_embed.Connects[0]-1];
        FilesData::Node & final_node = _fdat->Nodes[great_embed.Connects[1]-1];
        REAL x1 = init_node.X;
        REAL y1 = init_node.Y;
        REAL z1 = init_node.Z;
        REAL x2 = final_node.X;
        REAL y2 = final_node.Y;
        REAL z2 = final_node.Z;

        REAL length = sqrt(pow(x2-x1,2)+pow(y2-y1,2)+pow(z2-z1,2));
        REAL tinylength=0.01*length;
        REAL l = (x2-x1)/length; 
        REAL m = (y2-y1)/length; 
        REAL n = (z2-z1)/length;
        int init_elem   = 0;
        int prev_elem   = 0;
        int actual_elem = 0;
        int next_elem   = 0;
        int final_elem  = 0;
        //find the initial element
        for (int j=0; j<n_elems; j++)
            if (_data->Elements[j]->IsInside(x1+tinylength*l, y1+tinylength*m, z1+tinylength*n)) 
            {
                init_elem = j; break;
            }
        //find the final element
        for (int j=0; j<n_elems; j++)
            if (_data->Elements[j]->IsInside(x2-tinylength*l, y2-tinylength*m, z2-tinylength*n))
            {
                final_elem = j; break;
            }
        REAL xp=x1;
        REAL yp=y1;
        REAL zp=z1;
        REAL x=xp;
        REAL y=yp;
        REAL z=zp;
        actual_elem = prev_elem = init_elem;
        next_elem = 0;
        REAL step  = length;
        bool final = false;

        if (init_elem == final_elem)
        {
            x=x2; y=y2; z=z2; final=true;
        }
        
        //splitting bar:
        do{
            if (!final) {
                step *= 0.501;
                x    += step*l;
                y    += step*m;
                z    += step*n;
                for (int j=0; j<n_elems; j++)
                    if (_data->Elements[j]->IsInside(x, y, z))
                    {
                        actual_elem = j; break;
                    }
            }
            if (prev_elem == actual_elem)
            {
                REAL bf;
                if (!final)
                { // intersection is reached 
                    REAL r, s, t;
                    _data->Elements[next_elem]->InverseMap(x,y,z,r,s,t);
                    bf = fabs(_data->Elements[next_elem]->BoundDistance(r,s,t));
                }
                if ( final || bf <= tiny )  // intersection is reached
                {
                    //Allocate embedded element
                    Element  * elem     = AllocElement(S.Attributes[great_embed.IdxAtt].EleName);
                    Embedded * embedded = dynamic_cast<Embedded*>(elem);

                    //Set connectivities
                    Array<Node*> nodes;
                    size_t n_emb_nodes = elem->nNodes();
                    nodes.resize(n_emb_nodes);

                    // Initialize connectivities of element 
                    size_t n_div = n_emb_nodes - 1;
                    for (size_t j=0; j<n_emb_nodes; ++j) 
                    {
                        nodes[j] = new Node;
                        REAL _x = xp + j*(x-xp)/n_div; 
                        REAL _y = yp + j*(y-yp)/n_div; 
                        REAL _z = zp + j*(z-zp)/n_div; 
                        nodes[j]->Initialize(0, _x, _y, _z); 
                    }

                    // Initialize references of embedded element (embedded connects and tresspassed element)
                    embedded->Initialize(nodes, _data->Elements[actual_elem]);
                    
                    _data->Elements.push_back(embedded);
                    _fdat->Elements.push_back(_fdat->EmbElements[i_emb]);

                    if (final) break; //if this is the final element
                    else{ 
                        xp=x; yp=y; zp=z;
                        if (next_elem==final_elem)
                        {
                            x=x2; y=y2; z=z2;
                            actual_elem = prev_elem = next_elem;
                            final = true;
                        }
                        else
                        {
                            prev_elem = next_elem;
                            step = sqrt(pow(x-x2,2)+pow(y-y2,2)+pow(z-z2,2));
                        }
                    }
                }else{ //when intersection is forward
                    step=fabs(step); //advance
                }
            }else{ //when intersection is backward
                step=-fabs(step);//retrocession
                next_elem = actual_elem;
            }
        }while (true);
    }
} // }}}

inline void StagesManager::_clear_bry_conditions() // {{{
{
    for (size_t i=0; i<_data->Nodes.size(); ++i)
        _data->Nodes[i].ClearBryValues();
} // }}}

inline void StagesManager::_update_nodes_given_face_bry(int iStage) // {{{
{
    // Declare a reference to the iStage structure inside FilesData _fdat
    FilesData::Stage const & S = _fdat->Stages[iStage];

    // Loop along the number of faces
    for (size_t i=0; i<_fdat->Faces.size(); ++i)
    {
        // Index to the FaceBrys array inside a Stage structure.
        int idx_bry = _fdat->Faces[i].IdxBry;
        if (idx_bry!=-1) // Boundary condition prescribed
        {
            // Boundary condition structure
            FilesData::BryCond const & BC = S.FaceBrys[idx_bry];

            // Number of nodes of face [i]
            int n_face_nodes = _fdat->Faces[i].Connects.size(); 

            // Element that contains this face. Remember that inside _fdat OwnerElement numbers start at one (1)
            int k_element = _fdat->Faces[i].OwnerElement-1;

            // Declare a const pointer to the [k_element] element of _data structure
            FEM::Element * const ptrE = _data->Elements[k_element];

            // Temporary array of pointers to the nodes of the face [i]
            Array<FEM::Node*> a_ptr_face_nodes;
            a_ptr_face_nodes.resize(n_face_nodes);
            for (int j=0; j<n_face_nodes; ++j)
            {
                // Remember that inside Connects the node numbers start at one (1)
                int          j_node = _fdat->Faces[i].Connects[j]-1;
                a_ptr_face_nodes[j] = &_data->Nodes[j_node];
            }

            // Loop along prescribed boundary keys
            for (size_t j=0; j<BC.Names.size(); ++j)
            {
                // Check if this key (BC.Names[i]) is essential or not
                if (ptrE->IsEssential(BC.Names[j]))
                {
                    // Loop along the nodes of face [i]
                    for (int m=0; m<n_face_nodes; ++m)
                    {
                        // Get a reference to a DOFVarsStruct
                        FEM::Node::DOFVarsStruct & dof_var = a_ptr_face_nodes[m]->DOFVar(BC.Names[j]);

                        // Setup Essential nodal variable
                        dof_var.EssentialBry = BC.Values[j];
                        dof_var.IsEssenPresc = true;
                    }
                }
                else
                {
                    // Temporary place to put natural boundary name
                    String str_dof_name;

                    // Temporary array to put the nodal boundary values
                    LinAlg::Vector<REAL> a_nodal_values; // will be resized inside CalcFaceNodalValues() method
                    
                    // Calculate the nodal values for face [i] according to the type of the correspondent element
                    ptrE->CalcFaceNodalValues(BC.Names[j]     ,   // In:  Face DOF Name, ex.: "tx, ty, tz"
                                              BC.Values[j]    ,   // In:  Face DOF Value, ex.: 10 kPa, 20 kPa
                                              a_ptr_face_nodes,   // In:  Array of ptrs to face nodes
                                              str_dof_name    ,   // Out: Nodal DOF Name, ex.: "fx, fy, fz"
                                              a_nodal_values  );  // Out: Vector of nodal values

                    // Loop along the nodes of face [i]
                    for (int m=0; m<n_face_nodes; ++m)
                    {
                        // Get a reference to a DOFVarsStruct
                        FEM::Node::DOFVarsStruct & dof_var = a_ptr_face_nodes[m]->DOFVar(str_dof_name);

                        // Setup Natural nodal variable
                          dof_var.NaturalBry += a_nodal_values(m); // Upgrade (add) existent values
                        dof_var.IsEssenPresc  = false;
                    }
                }
            }
        }
        // else: Face does not have prescribed boundary condition
        //       => default case (ux=?,uy=?,fx=0,fy=0) is already set
    }
} // }}}

inline void StagesManager::_update_nodes_given_node_bry(int iStage) // {{{
{
    // Declare a reference to the iStage structure inside FilesData _fdat
    FilesData::Stage const & S = _fdat->Stages[iStage];

    // Loop along the number of nodes
    for (size_t i=0; i<_fdat->Nodes.size(); ++i)
    {
        // Index to the NodeBrys array inside a Stage structure.
        int idx_bry = _fdat->Nodes[i].IdxBry;
        if (idx_bry!=-1 && _data->Nodes[i].Refs()>0) // Boundary condition prescribed
        {
            // Boundary condition structure
            FilesData::BryCond const & BC = S.NodeBrys[idx_bry];

            // Loop along prescribed boundary keys
            for (size_t j=0; j<BC.Names.size(); ++j)
            {
                // Get a reference to a DOFVarsStruct
                FEM::Node::DOFVarsStruct & dof_var = _data->Nodes[i].DOFVar(BC.Names[j]);
                if (_data->Nodes[i].IsEssential(BC.Names[j]))
                {
                    dof_var.EssentialBry = BC.Values[j];
                    dof_var.IsEssenPresc = true;
                }
                else
                {
                      dof_var.NaturalBry += BC.Values[j]; // Upgrade (add) existent values
                    dof_var.IsEssenPresc  = false;
                }
            }
        }
        // else: Nodal does not have prescribed boundary condition
        //       => default case (ux=?,uy=?,fx=0,fy=0) is already set
    }
} // }}}

inline void StagesManager::_set_element_attributes(int iStage) // {{{
{
    // Declare a reference to the iStage structure inside FilesData _fdat
    FilesData::Stage const & S = _fdat->Stages[iStage];

    // Loop along the number of elements
    for (size_t i=0; i<_fdat->Elements.size(); ++i)
    {
        // Index to the Attributes array inside a Stage structure.
        int idx_att = _fdat->Elements[i].IdxAtt;

        // Attribute structure
        FilesData::Attribute const & A = S.Attributes[idx_att];

        // Set Active/Incactive state
        if (A.IsActive) _data->Elements[i]->Activate();
        else            _data->Elements[i]->Deactivate();

        // Set properties
        _data->Elements[i]->SetProperties(A.Properties);

        // Set Model: for [i] element, allocate a new model of type equal to Attribute.ModelName, ex.: Elastic, SubCam, etc.
        _data->Elements[i]->ReAllocateModel(A.ModelName, A.ModelPrms, _fdat->Elements[i].IniVals);
    }
} // }}}

inline void StagesManager::ActivateStage(int iStage) // {{{
{
    _current_stage_num = iStage;
    _clear_bry_conditions();
    _update_nodes_given_face_bry(iStage); // This must be first
    _update_nodes_given_node_bry(iStage); // This must be second => will override (increment) face brys
    _set_element_attributes(iStage);

#ifdef DO_DEBUG
    std::cout << "Stage " << iStage << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
    std::cout << (*_data) << std::endl;
#endif
    
} // }}}

}; // namespace FEM

#endif // MECHSYS_FEM_STAGESMANAGER_H

// vim:fdm=marker

---