Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/classes/objects/Materials/Matdamageice.cpp =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/classes/objects/Materials/Matdamageice.cpp (revision 0) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/classes/objects/Materials/Matdamageice.cpp (revision 13130) @@ -0,0 +1,896 @@ +/*!\file Matdamageice.c + * \brief: implementation of the Matdamageice object + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include +#include +#include "../../classes.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../shared/shared.h" +#include "../../../include/include.h" + +/*Matdamageice constructors and destructor*/ +/*FUNCTION Matdamageice::Matdamageice(){{{*/ +Matdamageice::Matdamageice(){ + this->inputs=NULL; + this->helement=NULL; + return; +} +/*}}}*/ +/*FUNCTION Matdamageice::Matdamageice(int id, int index, IoModel* iomodel, int num_vertices){{{*/ +Matdamageice::Matdamageice(int matice_mid,int index, IoModel* iomodel){ + + /*Intermediaries:*/ + int i; + int matice_eid; + + /*Initialize id*/ + this->mid=matice_mid; + + /*Initialize inputs*/ + this->inputs=new Inputs(); + + /*Initialize inputs from IoModel*/ + this->InputUpdateFromIoModel(index,iomodel); + + /*Hooks: */ + matice_eid=index+1; + this->helement=new Hook(&matice_eid,1); + + return; + +} +/*}}}*/ +/*FUNCTION Matdamageice::~Matdamageice(){{{*/ +Matdamageice::~Matdamageice(){ + delete helement; + delete inputs; + return; +} +/*}}}*/ + +/*Object virtual functions definitions:*/ +/*FUNCTION Matdamageice::Echo {{{*/ +void Matdamageice::Echo(void){ + + _printLine_("Matdamageice:"); + _printLine_(" mid: " << mid); + _printLine_(" inputs:"); + inputs->Echo(); + _printLine_(" element:"); + helement->Echo(); +} +/*}}}*/ +/*FUNCTION Matdamageice::DeepEcho {{{*/ +void Matdamageice::DeepEcho(void){ + + _printLine_("Matdamageice:"); + _printLine_(" mid: " << mid); + _printLine_(" inputs:"); + inputs->DeepEcho(); + _printLine_(" element:"); + helement->Echo(); +} +/*}}}*/ +/*FUNCTION Matdamageice::Id {{{*/ +int Matdamageice::Id(void){ return mid; } +/*}}}*/ +/*FUNCTION Matdamageice::MyRank {{{*/ +int Matdamageice::MyRank(void){ + extern int my_rank; + return my_rank; +} +/*}}}*/ +/*FUNCTION Matdamageice::ObjectEnum{{{*/ +int Matdamageice::ObjectEnum(void){ + + return MatdamageiceEnum; + +} +/*}}}*/ +/*FUNCTION Matdamageice::copy {{{*/ +Object* Matdamageice::copy() { + + /*Output*/ + Matdamageice* matice=NULL; + + /*Initialize output*/ + matice=new Matdamageice(); + + /*copy fields: */ + matice->mid=this->mid; + matice->helement=(Hook*)this->helement->copy(); + if(this->inputs) matice->inputs=(Inputs*)this->inputs->Copy(); + else matice->inputs=new Inputs(); + + return matice; +} +/*}}}*/ + +/*Matdamageice management*/ +/*FUNCTION Matdamageice::Configure {{{*/ +void Matdamageice::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(elementsin); +} +/*}}}*/ +/*FUNCTION Matdamageice::SetCurrentConfiguration {{{*/ +void Matdamageice::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){ + +} +/*}}}*/ +/*FUNCTION Matdamageice::GetA {{{*/ +IssmDouble Matdamageice::GetA(){ + /* + * A = 1/B^n + */ + + IssmDouble B,n; + + inputs->GetInputAverage(&B,MaterialsRheologyBEnum); + n=this->GetN(); + + return pow(B,-n); +} +/*}}}*/ +/*FUNCTION Matdamageice::GetB {{{*/ +IssmDouble Matdamageice::GetB(){ + + /*Output*/ + IssmDouble B; + + inputs->GetInputAverage(&B,MaterialsRheologyBEnum); + return B; +} +/*}}}*/ +/*FUNCTION Matdamageice::GetBbar {{{*/ +IssmDouble Matdamageice::GetBbar(){ + + /*Output*/ + IssmDouble Bbar; + + inputs->GetInputAverage(&Bbar,MaterialsRheologyBbarEnum); + return Bbar; +} +/*}}}*/ +/*FUNCTION Matdamageice::GetN {{{*/ +IssmDouble Matdamageice::GetN(){ + + /*Output*/ + IssmDouble n; + + inputs->GetInputAverage(&n,MaterialsRheologyNEnum); + return n; +} +/*}}}*/ +/*FUNCTION Matdamageice::GetZ {{{*/ +IssmDouble Matdamageice::GetZ(){ + + /*Output*/ + IssmDouble Z; + + inputs->GetInputAverage(&Z,MaterialsRheologyZEnum); + return Z; +} +/*}}}*/ +/*FUNCTION Matdamageice::GetZbar {{{*/ +IssmDouble Matdamageice::GetZbar(){ + + /*Output*/ + IssmDouble Zbar; + + inputs->GetInputAverage(&Zbar,MaterialsRheologyZbarEnum); + return Zbar; +} +/*}}}*/ +/*FUNCTION Matdamageice::GetVectorFromInputs{{{*/ +void Matdamageice::GetVectorFromInputs(Vector* vector,int input_enum){ + + /*Intermediaries*/ + Element *element= NULL; + + /*Recover element*/ + element=(Element*)helement->delivers(); + + /*Check that input_enum is a material input*/ + if (!IsInput(input_enum)) return; + + switch(element->ObjectEnum()){ + + case TriaEnum:{ + + /*Prepare index list*/ + int doflist1[3]; + for(int i=0;i<3;i++) doflist1[i]=((Tria*)element)->nodes[i]->GetVertexDof(); + + /*Get input (either in element or material)*/ + Input* input=inputs->GetInput(input_enum); + if(!input) _error_("Input " << EnumToStringx(input_enum) << " not found in material"); + + /*We found the enum. Use its values to fill into the vector, using the vertices ids: */ + input->GetVectorFromInputs(vector,&doflist1[0]);} + break; + + default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet"); + } +} +/*}}}*/ +/*FUNCTION Matdamageice::GetViscosity2d {{{*/ +void Matdamageice::GetViscosity2d(IssmDouble* pviscosity, IssmDouble* epsilon){ + /*From a string tensor and a material object, return viscosity, using Glen's flow law. + Z * B + viscosity= ------------------------------------------------------------------- + 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] + + where viscosity is the viscotiy, B the flow law parameter , (u,v) the velocity + vector, and n the flow law exponent. + + If epsilon is NULL, it means this is the first time SystemMatrices is being run, and we + return 10^14, initial viscosity. + */ + + /*output: */ + IssmDouble viscosity; + + /*input strain rate: */ + IssmDouble exx,eyy,exy; + + /*Intermediary: */ + IssmDouble A,e; + IssmDouble Btmp,B,n,Z; + + /*Get B and n*/ + Btmp=GetBbar(); + Z=GetZbar(); + n=GetN(); + B=Z*Btmp; + + if (n==1){ + /*Viscous behaviour! viscosity=B: */ + viscosity=B/2; + } + else{ + if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){ + viscosity=0.5*pow((IssmDouble)10,(IssmDouble)14); + } + else{ + /*Retrive strain rate components: */ + exx=*(epsilon+0); + eyy=*(epsilon+1); + exy=*(epsilon+2); + + /*Build viscosity: viscosity=B/(2*A^e) */ + A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy; + if(A==0){ + /*Maxiviscositym viscosity for 0 shear areas: */ + viscosity=2.5*pow(10.,17.); + } + else{ + e=(n-1)/(2*n); + viscosity=B/(2*pow(A,e)); + } + } + } + + /*Checks in debugging mode*/ + if(viscosity<=0) _error_("Negative viscosity"); + _assert_(B>0); + _assert_(n>0); + + /*Return: */ + *pviscosity=viscosity; +} +/*}}}*/ +/*FUNCTION Matdamageice::GetViscosity3d {{{*/ +void Matdamageice::GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* epsilon){ + + /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: + * + * B + * viscosity3d= ------------------------------------------------------------------- + * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] + * + * where mu is the viscotiy, B the flow law parameter , (u,v) the velocity + * vector, and n the flow law exponent. + * + * If epsilon is NULL, it means this is the first time Emg is being run, and we + * return g, initial viscosity. + */ + + /*output: */ + IssmDouble viscosity3d; + + /*input strain rate: */ + IssmDouble exx,eyy,exy,exz,eyz; + + /*Intermediaries: */ + IssmDouble A,e; + IssmDouble B,n,Z; + + /*Get B, Z and n*/ + n=GetN(); + Z=GetZ(); + B=Z*GetB(); + + if (n==1){ + /*Viscous behaviour! viscosity3d=B: */ + viscosity3d=B/2; + } + else{ + if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && + (epsilon[3]==0) && (epsilon[4]==0)){ + viscosity3d=0.5*pow((IssmDouble)10,(IssmDouble)14); + } + else{ + + /*Retrive strain rate components: */ + exx=*(epsilon+0); + eyy=*(epsilon+1); + exy=*(epsilon+2); + exz=*(epsilon+3); + eyz=*(epsilon+4); + + /*Build viscosity: viscosity3d=2*B/(2*A^e) */ + A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy; + if(A==0){ + /*Maxiviscosity3dm viscosity for 0 shear areas: */ + viscosity3d=2.25*pow((IssmDouble)10,(IssmDouble)17); + } + else{ + e=(n-1)/2/n; + + viscosity3d=B/(2*pow(A,e)); + } + } + } + + /*Checks in debugging mode*/ + if(viscosity3d<=0) _error_("Negative viscosity"); + _assert_(B>0); + _assert_(n>0); + + /*Assign output pointers:*/ + *pviscosity3d=viscosity3d; +} +/*}}}*/ +/*FUNCTION Matdamageice::GetViscosity3dStokes {{{*/ +void Matdamageice::GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon){ + /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: + * + * B + * viscosity3d= ------------------------------------------------------------------- + * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] + * + * where mu is the viscotiy, B the flow law parameter , (u,v) the velocity + * vector, and n the flow law exponent. + * + * If epsilon is NULL, it means this is the first time Emg is being run, and we + * return g, initial viscosity. + */ + + /*output: */ + IssmDouble viscosity3d; + + /*input strain rate: */ + IssmDouble exx,eyy,exy,exz,eyz,ezz; + + /*Intermediaries: */ + IssmDouble A,e; + IssmDouble B,n,Z; + IssmDouble eps0; + + /*Get B and n*/ + eps0=pow((IssmDouble)10,(IssmDouble)-27); + n=GetN(); + Z=GetZ(); + B=Z*GetB(); + + if (n==1){ + /*Viscous behaviour! viscosity3d=B: */ + viscosity3d=B/2; + } + else{ + if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && + (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){ + viscosity3d=0.5*pow((IssmDouble)10,(IssmDouble)14); + } + else{ + + /*Retrive strain rate components: */ + exx=*(epsilon+0); + eyy=*(epsilon+1); + ezz=*(epsilon+2); //not used + exy=*(epsilon+3); + exz=*(epsilon+4); + eyz=*(epsilon+5); + + /*Build viscosity: viscosity3d=B/(2*A^e) */ + A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy+pow(eps0,2); + if(A==0){ + /*Maxiviscosity3dm viscosity for 0 shear areas: */ + viscosity3d=2.25*pow((IssmDouble)10,(IssmDouble)17); + } + else{ + e=(n-1)/2/n; + viscosity3d=B/(2*pow(A,e)); + } + } + } + + /*Checks in debugging mode*/ + if(viscosity3d<=0) _error_("Negative viscosity"); + _assert_(B>0); + _assert_(n>0); + + /*Assign output pointers:*/ + *pviscosity3d=viscosity3d; +} +/*}}}*/ +/*FUNCTION Matdamageice::GetViscosityComplement {{{*/ +void Matdamageice::GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){ + /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: + * + * 1 + * viscosity= ------------------------------------------------------------------- + * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] + * + * If epsilon is NULL, it means this is the first time Gradjb is being run, and we + * return mu20, initial viscosity. + */ + + /*output: */ + IssmDouble viscosity_complement; + + /*input strain rate: */ + IssmDouble exx,eyy,exy; + + /*Intermediary value A and exponent e: */ + IssmDouble A,e; + IssmDouble B,n; + + /*Get B and n*/ + B=GetBbar(); + n=GetN(); + + if(epsilon){ + exx=*(epsilon+0); + eyy=*(epsilon+1); + exy=*(epsilon+2); + + /*Build viscosity: mu2=B/(2*A^e) */ + A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy; + if(A==0){ + /*Maximum viscosity_complement for 0 shear areas: */ + viscosity_complement=2.25*pow((IssmDouble)10,(IssmDouble)17); + } + else{ + e=(n-1)/(2*n); + + viscosity_complement=1/(2*pow(A,e)); + } + } + else{ + viscosity_complement=4.5*pow((IssmDouble)10,(IssmDouble)17); + } + + /*Checks in debugging mode*/ + _assert_(B>0); + _assert_(n>0); + _assert_(viscosity_complement>0); + + /*Return: */ + *pviscosity_complement=viscosity_complement; +} +/*}}}*/ +/*FUNCTION Matdamageice::GetViscosityZComplement {{{*/ +void Matdamageice::GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){ + /*Return viscosity derivative for control method d(mu)/dZ: + * + * B + * dviscosity= ------------------------------------------------------------------- + * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] + * + * If epsilon is NULL, it means this is the first time Gradjb is being run, and we + * return mu20, initial viscosity. + */ + + /*output: */ + IssmDouble viscosity_complement; + + /*input strain rate: */ + IssmDouble exx,eyy,exy; + + /*Intermediary value A and exponent e: */ + IssmDouble A,e; + IssmDouble B,n; + + /*Get B and n*/ + B=GetBbar(); + n=GetN(); + + if(epsilon){ + exx=*(epsilon+0); + eyy=*(epsilon+1); + exy=*(epsilon+2); + + /*Build viscosity: mu2=B/(2*A^e) */ + A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy; + if(A==0){ + /*Maximum viscosity_complement for 0 shear areas: */ + viscosity_complement=2.25*pow((IssmDouble)10,(IssmDouble)17); + } + else{ + e=(n-1)/(2*n); + + viscosity_complement=B/(2*pow(A,e)); + } + } + else{ + viscosity_complement=4.5*pow((IssmDouble)10,(IssmDouble)17); + } + + /*Checks in debugging mode*/ + _assert_(B>0); + _assert_(n>0); + _assert_(viscosity_complement>0); + + /*Return: */ + *pviscosity_complement=viscosity_complement; +} +/*}}}*/ +/*FUNCTION Matdamageice::GetViscosityDerivativeEpsSquare{{{*/ +void Matdamageice::GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){ + + /*output: */ + IssmDouble mu_prime; + IssmDouble mu,n,eff2; + + /*input strain rate: */ + IssmDouble exx,eyy,exy,exz,eyz; + + /*Get visocisty and n*/ + GetViscosity3d(&mu,epsilon); + n=GetN(); + + if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && + (epsilon[3]==0) && (epsilon[4]==0)){ + mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)14); + } + else{ + /*Retrive strain rate components: */ + exx=epsilon[0]; + eyy=epsilon[1]; + exy=epsilon[2]; + exz=epsilon[3]; + eyz=epsilon[4]; + eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz; + + mu_prime=(1-n)/(2*n) * mu/eff2; + } + + /*Assign output pointers:*/ + *pmu_prime=mu_prime; +} +/*}}}*/ +/*FUNCTION Matdamageice::GetViscosity2dDerivativeEpsSquare{{{*/ +void Matdamageice::GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){ + + /*output: */ + IssmDouble mu_prime; + IssmDouble mu,n,eff2; + + /*input strain rate: */ + IssmDouble exx,eyy,exy,exz; + + /*Get visocisty and n*/ + GetViscosity2d(&mu,epsilon); + n=GetN(); + + if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){ + mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)14); + } + else{ + /*Retrive strain rate components: */ + exx=epsilon[0]; + eyy=epsilon[1]; + exy=epsilon[2]; + eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy ; + + mu_prime=(1-n)/(2*n) * mu/eff2; + } + + /*Assign output pointers:*/ + *pmu_prime=mu_prime; +} +/*}}}*/ +/*FUNCTION Matdamageice::InputDuplicate{{{*/ +void Matdamageice::InputDuplicate(int original_enum,int new_enum){ + + /*Call inputs method*/ + if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum); + +} +/*}}}*/ +/*FUNCTION Matdamageice::InputUpdateFromVector(IssmDouble* vector, int name, int type) {{{*/ +void Matdamageice::InputUpdateFromVector(IssmDouble* vector, int name, int type){ + + /*Intermediaries*/ + Element *element = NULL; + + /*Recover element*/ + element=(Element*)helement->delivers(); + + /*Check that name is an element input*/ + if (!IsInput(name)) return; + + switch(type){ + + case VertexEnum: + + switch(element->ObjectEnum()){ + + case TriaEnum: { + IssmDouble values[3]; + for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->nodes[i]->GetVertexDof()]; + this->inputs->AddInput(new TriaP1Input(name,values)); + return; + } + default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet"); + } + default: _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); + } +} +/*}}}*/ +/*FUNCTION Matdamageice::InputUpdateFromVector(int* vector, int name, int type) {{{*/ +void Matdamageice::InputUpdateFromVector(int* vector, int name, int type){ + /*Nothing updated yet*/ +} +/*}}}*/ +/*FUNCTION Matdamageice::InputUpdateFromVector(bool* vector, int name, int type) {{{*/ +void Matdamageice::InputUpdateFromVector(bool* vector, int name, int type){ + /*Nothing updated yet*/ +} +/*}}}*/ +/*FUNCTION Matdamageice::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type) {{{*/ +void Matdamageice::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){ + + /*Intermediaries*/ + Element *element = NULL; + Parameters* parameters= NULL; + int dim; + + /*Recover element*/ + element=(Element*)helement->delivers(); + + /*Check that name is an element input*/ + if (!IsInput(name)) return; + + switch(type){ + + case VertexEnum: + + switch(element->ObjectEnum()){ + + case TriaEnum: { + IssmDouble values[3]; + for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->nodes[i]->GetSidList()]; //use sid list, to index into serial oriented vector + this->inputs->AddInput(new TriaP1Input(name,values)); + /*Special case for rheology B in 2D: Pourave land for this solution{{{*/ + if(name==MaterialsRheologyBEnum){ + /*Are we in 2D?:*/ + if(element->ObjectEnum()==TriaEnum){ + parameters=((Tria*)(element))->parameters; + } + else{ + parameters=((Penta*)(element))->parameters; + } + parameters->FindParam(&dim,MeshDimensionEnum); + if(dim==2){ + /*Dupliacte rheology input: */ + this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,values)); + } + } + /*}}}*/ + return; + } + default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet"); + } + default: _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); + } + + + +} +/*}}}*/ +/*FUNCTION Matdamageice::InputUpdateFromMatrixDakota(int* vector, int name, int type) {{{*/ +void Matdamageice::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols,int name, int type){ + /*Nothing updated yet*/ +} +/*}}}*/ +/*FUNCTION Matdamageice::InputUpdateFromVectorDakota(int* vector, int name, int type) {{{*/ +void Matdamageice::InputUpdateFromVectorDakota(int* vector, int name, int type){ + /*Nothing updated yet*/ +} +/*}}}*/ +/*FUNCTION Matdamageice::InputUpdateFromVectorDakota(bool* vector, int name, int type) {{{*/ +void Matdamageice::InputUpdateFromVectorDakota(bool* vector, int name, int type){ + /*Nothing updated yet*/ +} +/*}}}*/ +/*FUNCTION Matdamageice::InputUpdateFromConstant(IssmDouble constant, int name) {{{*/ +void Matdamageice::InputUpdateFromConstant(IssmDouble constant, int name){ + /*Nothing updated yet*/ +} +/*}}}*/ +/*FUNCTION Matdamageice::InputUpdateFromConstant(int constant, int name) {{{*/ +void Matdamageice::InputUpdateFromConstant(int constant, int name){ + /*Nothing updated yet*/ +} +/*}}}*/ +/*FUNCTION Matdamageice::InputUpdateFromConstant(bool constant, int name) {{{*/ +void Matdamageice::InputUpdateFromConstant(bool constant, int name){ + /*Nothing updated yet*/ +} +/*}}}*/ +/*FUNCTION Matdamageice::InputUpdateFromSolution{{{*/ +void Matdamageice::InputUpdateFromSolution(IssmDouble* solution){ + /*Nothing updated yet*/ +} +/*}}}*/ +/*FUNCTION Matdamageice::InputUpdateFromIoModel{{{*/ +void Matdamageice::InputUpdateFromIoModel(int index, IoModel* iomodel){ + + int i,j; + + int dim; + bool control_analysis; + int num_control_type; + + /*Fetch parameters: */ + iomodel->Constant(&dim,MeshDimensionEnum); + iomodel->Constant(&control_analysis,InversionIscontrolEnum); + if(control_analysis) iomodel->Constant(&num_control_type,InversionNumControlParametersEnum); + + /*if 2d*/ + if(dim==2){ + + /*Intermediaries*/ + const int num_vertices = 3; //Tria has 3 vertices + IssmDouble nodeinputs[num_vertices]; + IssmDouble cmmininputs[num_vertices]; + IssmDouble cmmaxinputs[num_vertices]; + + /*Get B*/ + if (iomodel->Data(MaterialsRheologyBEnum)) { + for(i=0;iData(MaterialsRheologyBEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; + this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,nodeinputs)); + } + + /*Get n*/ + if (iomodel->Data(MaterialsRheologyNEnum)) { + for(i=0;iData(MaterialsRheologyNEnum)[index]; + this->inputs->AddInput(new TriaP1Input(MaterialsRheologyNEnum,nodeinputs)); + } + + /*Get Z*/ + if (iomodel->Data(MaterialsRheologyZEnum)) { + for(i=0;iData(MaterialsRheologyZEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; + this->inputs->AddInput(new TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs)); + } + + /*Control Inputs*/ + #ifdef _HAVE_CONTROL_ + if (control_analysis && iomodel->Data(InversionControlParametersEnum)){ + for(i=0;i(iomodel->Data(InversionControlParametersEnum)[i])){ + case MaterialsRheologyBbarEnum: + if (iomodel->Data(MaterialsRheologyBEnum)){ + _assert_(iomodel->Data(MaterialsRheologyBEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); + for(j=0;jData(MaterialsRheologyBEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)]; + for(j=0;jData(InversionMinParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + case MaterialsRheologyZbarEnum: + if (iomodel->Data(MaterialsRheologyZEnum)){ + _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); + for(j=0;jData(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)]; + for(j=0;jData(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + + } + } + } + #endif + } + + /*if 3d*/ + #ifdef _HAVE_3D_ + else if(dim==3){ + + /*Intermediaries*/ + const int num_vertices = 6; //Penta has 6 vertices + IssmDouble nodeinputs[num_vertices]; + IssmDouble cmmininputs[num_vertices]; + IssmDouble cmmaxinputs[num_vertices]; + + /*Get B*/ + if (iomodel->Data(MaterialsRheologyBEnum)) { + for(i=0;iData(MaterialsRheologyBEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; + this->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,nodeinputs)); + } + + /*Get n*/ + if (iomodel->Data(MaterialsRheologyNEnum)) { + for(i=0;iData(MaterialsRheologyNEnum)[index]; + this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs)); + } + + /*Get Z*/ + if (iomodel->Data(MaterialsRheologyZEnum)) { + for(i=0;iData(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; + this->inputs->AddInput(new PentaP1Input(MaterialsRheologyZEnum,nodeinputs)); + } + + /*Control Inputs*/ + #ifdef _HAVE_CONTROL_ + if (control_analysis && iomodel->Data(InversionControlParametersEnum)){ + for(i=0;i(iomodel->Data(InversionControlParametersEnum)[i])){ + case MaterialsRheologyBbarEnum: + if (iomodel->Data(MaterialsRheologyBEnum)){ + _assert_(iomodel->Data(MaterialsRheologyBEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); + for(j=0;jData(MaterialsRheologyBEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)]; + for(j=0;jData(InversionMinParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + case MaterialsRheologyZbarEnum: + if (iomodel->Data(MaterialsRheologyZEnum)){ + _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); + for(j=0;jData(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)]; + for(j=0;jData(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + this->inputs->AddInput(new ControlInput(MaterialsRheologyZEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + } + } + } + #endif + } + #endif + else{ + _error_("Mesh type not supported yet!"); + } + + return; +} +/*}}}*/ +/*FUNCTION Matdamageice::IsInput{{{*/ +bool Matdamageice::IsInput(int name){ + if ( + name==MaterialsRheologyBEnum || + name==MaterialsRheologyBbarEnum || + name==MaterialsRheologyNEnum || + name==MaterialsRheologyZEnum || + name==MaterialsRheologyZbarEnum + ){ + return true; + } + else return false; +} +/*}}}*/ Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/classes/objects/Materials/Matdamageice.h =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/classes/objects/Materials/Matdamageice.h (revision 0) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/classes/objects/Materials/Matdamageice.h (revision 13130) @@ -0,0 +1,70 @@ +/*!\file Matdamageice.h + * \brief: header file for matice object + */ + +#ifndef MATDAMAGEICE_H_ +#define MATDAMAGEICE_H_ + +/*Headers:*/ +/*{{{*/ +#include "./Material.h" +class IoModel; +/*}}}*/ + +class Matdamageice: public Material{ + + private: + int mid; + Hook* helement; + + public: + /*Matdamageice constructors, destructors: {{{*/ + Matdamageice(); + Matdamageice(int mid,int i, IoModel* iomodel); + ~Matdamageice(); + /*}}}*/ + /*Object virtual functions definitions:{{{ */ + void Echo(); + void DeepEcho(); + int Id(); + int MyRank(); + int ObjectEnum(); + Object* copy(); + /*}}}*/ + /*Update virtual functions definitions: {{{*/ + void InputUpdateFromVector(IssmDouble* vector, int name, int type); + void InputUpdateFromVector(int* vector, int name, int type); + void InputUpdateFromVector(bool* vector, int name, int type); + void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrow, int ncols, int name, int type); + void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type); + void InputUpdateFromVectorDakota(int* vector, int name, int type); + void InputUpdateFromVectorDakota(bool* vector, int name, int type); + void InputUpdateFromConstant(IssmDouble constant, int name); + void InputUpdateFromConstant(int constant, int name); + void InputUpdateFromConstant(bool constant, int name); + void InputUpdateFromSolution(IssmDouble* solution); + void InputUpdateFromIoModel(int index, IoModel* iomodel); + /*}}}*/ + /*Material virtual functions resolution: {{{*/ + void InputDuplicate(int original_enum,int new_enum); + void Configure(Elements* elements); + void GetVectorFromInputs(Vector* vector,int input_enum); + void SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin); + void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon); + void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon); + void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon); + void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); + void GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); + void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); + void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); + IssmDouble GetA(); + IssmDouble GetB(); + IssmDouble GetBbar(); + IssmDouble GetN(); + IssmDouble GetZ(); + IssmDouble GetZbar(); + bool IsInput(int name); + /*}}}*/ +}; + +#endif /* _MATICE_H_ */