ysurf.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_YSURF_H
#define MECHSYS_YSURF_H

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

#include <string>

// Tensors
#include "tensors/tensor1.h"
#include "tensors/tensor2.h"
#include "tensors/tensoperators.h"

// VTK
#include "vtkActor.h"

// MechSys
#include "vtkwrap/structgrid.h"
#include "vtkwrap/hedgehog.h"
#include "vtkwrap/sgridoutline.h"
#include "vtkwrap/sgridisosurf.h"
#include "vtkwrap/vtkwin.h"
#include "vtkwrap/axes.h"
#include "vtkwrap/arrow.h"
#include "vtkwrap/colors.h"
#include "vtkwrap/cutclip.h"
#include "vtkwrap/vtkwin.h"
#include "util/util.h"
#include "util/string.h"
#include "constmod/failcrit.h"
#include "numerical/meshgrid.h"
#include "tensors/functions.h"

using Tensors::Sin3Th;
using Tensors::Hid2Sig;
using Util::ToRad;
using Util::Sq3;
using Util::Sq6;

class YSurf
{
    static REAL BIGNUM;
    static REAL SMMINSIG;
    static REAL Q_MIN;
public:
    // Enum
    enum Type {YSurf_DP, YSurf_SM, YSurf_AR};

    // Constructor & Destructor
     YSurf(Type type, FailCrit const & FCrit, REAL MinP=1.0e-3);
    ~YSurf();

    // Config Methdos
    YSurf & SetAlpha (REAL           Alpha) { _alpha     = Alpha;        return (*this); }
    YSurf & SetBeta  (REAL            Beta) { _beta      = Beta;         return (*this); }
    YSurf & SetL     (REAL               L) { _L         = L;            return (*this); }
    YSurf & SetPc    (REAL              Pc) { _Pc        = Pc;           return (*this); }
    YSurf & SetColor (String const & Color) { _ycolor    = Color;        return (*this); }
    YSurf & SetOpac  (REAL         Opacity) { _yopac     = Opacity;      return (*this); }
    YSurf & NormHH   (bool     NormalizeHH) { _hh_norm   = NormalizeHH;  return (*this); }
    YSurf & HHScale  (REAL     ScaleFactor) { _hh_sf     = ScaleFactor;  return (*this); }
    YSurf & HHyfTol  (REAL    YieldFuncTol) { _hh_yf_tol = YieldFuncTol; return (*this); }

    // Config Cut Part of Surfaces
    YSurf & CPartWire  (bool         UseWire) { _cpart_wire  = UseWire; return (*this); }
    YSurf & CPartColor (char const * Color  ) { _cpart_color = Color;   return (*this); }
    YSurf & CPartOpac  (REAL         Opacity) { _cpart_opac  = Opacity; return (*this); }

    // Methods
    REAL Func(REAL SI, REAL SII, REAL SIII, REAL * Grad=NULL); // Grad[3]
    String Name() const;

    // VTK
    vtkActor * GenIsoSurf    (int nPoints);
    vtkActor * GenHedgeHog   (int nPoints);
    CutClip  * AddOctPlane   (REAL p, bool Positive=true);
    void       AddActorsTo   (VTKWin & Win);
    void       DelActorsFrom (VTKWin & Win);
    vtkActor * GetHHActor    () { if (_hh!=NULL) { return _hh->GetActor(); } else { return NULL; } }
    void       DelHedgeHog   () { if (_hh!=NULL) { delete _hh; _hh=NULL; } }

private:
    // Data
    Type             _type;
    FailCrit const & _fcrit;
    REAL             _alpha;
    REAL             _beta;
    REAL             _L;
    REAL             _Pc;
    REAL             _min_P;
    String           _ycolor;
    REAL             _yopac;
    SGridIsoSurf   * _ysgiso;
    HedgeHog       * _hh;
    CutClip        * _cc;
    bool             _cpart_wire;
    String           _cpart_color;
    REAL             _cpart_opac;
    bool             _hh_norm;
    REAL             _hh_sf;
    REAL             _hh_yf_tol;

}; // class YSurf

