source: issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp@ 20687

/*!\file Matestar.c
 * \brief: implementation of the Matestar object
 */

#ifdef HAVE_CONFIG_H
        #include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "./Matestar.h"
#include "./Materials.h"
#include "../Inputs/Input.h"
#include "../Inputs/Inputs.h"
#include "../Inputs/TriaInput.h"
#include "../Inputs/PentaInput.h"
#include "../Inputs/ControlInput.h"
#include "../Elements/Element.h"
#include "../Elements/Tria.h"
#include "../Elements/Penta.h"
#include "../Params/Parameters.h"
#include "../Vertex.h"
#include "../Hook.h"
#include "../Node.h"
#include "../IoModel.h"
#include "../../shared/shared.h"

/*Matestar constructors and destructor*/
Matestar::Matestar(){/*{{{*/
        this->helement=NULL;
        this->element=NULL;
        return;
}
/*}}}*/
Matestar::Matestar(int matestar_mid,int index, IoModel* iomodel){/*{{{*/

        /*Intermediaries:*/
        int    matestar_eid;

        /*Initialize id*/
        this->mid=matestar_mid;

        /*Hooks: */
        matestar_eid=index+1;
        this->helement=new Hook(&matestar_eid,1);
        this->element=NULL;

        return;
}
/*}}}*/
Matestar::~Matestar(){/*{{{*/
        delete helement;
        return;
}
/*}}}*/

/*Object virtual functions definitions:*/
Object*   Matestar::copy() {/*{{{*/

        /*Output*/
        Matestar* matestar=NULL;

        /*Initialize output*/
        matestar=new Matestar();

        /*copy fields: */
        matestar->mid=this->mid;
        matestar->helement=(Hook*)this->helement->copy();
        matestar->element =(Element*)this->helement->delivers();

        return matestar;
}
/*}}}*/
Material* Matestar::copy2(Element* element_in) {/*{{{*/

        /*Output*/
        Matestar* matestar=NULL;

        /*Initialize output*/
        matestar=new Matestar();

        /*copy fields: */
        matestar->mid=this->mid;
        matestar->helement=(Hook*)this->helement->copy();
        matestar->element =element_in;

        return matestar;
}
/*}}}*/
void      Matestar::DeepEcho(void){/*{{{*/

        _printf_("Matestar:\n");
        _printf_("   mid: " << mid << "\n");
        _printf_("   element:\n");
        helement->Echo();
}               
/*}}}*/
void      Matestar::Echo(void){/*{{{*/

        _printf_("Matestar:\n");
        _printf_("   mid: " << mid << "\n");
        _printf_("   element:\n");
        helement->Echo();
}
/*}}}*/
int       Matestar::Id(void){ return mid; }/*{{{*/
/*}}}*/
int       Matestar::ObjectEnum(void){/*{{{*/

        return MatestarEnum;

}
/*}}}*/
void      Matestar::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/

        if(marshall_direction==MARSHALLING_BACKWARD)helement=new Hook(); 
        
        MARSHALLING_ENUM(MatestarEnum);
        MARSHALLING(mid);
        this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
        this->element=(Element*)this->helement->delivers();

}
/*}}}*/

