Changeset 11375

Timestamp:

02/08/12 17:17:19 (13 years ago)

Author:

cborstad

Message:

tria changes for new damage variable

Location:

issm/trunk-jpl-damage/src/c/objects/Elements

Files:

• Tria.cpp (modified) (11 diffs)
• Tria.h (modified) (1 diff)

issm/trunk-jpl-damage/src/c/objects/Elements/Tria.cpp

@@ -1440 +1440 @@
 
         /*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 +1544 @@
                                         break;
                                 case MaterialsRheologyBbarEnum:
+                                case MaterialsRheologyZbarEnum:
                                         /*Matice will take care of it*/ break;
                                 default:
@@ -1732 +1733 @@
 
                         /*update input*/
-                        if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum){
+                        if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZbarEnum || name==MaterialsRheologyZEnum ){
                                 matice->inputs->AddInput(new TriaP1Input(name,values));
                         }
@@ -2732 +2733 @@
                         *presponse=this->matice->GetBbar();
                         break;
+                case MaterialsRheologyZbarEnum:
+                        *presponse=this->matice->GetZbar();
+                        break;
                 case VelEnum:
 
@@ -3244 +3248 @@
         for(int i=0;i<num_controls;i++){
 
-                if(control_type[i]==MaterialsRheologyBbarEnum){
+                if(control_type[i]==MaterialsRheologyBbarEnum || control_type[i]==MaterialsRheologyZbarEnum){
                         input=(Input*)matice->inputs->GetInput(control_type[i]); _assert_(input);
                 }
@@ -3271 +3275 @@
         Input* input=NULL;
 
-        if(enum_type==MaterialsRheologyBbarEnum){
+        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
                 input=(Input*)matice->inputs->GetInput(enum_type);
         }
@@ -3289 +3293 @@
         Input* input=NULL;
 
-        if(enum_type==MaterialsRheologyBbarEnum){
+        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
                 input=(Input*)matice->inputs->GetInput(enum_type);
         }
@@ -3308 +3312 @@
         Input* input=NULL;
 
-        if(enum_type==MaterialsRheologyBbarEnum){
+        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
                 input=(Input*)matice->inputs->GetInput(enum_type);
         }
@@ -3338 +3342 @@
                 case MaterialsRheologyBbarEnum:
                         GradjBMacAyeal(gradient);
+                        break;
+                case MaterialsRheologyZbarEnum:
+                        GradjZMacAyeal(gradient);
                         break;
                 case BalancethicknessThickeningRateEnum:
@@ -3423 +3430 @@
 }
 /*}}}*/
+/*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();ig<gauss->end();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;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(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 +3511 @@
                 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;i<NUMVERTICES;i++) grad[i]+=-viscosity_complement*thickness*(
+                                        (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
+                                        )*Jdet*gauss->weight*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();ig<gauss->end();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);

issm/trunk-jpl-damage/src/c/objects/Elements/Tria.h

@@ -142 +142 @@
                 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);