REAL YSurf::BIGNUM   = 1.0e+10;
REAL YSurf::SMMINSIG = 1.0e-4;
REAL YSurf::Q_MIN    = 1.0e-7;




inline YSurf::YSurf(Type type, FailCrit const & FCrit, REAL MinP) // {{{
    : _type           (type),
      _fcrit          (FCrit), 
      _alpha          (1.0),
      _beta           (0.0),
      _L              (1.0),
      _Pc             (1.0),
      _min_P          (MinP), 
      _ycolor         (String("blue_light")),
      _yopac          (1.0),
      _ysgiso         (NULL), 
      _hh             (NULL),
      _cc             (NULL),
      _cpart_wire     (false),
      _cpart_color    (String("yellow")),
      _cpart_opac     (0.8),
      _hh_norm        (false),
      _hh_sf          (0.05),
      _hh_yf_tol      (0.05)
{
} // }}}

inline YSurf::~YSurf() // {{{
{
    if (_ysgiso !=NULL) delete _ysgiso;
    if (_hh     !=NULL) delete _hh;
    if (_cc     !=NULL) delete _cc;
} // }}}

inline REAL YSurf::Func(REAL SI, REAL SII, REAL SIII, REAL * Grad) // {{{
{
    REAL P,Q,mu;
    REAL      s[3] = {SI,SII,SIII};
    REAL      S[3] = {(2.0*SI-SII-SIII)/3.0, (2.0*SII-SIII-SI)/3.0, (2.0*SIII-SI-SII)/3.0};
    REAL  dP_ds[3] = {0,0,0};
    REAL  dQ_ds[3] = {0,0,0};
    REAL   dmu_dth = 0.0;
    REAL dth_ds[3] = {0,0,0};
    switch (_type)
    {
        case YSurf_DP:
        {
            // Function
            P  = (SI+SII+SIII)/Sq3();
            Q  = sqrt(SI*SI + SII*SII + SIII*SIII - P*P);
            mu = _fcrit.kDP();
            // Grads
            if (Grad!=NULL)
            {
                for (int k=0; k<3; ++k) dP_ds[k] = 1.0/Sq3();  if (Q>Q_MIN)
                for (int k=0; k<3; ++k) dQ_ds[k] = S[k]/Q;
            }
            break;
        }
        case YSurf_SM:
        {
            // Function
            if (SI   <= SMMINSIG) SI  =SMMINSIG;
            if (SII  <= SMMINSIG) SII =SMMINSIG;
            if (SIII <= SMMINSIG) SIII=SMMINSIG;
            REAL I1 = SI+SII+SIII;
            REAL I2 = SI*SII + SII*SIII + SIII*SI;
            REAL I3 = SI*SII*SIII;
                 P  = sqrt(I3/I2)*(sqrt(SI)+sqrt(SII)+sqrt(SIII));
                 Q  = sqrt(SI*SI + SII*SII + SIII*SIII - P*P);
                 mu = _fcrit.kSM();
            // Grads
            if (Grad!=NULL)
            {
                REAL t[3] = {sqrt(s[0]), sqrt(s[1]), sqrt(s[2])};
                REAL  srz = sqrt(I3/I2);
                for (int k=0; k<3; ++k) dP_ds[k] = 0.5*P*(1.0/s[k]-(I1-s[k])/I2) + 0.5*srz/t[k];  if (Q>Q_MIN)
                for (int k=0; k<3; ++k) dQ_ds[k] = (s[k] - P*dP_ds[k])/Q;
            }
            break;
        }
        case YSurf_AR:
        {
            REAL sin3th = Sin3Th(SI,SII,SIII);
            // Function
            P  = (SI+SII+SIII)/Sq3();
            Q  = sqrt(SI*SI + SII*SII + SIII*SIII - P*P);
            mu = _fcrit.AR_Calc_M(sin3th);
            // Grads
            if (Grad!=NULL)
            {
                for (int k=0; k<3; ++k) dP_ds[k] = 1.0/Sq3();  if (Q>Q_MIN)
                for (int k=0; k<3; ++k) dQ_ds[k] = S[k]/Q;
                if (Q>Q_MIN)
                {
                    REAL cos3th = sqrt(1.0-pow(sin3th,2.0));
                    if (cos3th>Q_MIN)
                    {
                               dmu_dth = (0.75*mu*(1.0-_fcrit.AR_w())*cos3th)/(1.0+_fcrit.AR_w()-(1.0-_fcrit.AR_w())*sin3th);
                        REAL     SS[3] = {S[0]*S[0],S[1]*S[1],S[2]*S[2]};
                        REAL      p_SS = (SS[0]+SS[1]+SS[2])/3.0;
                        REAL dev_SS[3] = {SS[0]-p_SS,SS[1]-p_SS,SS[2]-p_SS};   for (int k=0; k<3; ++k)
                             dth_ds[k] = (1.0/(Q*Q*cos3th)) * (Sq6()*dev_SS[k]/Q - sin3th*S[k]);
                    }
                }
            }
            break;
        }
        default:
            throw "YSurf::Func Invalid Type";
    }
    // Grads
    REAL T = _alpha*exp(_beta*P);
    if (Grad!=NULL)
    {
        REAL  dT_dP = _beta*T;
        REAL  dF_dP = 2.0*(P-_Pc)+pow(Q/mu,2.0)*dT_dP;
        REAL  dF_dQ = 2.0*T*Q/(mu*mu);
        REAL dF_dmu = -2.0*T*Q*Q/pow(mu,3.0);
        for (int k=0; k<3; ++k) Grad[k] = dF_dP*dP_ds[k] + dF_dQ*dQ_ds[k] + dF_dmu*dmu_dth*dth_ds[k];
    }
    // Function
    return pow(P-_Pc,2.0) - _L*_L + T*pow(Q/mu,2.0);
} // }}}

