modelwrapper.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_MODELWRAPPER_H
#define MECHSYS_MODELWRAPPER_H

// MechSys - models
#include "models/equilibmodel.h"
#include "models/coupled/coupledmodel.h"
#include "models/linearelastic.h"
#include "models/camclay.h"
#include "models/subcam.h"
#include "models/subtij.h"
#include "models/subcamx.h"
#include "models/barcelonax.h"
#include "models/subbar.h"

#ifdef MODELWRAPPER_WITH_VTK
  #include "vtkwrap/vtkwin.h"
  #include "vtkwrap/structgrid.h"
  #include "vtkwrap/sgridisosurf.h"
  #include "vtkwrap/colors.h"
  #include "vtkwrap/plane.h"
  #include "vtkwrap/hedgehog.h"
  #include "vtkwrap/sphere.h"
  #include "tensors/tensors.h"
  #include "tensors/functions.h"
  #include "numerical/meshgrid.h"
#endif // MODELWRAPPER_WITH_VTK

class ModelWrapper
{
public:
    // Constructor
    ModelWrapper(String const & Name, Array<REAL> const & Prms, Array<REAL> const & IniData);

    // Destructor
    ~ModelWrapper();

    // Methods
    void Actualize(Tensor2 const & DSig, REAL const & DPp, Tensor2 & DEps, REAL & DnSr)
    { (_cmodel==NULL ? _emodel->Actualize(DSig,DEps) : _cmodel->Actualize(DSig,DPp,DEps,DnSr)); }

    // Access methods
    //
    String          Name() const { return (_cmodel==NULL ? _emodel->Name() : _cmodel->Name()); }
    Tensor2 const & Sig () const { return (_cmodel==NULL ? _emodel->Sig () : _cmodel->Sig ()); }
    Tensor2 const & Eps () const { return (_cmodel==NULL ? _emodel->Eps () : _cmodel->Eps ()); }
    REAL            Pp  () const { return (_cmodel==NULL ?      0.0        : _cmodel->Pp  ()); }

    void InternalStateValues (Array<REAL>   & IntStateVals ) const
    { (_cmodel==NULL ? _emodel->InternalStateValues(IntStateVals ) : _cmodel->InternalStateValues(IntStateVals));  }

    void InternalStateNames  (Array<String> & IntStateNames) const
    { (_cmodel==NULL ? _emodel->InternalStateNames (IntStateNames) : _cmodel->InternalStateNames (IntStateNames)); }

#ifdef MODELWRAPPER_WITH_VTK
    // VTK
    void AddInipsqActorsTo  (VTKWin & win, int nPts, bool ShowHH, REAL Opac, REAL & p_Max);
    void AddIniStateActorsTo(VTKWin & win, int nPts, bool ShowHH, REAL Opac, REAL & p_Max);
#endif // MODELWRAPPER_WITH_VTK

private:
    // Data
    CoupledModel * _cmodel;
    EquilibModel * _emodel;
    Array<REAL>    _prms;

}; // class ModelWrapper

#define ALLOCMODELFUNC(NAME,PRMS,INIDATA, PTR) (PTR = new ModelWrapper(NAME,PRMS,INIDATA))
#define USE_COUPLEDMODEL
#include "labtestsim/lts.h"
#include "labtestsim/inputdata.h"