/*Matestar management*/
void  Matestar::Configure(Elements* elementsin){/*{{{*/

        /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective 
         * datasets, using internal ids and offsets hidden in hooks: */
        helement->configure((DataSet*)elementsin);
        this->element  = (Element*)helement->delivers();
}
/*}}}*/
void  Matestar::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/

}
/*}}}*/
IssmDouble Matestar::GetA(){/*{{{*/
        _error_("not implemented yet");
}
/*}}}*/
IssmDouble Matestar::GetAbar(){/*{{{*/
        _error_("not implemented yet");
}
/*}}}*/
IssmDouble Matestar::GetB(){/*{{{*/
        _error_("not implemented yet");
}
/*}}}*/
IssmDouble Matestar::GetBbar(){/*{{{*/

        _error_("not implemented yet");
}
/*}}}*/
IssmDouble Matestar::GetN(){/*{{{*/
        _error_("not implemented yet");
}
/*}}}*/
IssmDouble Matestar::GetD(){/*{{{*/
        _error_("not implemented yet");
}
/*}}}*/
IssmDouble Matestar::GetDbar(){/*{{{*/

        _error_("not implemented yet");
}
/*}}}*/
bool Matestar::IsDamage(){/*{{{*/

        _error_("not implemented yet");
}
/*}}}*/
void  Matestar::GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
        _error_("not implemented yet");
        /*From a string tensor and a material object, return viscosity, using Glen's flow law.
                                                                (1-D) B
          viscosity= -------------------------
                                                  2 eps_eff ^[(n-1)/n]

          where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
          and n the flow law exponent.

          If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we 
          return 10^14, initial viscosity.
          */

        /*output: */
        IssmDouble viscosity;

        /*Intermediary: */
        IssmDouble B,D=0.,n;

        /*Get B and n*/
        B=GetB(); _assert_(B>0.);
        n=GetN(); _assert_(n>0.);

        if (n==1.){
                /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
                viscosity=(1.-D)*B/2.;
        }
        else{

                /*if no strain rate, return maximum viscosity*/
                if(eps_eff==0.){
                        viscosity = 1.e+14/2.;
                        //viscosity = B;
                        //viscosity=2.5*pow(10.,17);
                }

                else{
                        viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n));
                }
        }

        /*Checks in debugging mode*/
        if(viscosity<=0) _error_("Negative viscosity");

        /*Return: */
        *pviscosity=viscosity;
}
/*}}}*/
void  Matestar::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){/*{{{*/
        _error_("not implemented yet");
}
/*}}}*/
void  Matestar::GetViscosity_D(IssmDouble* pdmudD,IssmDouble eps_eff){/*{{{*/
         _error_("not implemented yet");
}
/*}}}*/
void  Matestar::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
        _error_("not implemented yet");
}
/*}}}*/
void  Matestar::GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
        _error_("not implemented yet");
}
/*}}}*/
void  Matestar::GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
        _error_("not implemented yet");
}
/*}}}*/
void  Matestar::GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/
        _error_("not implemented yet");
}
/*}}}*/
void  Matestar::GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/
        _error_("not implemented yet");
}
/*}}}*/
void  Matestar::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/

}
/*}}}*/
void  Matestar::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/

}
/*}}}*/
void  Matestar::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols,int name, int type){/*{{{*/
        /*Nothing updated yet*/
}
/*}}}*/
void  Matestar::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/
        /*Nothing updated yet*/
}
/*}}}*/
void  Matestar::InputUpdateFromConstant(int constant, int name){/*{{{*/
        /*Nothing updated yet*/
}
/*}}}*/
void  Matestar::InputUpdateFromConstant(bool constant, int name){/*{{{*/
        /*Nothing updated yet*/
}
/*}}}*/
void  Matestar::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/

        /*Intermediaries*/
        IssmDouble vx,vy,vz;
        IssmDouble dvx[3],dvy[3],dvz[3];
        IssmDouble ko,Ec,Es;

        /*Get velocity derivatives in all directions*/
        _assert_(dim>1);
        _assert_(vx_input);
        vx_input->GetInputValue(&vx,gauss);
        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
        _assert_(vy_input);
        vy_input->GetInputValue(&vy,gauss);
        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
        if(dim==3){
                _assert_(vz_input);
                vz_input->GetInputValue(&vz,gauss);
                vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
        }
        else{
                vz = 0.;
                dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
        }

        /*Get enhancement factors and ko*/
        Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
        Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
        Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
        ec_input->GetInputValue(&Ec,gauss);
        es_input->GetInputValue(&Es,gauss);
        ko_input->GetInputValue(&ko,gauss);

        /*Compute viscosity*/
        *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
}
/*}}}*/
void  Matestar::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
        this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
}/*}}}*/
void  Matestar::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/

        /*Intermediaries*/
        IssmDouble vx,vy,vz;
        IssmDouble dvx[3],dvy[3],dvz[3];
        IssmDouble ko,Ec,Es;

        /*Get velocity derivatives in all directions*/
        _assert_(dim==2 || dim==3);
        _assert_(vx_input);
        vx_input->GetInputValue(&vx,gauss);
        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
        if(dim==3){
                _assert_(vy_input);
                vy_input->GetInputValue(&vy,gauss);
                vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
        }
        else{
                dvx[2] = 0.;
                vy = 0.;
                dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
        }
        vz = 0.;
        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];

        /*Get enhancement factors and ko*/
        Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
        Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
        Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
        ec_input->GetInputValue(&Ec,gauss);
        es_input->GetInputValue(&Es,gauss);
        ko_input->GetInputValue(&ko,gauss);

        /*Compute viscosity*/
        *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
}/*}}}*/
void  Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
        _error_("not implemented yet");
}/*}}}*/
void  Matestar::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
        _error_("not implemented yet");
}/*}}}*/
void  Matestar::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
        /*Intermediaries*/
        IssmDouble vx,vy,vz;
        IssmDouble dvx[3],dvy[3],dvz[3];
        IssmDouble ko,Ec,Es;

        /*Get velocity derivatives in all directions*/
        _assert_(dim==1 || dim==2);
        _assert_(vx_input);
        vx_input->GetInputValue(&vx,gauss);
        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
        if(dim==2){
                _assert_(vy_input);
                vy_input->GetInputValue(&vy,gauss);
                vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
        }
        else{
                dvx[1] = 0.;
                dvx[2] = 0.;
                vy = 0.;
                dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
        }
        vz = 0.;
        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];

        /*Get enhancement factors and ko*/
        Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcbarEnum); _assert_(ec_input);
        Input* es_input = element->inputs->GetInput(MaterialsRheologyEsbarEnum); _assert_(es_input);
        Input* ko_input = element->inputs->GetInput(MaterialsRheologyKobarEnum); _assert_(ko_input);
        ec_input->GetInputValue(&Ec,gauss);
        es_input->GetInputValue(&Es,gauss);
        ko_input->GetInputValue(&ko,gauss);

        /*Compute viscosity*/
        *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
}/*}}}*/
void  Matestar::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
        _error_("not implemented yet");
}/*}}}*/
void  Matestar::ResetHooks(){/*{{{*/

        this->element=NULL;

        /*Get Element type*/
        this->helement->reset();

}
/*}}}*/
IssmDouble Matestar::GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/

        /*Intermediaries*/
        IssmDouble viscosity;
        IssmDouble E,lambdas;
        IssmDouble epso,epsprime_norm;
        IssmDouble vmag,dvmag[3];

        /*Calculate velocity magnitude and its derivative*/
        vmag = sqrt(vx*vx+vy*vy+vz*vz);
        if(vmag<1e-12){
                dvmag[0]=0;
                dvmag[1]=0;
                dvmag[2]=0;
        }
        else{
                dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]);
                dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]);
                dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
        }

        EstarStrainrateQuantities(&epso,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
        lambdas=EstarLambdaS(epso,epsprime_norm);

        /*Get total enhancement factor E(lambdas)*/
        E = Ec + (Es-Ec)*lambdas*lambdas;

        /*Compute viscosity*/
        _assert_(E>0.); 
        _assert_(ko>0.); 
        _assert_(epso>0.); 
        viscosity = 1./(2*pow(ko*E*epso*epso,1./3.));

        /*Assign output pointer*/
        return viscosity;
}
/*}}}*/