inline String YSurf::Name() const // {{{
{
    switch (_type)
    {
        case YSurf_DP: { return String("Drucker-Prager Shaped"      ); }
        case YSurf_SM: { return String("SMP (Novo) Shaped"          ); }
        case YSurf_AR: { return String("Argyris-Sheng et al. Shaped"); }
        default:
            throw "YSurf::Name: Invalid YSurf Type";
    }
} // }}}

inline vtkActor * YSurf::GenIsoSurf(int nPoints) // {{{
{
    // Delete old objects
    if (_ysgiso!=NULL) delete _ysgiso;

    // Create a temporary grid around hydrostatic axis
    REAL minP=_min_P; REAL maxP=_fcrit.Maxp();
    REAL minQ=0.0;    REAL maxQ=_fcrit.Maxp();
    REAL minT=0.0;    REAL maxT=ToRad(360.0);
    MeshGrid mg(minP,maxP,nPoints, minQ,maxQ,nPoints, minT,maxT,nPoints);
    REAL const * P = mg.X();
    REAL const * Q = mg.Y();
    REAL const * T = mg.Z();
    int       Size = mg.Length();

    // Rotate (convert) P,Q,T values into SI,SII,SIII values
    REAL * SI   = new REAL [Size];
    REAL * SII  = new REAL [Size];
    REAL * SIII = new REAL [Size];
    REAL * YF   = new REAL [Size];
    Hid2Sig(P,Q,T, SI,SII,SIII, Size);
    
    // Calculate yield function
    for (int i=0; i<Size; ++i)
        YF[i] = Func(SI[i], SII[i], SIII[i]);

    // Create a VTK structured grid
    StructGrid ysg(SII,nPoints, SIII,nPoints, SI,nPoints, YF);

    // LookupTable for colors
    vtkLookupTable * lt = vtkLookupTable::New();
    lt->SetNumberOfColors(1);
    lt->Build();
    lt->SetTableValue(0,CLR[_ycolor.GetSTL()].C);

    // Generate the Surfaces
    _ysgiso = new SGridIsoSurf (ysg.GetGrid(), 0.0, lt); // for F==0
    _ysgiso->GetActor()->GetProperty()->SetOpacity (_yopac);

    // Clean up
    delete [] SI;
    delete [] SII;
    delete [] SIII;
    delete [] YF;

    // Return
    return _ysgiso->GetActor();
} // }}}

