/home/dorival/mechsys/lib/models/model.h

Go to the documentation of this file.

/************************************************************************
 * MechSys - Open Library for Mechanical Systems                        *
 * Copyright (C) 2005 Dorival M. Pedroso, Raul Durand                   *
 * Copyright (C) 2009 Sergio Galindo                                    *
 *                                                                      *
 * This program 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 3 of the License, or    *
 * any later version.                                                   *
 *                                                                      *
 * This program 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 this program. If not, see <http://www.gnu.org/licenses/>  *
 ************************************************************************/

#ifndef MECHSYS_MODEL_H
#define MECHSYS_MODEL_H

// Std Lib
#include <iostream>
#include <fstream>

// MechSys
#include <mechsys/geomtype.h>
#include <mechsys/util/maps.h>
#include <mechsys/util/fatal.h>
#include <mechsys/linalg/matvec.h>
#include <mechsys/numerical/odesolver.h>

class State
{
public:
    State (int NDim) {}
    virtual ~State () {}
    virtual void   Init    (SDPair const & Ini, size_t NIvs=0) =0;
    virtual void   Backup  () =0;
    virtual void   Restore () =0;
    virtual size_t PckSize ()                        const =0; 
    virtual void   Pack    (Array<double>       & V) const =0; 
    virtual void   Unpack  (Array<double> const & V)       =0; 
    virtual void   Output  (SDPair & KeysVals)       const =0; 
};

#include <mechsys/models/equilibstate.h>

class Model
{
public:
    // static
    static Vec_t I;     
    static Mat_t Psd;   
    static Mat_t IdyI;  

    // Constructor & Destructor
    Model (int NDim, SDPair const & Prms, size_t NIvs=0, char const * Name="__unnamed_model__"); 
    virtual ~Model () {}

    // Methods
    virtual void   InitIvs      (SDPair const & Ini, State * Sta)                              const =0;
    virtual void   Rate         (State const * Sta, Vec_t const & DEpsDt, Vec_t & DSigDt, Vec_t & DIvsDt) const { throw new Fatal("Model::Rate: This method is not available in this model (%s)",Name.CStr()); }
    virtual void   TgIncs       (State const * Sta, Vec_t & DEps, Vec_t & DSig, Vec_t & DIvs)  const { throw new Fatal("Model::TgIncs: This method is not available in this model (%s)",Name.CStr()); }
    virtual void   UpdateSta    (ATensor2 const & F, State * Sta)                              const { throw new Fatal("Model::UpdateSta: This method is not available in this model (%s)",Name.CStr()); }
    virtual void   InvTgIncs    (State const * Sta, Vec_t & DSig, Vec_t & DEps, Vec_t & DIvs)  const { throw new Fatal("Model::InvTgIncs: This method is not available in this model (%s)",Name.CStr()); }
    virtual void   Stiffness    (State const * Sta, Mat_t & D)                                 const { throw new Fatal("Model::Stiffness: This method is not available in this model (%s)",Name.CStr()); }
    virtual void   Stiffness    (State const * Sta, Mat_t & D, Vec_t & Dw)                     const { throw new Fatal("Model::Stiffness(Dw): This method is not available in this model (%s)",Name.CStr()); }
    virtual void   Hydraulic    (State const * Sta, Mat_t & Kw, double & ChiW, double & InvQs) const { throw new Fatal("Model::Hydraulic: This method is not available in this model (%s)",Name.CStr()); }
    virtual bool   LoadCond     (State const * Sta, Vec_t const & DEps, double & alpInt)       const { alpInt=-1; return false; }
    virtual size_t CorrectDrift (State * Sta)                                                  const { return 0; }
    virtual size_t CorrectDrift (State * Sta, double & pw)                                     const { return 0; }
    virtual void   UpdatePath   (State const * Sta, Vec_t const & DEps, Vec_t const & DSig)    const {}

    // Methods: HydroMech
    virtual void   TgIncs       (State const * Sta, double pw, Vec_t & DEps, double Dpw, Vec_t & DSig, Vec_t & DIvs)  const { throw new Fatal("Model::TgIncs(HydroMech): This method is not available in this model (%s)",Name.CStr()); }
    virtual void   Stiffness    (State const * Sta, double pw, Mat_t & D, Vec_t & d)                                  const { throw new Fatal("Model::Stiffness(HydroMech): This method is not available in this model (%s)",Name.CStr()); }
    virtual bool   LoadCond     (State const * Sta, double pw, Vec_t const & DEps, double Dpw)                        const { throw new Fatal("Model::LoadCond(HydroMech): This method is not available in this model (%s)",Name.CStr()); }

    // Data
    int           NDim;    
    SDPair        Prms;    
    GeomType      GTy;     
    String        Name;    
    size_t        NCps;    
    size_t        NIvs;    
    Array<String> IvNames; 
    bool          HMCoup;  

    // Data to be set by derived classes
    bool UseUpdateSta; 