inline ModelWrapper::ModelWrapper(String const & Name, Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
    : _cmodel (NULL),
      _emodel (NULL),
      _prms   (Prms)
{
    // Discover if Name is a coupled model or not
    bool is_coupled = false;
    if (EquilibModelFactory.find(Name)==EquilibModelFactory.end())
    {
        if (CoupledModelFactory.find(Name)==CoupledModelFactory.end())
            throw new Fatal(_("ModelWrapper::ModelWrapper: Could not find %s model neither inside EquilibModelFactory nor CoupledModelFactory"),(Name.GetSTL().c_str()));
        is_coupled = true;
    }

    // Allocate the model
    if (is_coupled) _cmodel = AllocCoupledModel(Name,Prms,IniData);
    else
    {
        // Set inidata array
        Array<REAL> inidata;
        inidata.push_back(IniData[0]); // X stress component
        inidata.push_back(IniData[1]); // Y stress component
        inidata.push_back(IniData[2]); // Z stress component
        inidata.push_back(IniData[3]); // vini initial specific volume
        inidata.push_back(IniData[4]); // OCR Over consolidation rate
        _emodel = AllocEquilibModel(Name,Prms,inidata);
    }
} // }}}

inline ModelWrapper::~ModelWrapper() // {{{
{
    if (_cmodel!=NULL) delete _cmodel;
    if (_emodel!=NULL) delete _emodel;
} // }}}

#ifdef MODELWRAPPER_WITH_VTK

inline void ModelWrapper::AddInipsqActorsTo(VTKWin & win, int nPts, bool ShowHH, REAL Opac, REAL & p_Max) // {{{
{
    using Tensors::Stress_p_q_S_sin3th;
    using Tensors::Tensor2;
    using Tensors::Norm;
    using Tensors::Hid2Sig;
    using Tensors::I;

    // Initial state
    Tensor2 sig_ini =  this->Sig();
    REAL    suc_ini = -this->Pp();
    Array<REAL> isv;   this->InternalStateValues(isv);

    if (this->Name()=="BarcelonaX")
    { // {{{

        // Initial data
        REAL k  = _prms[6];
        REAL z0 = isv[0];
        REAL z1 = isv[1];
        REAL M  = dynamic_cast<BarcelonaX*>(_cmodel)->GetM(/*Theta*/Util::ToRad(30.0));
        BarcelonaX::IntVars z; z=z0,z1;
        Tensor2 sig;
        Tensor2 dfdsig;
        REAL    dfdsuc;

        std::cout << "BarcelonaX:  z0=" << z0 << ", z1=" << z1 << std::endl;

        // MeshGrid
        REAL cp  = 1.3;
        REAL cs  = 1.3;
        REAL cq  = 1.2;
        REAL p0m = dynamic_cast<BarcelonaX*>(_cmodel)->Calc_p0(z1,z0);                 // p0 max
        REAL qm  = M*sqrt((z0/2.0+k*z1)*(p0m-z0/2.0)); // q max
        REAL pm  = p0m;                                // p max
        REAL p_min=-cp*k*z1;      p_Max=cp*pm;
        REAL s_min=0.0;      REAL s_max=cs*z1;
        REAL q_min=-cq*qm;   REAL q_max=cq*qm;
        MeshGrid mg(p_min ,p_Max ,nPts,
                    s_min ,s_max ,nPts,
                    q_min ,q_max ,nPts);
        int size = mg.Length();

        // Constants
        REAL        sq2 = sqrt(2.0);
        REAL        sq3 = sqrt(3.0);
        REAL F_ZERO_TOL = 1.0e-2;
        REAL          L = sqrt(z0*z0+z1*z1);

        // Calculate Yield Function Values and Gradients
        double * F = new double [size];
        StructGrid::VectorTuple * T = new StructGrid::VectorTuple [size];
        for (int i=0; i<size; ++i)
        { // {{{
            // Grid values (position)
            REAL   p = mg.X()[i];
            REAL suc = mg.Y()[i];
            REAL   q = mg.Z()[i];
            REAL  th = 0.0;
            REAL s1,s2,s3; // principal values of stress
            Hid2Sig(p,q,th, s1,s2,s3);
            sig=s1,s2,s3,0,0,0;

            // Yield function value
            F[i] = dynamic_cast<BarcelonaX*>(_cmodel)->YieldFunc(sig,z,suc);

            // Derivatives
            dynamic_cast<BarcelonaX*>(_cmodel)->Calc_dFdSig_dFdSuc(sig,z,suc, dfdsig,dfdsuc); // V=dfdsig, S=dfdsuc
            REAL     trV = (dfdsig(0)+dfdsig(1)+dfdsig(2));   // = tr(V) = tr(dfdsig)
            Tensor2 devV = dfdsig - (trV/3.0)*I; 
            REAL    dfdp = trV;
            REAL    dfdq = 2.0*sq3*Norm(devV)/(3.0*sq2);
            if (q<0.0) dfdq=-dfdq;

            // Normal vectors
            REAL Grad[3];
            Grad[0] = dfdp;
            Grad[1] = 0.0; //dfdsuc;
            Grad[2] = dfdq;

            // Normalize
            T[i].vx=Grad[0]; T[i].vy=Grad[1]; T[i].vz=Grad[2];
            if (fabs(F[i])/L<=F_ZERO_TOL)
            {
                REAL norm=1.0;
                if (true) norm = sqrt(Grad[0]*Grad[0]+Grad[1]*Grad[1]+Grad[2]*Grad[2]);
                if (norm/z0>1.0e-5)
                { T[i].vx=Grad[0]/norm; T[i].vy=Grad[1]/norm; T[i].vz=Grad[2]/norm; }
                else
                { T[i].vx=0.0; T[i].vy=0.0; T[i].vz=0.0; }
            }
            else
            { T[i].vx=0.0; T[i].vy=0.0; T[i].vz=0.0; }
        } // }}}

        // Structured grid
        StructGrid sg   (mg.X(),nPts, mg.Y(),nPts, mg.Z(),nPts, F, 0); // p,suc,q, F,T
        StructGrid sg_hh(mg.X(),nPts, mg.Y(),nPts, mg.Z(),nPts, 0, T); // p,suc,q, F,T
        //sg.WriteFile("tmp.sgrid.vtk");

        // LookupTable of colors
        vtkLookupTable * lt = vtkLookupTable::New();
        lt->SetNumberOfColors(10);
        lt->Build();
        for (int i=0; i<10; ++i) lt->SetTableValue(i,CLR["salmon"].C);

        // Generate the IsoSurface
        SGridIsoSurf iso(sg.GetGrid(), /*F*/0.0, lt);
        iso.GetActor()->GetProperty()->SetOpacity(Opac);
        win.AddActor(iso.GetActor());

        /*
        // Vertical plane
        Plane pl(p_min,z1,q_min, 0,1,0, "peacock", 0.3);
        pl.GetObject()->SetPoint1(p_min,z1,q_max);
        pl.GetObject()->SetPoint2(p_Max,z1,q_min);
        //pl.GetObject()->SetPoint2(p_Max/2.0,z1,q_min);
        //win.AddActor(pl.GetActor());
        */

        // Normals
        if (ShowHH)
        {
            REAL normals_sf = 0.2*L;
            HedgeHog * hh = new HedgeHog (sg_hh.GetGrid(), normals_sf);
            hh->GetActor()->GetProperty()->SetColor(CLR["light_blue"].C);
            win.AddActor(hh->GetActor());
            delete hh;
        }

        // Stress/suction points
        REAL p_ini,q_ini,sin3th_ini;
        Tensor2 S_ini;
        Stress_p_q_S_sin3th(sig_ini,p_ini,q_ini,S_ini,sin3th_ini);
        Sphere sp_ini(p_ini,suc_ini,q_ini, 0.03*p_Max, "yellow");
        win.AddActor(sp_ini.GetActor());

        std::cout << "BarcelonaX:  p_ini=" << p_ini << ", q_ini=" << q_ini << ", suc_ini=" << suc_ini << std::endl;

        // Clean up
        lt->Delete();
        delete [] F;

    } // BarcelonaX // }}}

    else if (this->Name()=="SubBar")
    { // {{{

        REAL k  = _prms[6];
        REAL z0 = isv[0];
        REAL z1 = isv[1];
        REAL z2 = isv[2];
        REAL z3 = isv[3];
        REAL M  = dynamic_cast<SubBar*>(_cmodel)->GetM(/*Theta*/Util::ToRad(30.0));
        SubBar::IntVars z; z=z0,z1,z2,z3;
        Tensor2 sig;
        Tensor2 dfdsig;
        REAL    dfdsuc;

        std::cout << "VTKSubBar:  z0=" << z0 << ", z1=" << z1 << ", z2=" << z2 << ", z3=" << z3 << std::endl;

        // MeshGrid
        SubBar::IntVars z_tmp; z_tmp=z2,z3,0,0;
        REAL cp  = 1.1;
        REAL cs  = 1.3;
        REAL cq  = 1.1;
        REAL p0m = dynamic_cast<SubBar*>(_cmodel)->Calc_p0(z_tmp,z3);             // p0 max
        REAL qm  = M*sqrt((z2/2.0+k*z3)*(p0m-z2/2.0)); // q max
        REAL pm  = p0m;                                // p max
        REAL p_min=-cp*k*z3;      p_Max=cp*pm;
        REAL s_min=0.0;      REAL s_max=cs*z3;
        REAL q_min=-cq*qm;   REAL q_max=cq*qm;
        MeshGrid mg(p_min ,p_Max ,nPts,
                    s_min ,s_max ,nPts,
                    q_min ,q_max ,nPts);
        int size = mg.Length();

        // Constants
        REAL        sq2 = sqrt(2.0);
        REAL        sq3 = sqrt(3.0);
        REAL F_ZERO_TOL = 1.0e-2;
        REAL          L = sqrt(z0*z0+z1*z1);

        // Calculate Yield Function Values and Gradients
        double * F_sl = new double [size]; // Subloading
        double * F_no = new double [size]; // Normal
        StructGrid::VectorTuple * T = new StructGrid::VectorTuple [size];
        for (int i=0; i<size; ++i)
        { // {{{
            // Grid values (position)
            REAL   p = mg.X()[i];
            REAL suc = mg.Y()[i];
            REAL   q = mg.Z()[i];
            REAL  th = 0.0;
            REAL s1,s2,s3; // principal values of stress
            Hid2Sig(p,q,th, s1,s2,s3);
            sig=s1,s2,s3,0,0,0;

            // Yield function value
            F_sl[i] = dynamic_cast<SubBar*>(_cmodel)->YieldFunc       (sig,z,suc);
            F_no[i] = dynamic_cast<SubBar*>(_cmodel)->YieldFunc_Normal(sig,z,suc);

            // Derivatives
            dynamic_cast<SubBar*>(_cmodel)->Calc_dFdSig_dFdSuc(sig,z,suc, dfdsig,dfdsuc); // V=dfdsig, S=dfdsuc
            REAL     trV = (dfdsig(0)+dfdsig(1)+dfdsig(2));   // = tr(V) = tr(dfdsig)
            Tensor2 devV = dfdsig - (trV/3.0)*I; 
            REAL    dfdp = trV;
            REAL    dfdq = 2.0*sq3*Norm(devV)/(3.0*sq2);
            if (q<0.0) dfdq = -dfdq;

            // Normal vectors
            REAL Grad[3];
            Grad[0] = dfdp;
            Grad[1] = dfdsuc;
            Grad[2] = dfdq;

            // Normalize
            T[i].vx=Grad[0]; T[i].vy=Grad[1]; T[i].vz=Grad[2];
            if (fabs(F_sl[i])/L<=F_ZERO_TOL)
            {
                REAL norm=1.0;
                if (true) norm = sqrt(Grad[0]*Grad[0]+Grad[1]*Grad[1]+Grad[2]*Grad[2]);
                if (norm/z0>1.0e-5)
                { T[i].vx=Grad[0]/norm; T[i].vy=Grad[1]/norm; T[i].vz=Grad[2]/norm; }
                else
                { T[i].vx=0.0; T[i].vy=0.0; T[i].vz=0.0; }
            }
            else
            { T[i].vx=0.0; T[i].vy=0.0; T[i].vz=0.0; }
        } // }}}

        // Structured grid
        StructGrid sg_sl(mg.X(),nPts, mg.Y(),nPts, mg.Z(),nPts, F_sl, T); // p,suc,q, F,T
        StructGrid sg_no(mg.X(),nPts, mg.Y(),nPts, mg.Z(),nPts, F_no, T); // p,suc,q, F,T

        // LookupTable of colors
        vtkLookupTable * lt_sl = vtkLookupTable::New();
        vtkLookupTable * lt_no = vtkLookupTable::New();
        lt_sl->SetNumberOfColors(10);
        lt_sl->Build();
        lt_no->SetNumberOfColors(10);
        lt_no->Build();
        for (int i=0; i<10; ++i) lt_sl->SetTableValue(i,CLR["forest_green"].C);
        for (int i=0; i<10; ++i) lt_no->SetTableValue(i,CLR["peacock"].C);

        // Generate the IsoSurface
        SGridIsoSurf iso_sl(sg_sl.GetGrid(), /*F*/0.0, lt_sl);
        SGridIsoSurf iso_no(sg_no.GetGrid(), /*F*/0.0, lt_no);
        iso_sl.GetActor()->GetProperty()->SetOpacity(Opac*0.7);
        iso_no.GetActor()->GetProperty()->SetOpacity(Opac);
        win.AddActor(iso_sl.GetActor());
        win.AddActor(iso_no.GetActor());

        /*
        // Vertical plane
        Plane pl(p_min,z3,q_min, 0,1,0, "wheat", 0.5);
        pl.GetObject()->SetPoint1(p_min,z3,q_max);
        pl.GetObject()->SetPoint2(p_Max,z3,q_min);
        //pl.GetObject()->SetPoint2(p_Max/2.0,z1,q_min);
        win.AddActor(pl.GetActor());
        */

        // Stress/suction points
        REAL p_ini,q_ini,sin3th_ini;
        Tensor2 S_ini;
        Stress_p_q_S_sin3th(sig_ini,p_ini,q_ini,S_ini,sin3th_ini);
        Sphere sp_ini(p_ini,suc_ini,q_ini, 0.03*p_Max, "yellow");
        win.AddActor(sp_ini.GetActor());

        std::cout << "VTKSubBar:  p_ini=" << p_ini << ", q_ini=" << q_ini << ", suc_ini=" << suc_ini << std::endl;

        /*
        // Normals
        REAL normals_sf = 0.3*L;
        HedgeHog * hh_sl = new HedgeHog (sg_sl.GetGrid(), normals_sf);
        win.AddActor(hh_sl->GetActor());
        delete hh_sl;
        */

        // Clean up
        delete [] F_sl;
        delete [] F_no;
        lt_sl -> Delete();
        lt_no -> Delete();

    } // SubBar // }}}

} // }}}

inline void ModelWrapper::AddIniStateActorsTo(VTKWin & win, int nPts, bool ShowHH, REAL Opac, REAL & p_Max) // {{{
{
} // }}}

#endif // MODELWRAPPER_WITH_VTK

#endif // MECHSYS_MODELWRAPPER_H

// vim:fdm=marker

---