inline vtkActor * YSurf::GenHedgeHog(int nPoints) // {{{
{
    // Delete old objects
    if (_hh!=NULL) delete _hh;

    // Create a temporary grid around hydrostatic axis
    REAL minP=_min_P; REAL maxP=_fcrit.Maxp();
    REAL minQ=0.0;    REAL maxQ=_fcrit.Maxp();
    REAL minT=0.0;    REAL maxT=ToRad(360.0);
    MeshGrid mg(minP,maxP,nPoints, minQ,maxQ,nPoints, minT,maxT,nPoints);
    REAL const * P = mg.X();
    REAL const * Q = mg.Y();
    REAL const * T = mg.Z();
    int       Size = mg.Length();

    // Rotate (convert) P,Q,T values into SI,SII,SIII values
    REAL * SI   = new REAL [Size];
    REAL * SII  = new REAL [Size];
    REAL * SIII = new REAL [Size];
    REAL * YF   = new REAL [Size];
    Hid2Sig(P,Q,T, SI,SII,SIII, Size);
    
    // Calculate yield and gradients
    StructGrid::VectorTuple * V = new StructGrid::VectorTuple [Size];
    for (int i=0; i<Size; ++i)
    {
        REAL Grad[3];
        YF[i] = Func(SI[i], SII[i], SIII[i], Grad);
        if (fabs(YF[i])<=_hh_yf_tol)
        {
            REAL norm=1.0;
            if (_hh_norm) norm = sqrt(Grad[0]*Grad[0]+Grad[1]*Grad[1]+Grad[2]*Grad[2]);
            if (norm>Q_MIN)
            {
                V[i].vx=Grad[1]/norm;
                V[i].vy=Grad[2]/norm;
                V[i].vz=Grad[0]/norm;
            }
            else
            {
                V[i].vx=0.0;
                V[i].vy=0.0;
                V[i].vz=0.0; 
            }
        }
        else
        {
            V[i].vx=0.0;
            V[i].vy=0.0;
            V[i].vz=0.0; 
        }
    }

    // Create a VTK structured grid
    StructGrid ysg(SII,nPoints, SIII,nPoints, SI,nPoints, YF, V);

    // Create a Hedgehog
    _hh = new HedgeHog (ysg.GetGrid(), _hh_sf);

    // Clean up
    delete [] SI;
    delete [] SII;
    delete [] SIII;
    delete [] YF;
    delete [] V;

    // Return
    return _hh->GetActor();
} // }}}

inline CutClip * YSurf::AddOctPlane(REAL p, bool Positive) // {{{
{
    REAL n=1.0;
    if (!Positive) n=-1.0;
    _cc = new CutClip;
    //p,p,p,n,n,n, _ysgiso->GetIsoSurf()->GetOutputPort(), _ycolor, _yopac, _ycolor, _yopac, _cpart_wire, _cpart_color, _cpart_opac);
    return _cc;
} // }}}

inline void YSurf::AddActorsTo(VTKWin & Win) // {{{
{
    if (_ysgiso!=NULL) Win.AddActor(_ysgiso->GetActor());
    if (_hh    !=NULL) Win.AddActor(_hh    ->GetActor());
    if (_cc    !=NULL) _cc->AddActorsTo(Win);
} // }}}

inline void YSurf::DelActorsFrom(VTKWin & Win) // {{{
{
    if (_ysgiso!=NULL) Win.DelActor(_ysgiso->GetActor());
    if (_hh    !=NULL) Win.DelActor(_hh    ->GetActor());
    if (_cc    !=NULL) _cc->DelActorsFrom(Win);
} // }}}

#endif // MECHSYS_YSURF_H

// vim:fdm=marker

---