    // General constants
    double Grav;     
    double GamW;     
    double Rho;      
    double RhoS;     
    double Por;      
    double GamNat;   
    double GamSat;   

#define STRESSUPDATE_DECLARE
    #include <mechsys/models/stressupdate.h>
    mutable StressUpdate SUp;
#undef STRESSUPDATE_DECLARE
};

Vec_t Model::I;
Mat_t Model::Psd;
Mat_t Model::IdyI;




inline Model::Model (int TheNDim, SDPair const & ThePrms, size_t NIv, char const * TheName)
    : NDim(TheNDim), Prms(ThePrms), GTy(SDPairToGType(ThePrms,(TheNDim==3?"d3d":"d2d"))), 
      Name(TheName), NCps(2*NDim),  NIvs(NIv), HMCoup(false), UseUpdateSta(false)
{
    // constants
    Grav   = (Prms.HasKey("grav")   ? Prms("grav")   : -1);
    GamW   = (Prms.HasKey("gamW")   ? Prms("gamW")   : -1);
    Rho    = (Prms.HasKey("rho")    ? Prms("rho")    : -1);
    RhoS   = (Prms.HasKey("rhoS")   ? Prms("rhoS")   : -1);
    Por    = (Prms.HasKey("por")    ? Prms("por")    : -1);
    GamNat = (Prms.HasKey("gamNat") ? Prms("gamNat") : -1);
    GamSat = (Prms.HasKey("gamSat") ? Prms("gamSat") : -1);
    if (Prms.HasKey("gamNat") && !Prms.HasKey("gamSat")) GamSat = GamNat;
    if (Prms.HasKey("gamSat") && !Prms.HasKey("gamNat")) GamNat = GamSat;
    if (Prms.HasKey("rhoS"))
    {
        double rho_w = Prms("gamW") / Prms("grav");
        double n     = Prms("por");
        Rho = n*rho_w + (1.0-n)*Prms("rhoS");
        if (!Prms.HasKey("gamNat")) GamNat = Rho * Grav;
        if (!Prms.HasKey("gamSat")) GamSat = GamNat;
    }

    // stress update
    SUp.SetModel   (this);
    IvNames.Resize (NIv);

    // static variables
    if (size(I)==0)
    {
        I.change_dim(NCps);
        I(0) = 1.;
        I(1) = 1.;
        I(2) = 1.;
        Calc_Psd  (NCps, Psd);
        Calc_IdyI (NCps, IdyI);
    }
}

std::ostream & operator<< (std::ostream & os, Model const & D)
{
    os << D.Name << " " << D.NDim << "D " << GTypeToStr(D.GTy) << std::endl;
    os << "  Prms = {" << D.Prms << "}\n";
    os << "  Gravity                               : grav   = " << D.Grav   << std::endl;
    os << "  Unit weight of water                  : gamW   = " << D.GamW   << std::endl;
    os << "  Density of the material               : rho    = " << D.Rho    << std::endl;
    os << "  Density of solids                     : rhoS   = " << D.RhoS   << std::endl;
    os << "  Porosity of material                  : por    = " << D.Por    << std::endl;
    os << "  Natural unit weight of geo-material   : gamNat = " << D.GamNat << std::endl;
    os << "  Saturated unit weight of geo-material : gamSat = " << D.GamSat << std::endl;
    return os;
}


#define STRESSUPDATE_IMPLEMENT
  #include <mechsys/models/stressupdate.h>
#undef STRESSUPDATE_IMPLEMENT




SDPair MODEL; 

typedef std::map<String,Array<String> > Str2ArrayStr_t; 

Str2ArrayStr_t MODEL_PRM_NAMES; 
Str2ArrayStr_t MODEL_IVS_NAMES; 

typedef Model * (*ModelMakerPtr)(int NDim, SDPair const & Prms, Model const * AnotherMdl);

typedef std::map<String, ModelMakerPtr> ModelFactory_t;

ModelFactory_t ModelFactory;

Model * AllocModel(String const & Name, int NDim, SDPair const & Prms, Model const * Mdl)
{
    ModelFactory_t::iterator it = ModelFactory.find(Name);
    if (it==ModelFactory.end()) throw new Fatal("AllocModel: '%s' is not available", Name.CStr());

    Model * ptr = (*it->second)(NDim,Prms,Mdl);

    return ptr;
}

Model * AllocModel(double IDinMODEL, int NDim, SDPair const & Prms, Model const * Mdl)
{
    String model_name;
    MODEL.Val2Key (IDinMODEL, model_name);
    return AllocModel (model_name, NDim, Prms, Mdl);
}

Array<String> MODEL_CTE_NAMES; 

int ModelRegister()
{
    MODEL_CTE_NAMES.Resize(7);
    MODEL_CTE_NAMES = "grav", "gamW", "rho", "rhoS", "por", "gamNat", "gamSat";
    return 0;
}

int __Model_dummy_int = ModelRegister();




#ifdef USE_BOOST_PYTHON
double PyMODEL (BPy::str const & Key) { return MODEL(BPy::extract<char const *>(Key)()); }
#endif

#endif // MECHSYS_MODEL_H