Index: /issm/trunk-jpl-damage/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Elements/Tria.cpp (revision 11374) +++ /issm/trunk-jpl-damage/src/c/objects/Elements/Tria.cpp (revision 11375) @@ -1440,5 +1440,5 @@ /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ - if (enum_type==MaterialsRheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type); + if (enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum) input=this->matice->inputs->GetInput(enum_type); else input=this->inputs->GetInput(enum_type); //if (!input) _error_("Input %s not found in tria->inputs",EnumToStringx(enum_type)); @@ -1544,4 +1544,5 @@ break; case MaterialsRheologyBbarEnum: + case MaterialsRheologyZbarEnum: /*Matice will take care of it*/ break; default: @@ -1732,5 +1733,5 @@ /*update input*/ - if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum){ + if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZbarEnum || name==MaterialsRheologyZEnum ){ matice->inputs->AddInput(new TriaP1Input(name,values)); } @@ -2732,4 +2733,7 @@ *presponse=this->matice->GetBbar(); break; + case MaterialsRheologyZbarEnum: + *presponse=this->matice->GetZbar(); + break; case VelEnum: @@ -3244,5 +3248,5 @@ for(int i=0;iinputs->GetInput(control_type[i]); _assert_(input); } @@ -3271,5 +3275,5 @@ Input* input=NULL; - if(enum_type==MaterialsRheologyBbarEnum){ + if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ input=(Input*)matice->inputs->GetInput(enum_type); } @@ -3289,5 +3293,5 @@ Input* input=NULL; - if(enum_type==MaterialsRheologyBbarEnum){ + if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ input=(Input*)matice->inputs->GetInput(enum_type); } @@ -3308,5 +3312,5 @@ Input* input=NULL; - if(enum_type==MaterialsRheologyBbarEnum){ + if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ input=(Input*)matice->inputs->GetInput(enum_type); } @@ -3338,4 +3342,7 @@ case MaterialsRheologyBbarEnum: GradjBMacAyeal(gradient); + break; + case MaterialsRheologyZbarEnum: + GradjZMacAyeal(gradient); break; case BalancethicknessThickeningRateEnum: @@ -3423,4 +3430,44 @@ } /*}}}*/ +/*FUNCTION Tria::GradjZGradient{{{1*/ +void Tria::GradjZGradient(Vec gradient, int weight_index){ + + int i,ig; + int doflist1[NUMVERTICES]; + double Jdet,weight; + double xyz_list[NUMVERTICES][3]; + double dbasis[NDOF2][NUMVERTICES]; + double dk[NDOF2]; + double grade_g[NUMVERTICES]={0.0}; + GaussTria *gauss=NULL; + + /*Retrieve all inputs we will be needing: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + GetDofList1(&doflist1[0]); + Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input); + Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); + weights_input->GetInputValue(&weight,gauss,weight_index); + + /*Build alpha_complement_list: */ + rheologyz_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for (i=0;iweight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]); + } + VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); + + /*Clean up and return*/ + delete gauss; +} +/*}}}*/ /*FUNCTION Tria::GradjBMacAyeal{{{1*/ void Tria::GradjBMacAyeal(Vec gradient){ @@ -3464,4 +3511,61 @@ this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(basis,gauss); + + /*standard gradient dJ/dki*/ + for (i=0;iweight*basis[i]; + } + + VecSetValues(gradient,NUMVERTICES,doflist,(const double*)grad,ADD_VALUES); + + /*clean-up*/ + delete gauss; +} +/*}}}*/ +/*FUNCTION Tria::GradjZMacAyeal{{{1*/ +void Tria::GradjZMacAyeal(Vec gradient){ + + /*Intermediaries*/ + int i,ig; + int doflist[NUMVERTICES]; + double vx,vy,lambda,mu,thickness,Jdet; + double viscosity_complement; + double dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dZ[NDOF2]; + double xyz_list[NUMVERTICES][3]; + double basis[3],epsilon[3]; + double grad[NUMVERTICES]={0.0}; + GaussTria *gauss = NULL; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + GetDofList1(&doflist[0]); + + /*Retrieve all inputs*/ + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(4); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + thickness_input->GetInputValue(&thickness,gauss); + rheologyz_input->GetInputDerivativeValue(&dZ[0],&xyz_list[0][0],gauss); + vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); + vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); + adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss); + adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss); + + this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + matice->GetViscosityZComplement(&viscosity_complement,&epsilon[0]); GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); Index: /issm/trunk-jpl-damage/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Elements/Tria.h (revision 11374) +++ /issm/trunk-jpl-damage/src/c/objects/Elements/Tria.h (revision 11375) @@ -142,5 +142,7 @@ void Gradj(Vec gradient,int control_type); void GradjBGradient(Vec gradient,int weight_index); + void GradjZGradient(Vec gradient,int weight_index); void GradjBMacAyeal(Vec gradient); + void GradjZMacAyeal(Vec gradient); void GradjDragMacAyeal(Vec gradient); void GradjDragStokes(Vec gradient);