source: issm/oecreview/Archive/11431-11450/ISSM-11447-11448.diff@ 11515

proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Tria.cpp

@@ -335 +335 @@
         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
         switch(analysis_type){
                 #ifdef _HAVE_DIAGNOSTIC_
-                case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum:
+                case DiagnosticHorizAnalysisEnum:
                         Ke=CreateKMatrixDiagnosticMacAyeal();
                         break;
+                case AdjointHorizAnalysisEnum:
+                        Ke=CreateKMatrixAdjointMacAyeal();
+                        break;
                 case DiagnosticHutterAnalysisEnum:
                         Ke=CreateKMatrixDiagnosticHutter();
                         break;
@@ -3062 +3065 @@
         return pe;
 }
 /*}}}*/
-/*FUNCTION Tria::CreateJacobianDiagnosticPattyn{{{1*/
+/*FUNCTION Tria::CreateJacobianDiagnosticMacayeal{{{1*/
 ElementMatrix* Tria::CreateJacobianDiagnosticMacayeal(void){
 
         /*Constants*/
@@ -3863 +3866 @@
                 rheologyb_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
 
                 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
-                //Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
+                Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
         }
 
         /*Clean up and return*/
@@ -4735 +4738 @@
                 drag_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
 
                 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
-                //Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
+                Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
         }
 
         /*Clean up and return*/
@@ -4765 +4768 @@
         return Ke;
 }
 /*}}}*/
+/*FUNCTION Tria::CreateKMatrixAdjointMacAyeal{{{1*/
+ElementMatrix* Tria::CreateKMatrixAdjointMacAyeal(void){
+
+        /*Constants*/
+        const int    numdof=NDOF2*NUMVERTICES;
+
+        /*Intermediaries */
+        int        i,j,ig;
+        bool       incomplete_adjoint;
+        double     xyz_list[NUMVERTICES][3];
+        double     Jdet,thickness;
+        double     eps1dotdphii,eps1dotdphij;
+        double     eps2dotdphii,eps2dotdphij;
+        double     mu_prime;
+        double     epsilon[3];/* epsilon=[exx,eyy,exy];*/
+        double     eps1[2],eps2[2];
+        double     phi[NUMVERTICES];
+        double     dphi[2][NUMVERTICES];
+        GaussTria *gauss=NULL;
+
+        /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/
+        parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
+        ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal();
+        if(incomplete_adjoint) return Ke;
+
+        /*Retrieve all inputs and parameters*/
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
+        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
+        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_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(&dphi[0][0],&xyz_list[0][0],gauss);
+
+                thickness_input->GetInputValue(&thickness, gauss);
+                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+                matice->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);
+                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
+                eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
+
+                for(i=0;i<3;i++){
+                        for(j=0;j<3;j++){
+                                eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i];
+                                eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j];
+                                eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i];
+                                eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j];
+
+                                Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
+                                Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
+                                Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
+                                Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
+                        }
+                }
+        }
+
+        /*Transform Coordinate System*/
+        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
+
+        /*Clean up and return*/
+        delete gauss;
+        //Ke->Transpose();
+        return Ke;
+}
+/*}}}*/
 /*FUNCTION Tria::InputUpdateFromSolutionAdjointHoriz {{{1*/
 void  Tria::InputUpdateFromSolutionAdjointHoriz(double* solution){
 

proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Tria.h

@@ -216 +216 @@
 
                 #ifdef _HAVE_CONTROL_
                 ElementMatrix* CreateKMatrixAdjointBalancethickness(void);
+                ElementMatrix* CreateKMatrixAdjointMacAyeal(void);
                 ElementVector* CreatePVectorAdjointHoriz(void);
                 ElementVector* CreatePVectorAdjointStokes(void);
                 ElementVector* CreatePVectorAdjointBalancethickness(void);

proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Penta.cpp

@@ -584 +584 @@
         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
         switch(analysis_type){
                 #ifdef _HAVE_DIAGNOSTIC_
-                case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum:
+                case DiagnosticHorizAnalysisEnum:
                         Ke=CreateKMatrixDiagnosticHoriz(); De=CreateDVectorDiagnosticHoriz();
                         break;
+                case AdjointHorizAnalysisEnum:
+                        Ke=CreateKMatrixAdjointHoriz();
+                        break;
                 case DiagnosticHutterAnalysisEnum:
                         Ke=CreateKMatrixDiagnosticHutter();
                         break;
@@ -4442 +4445 @@
         ((ControlInput*)input)->SetGradient(grad_input);
 
 }/*}}}*/
+/*FUNCTION Penta::CreateKMatrixAdjointHoriz{{{1*/
+ElementMatrix* Penta::CreateKMatrixAdjointHoriz(void){
+
+        int approximation;
+        inputs->GetInputValue(&approximation,ApproximationEnum);
+
+        switch(approximation){
+                case MacAyealApproximationEnum:
+                        return CreateKMatrixAdjointMacAyeal2d();
+                case PattynApproximationEnum:
+                        return CreateKMatrixAdjointPattyn();
+                case StokesApproximationEnum:
+                        return CreateKMatrixAdjointStokes();
+                case NoneApproximationEnum:
+                        return NULL;
+                default:
+                        _error_("Approximation %s not supported yet",EnumToStringx(approximation));
+        }
+}
+/*}}}*/
+/*FUNCTION Penta::CreateKMatrixAdjointMacAyeal2d{{{1*/
+ElementMatrix* Penta::CreateKMatrixAdjointMacAyeal2d(void){
+
+        /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the 
+          bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build 
+          the stiffness matrix. */
+        if (!IsOnBed()) return NULL;
+
+        /*Depth Averaging B*/
+        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
+
+        /*Call Tria function*/
+        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
+        ElementMatrix* Ke=tria->CreateKMatrixAdjointMacAyeal();
+        delete tria->matice; delete tria;
+
+        /*Delete B averaged*/
+        this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum);
+
+        /*clean up and return*/
+        return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateKMatrixAdjointPattyn{{{1*/
+ElementMatrix* Penta::CreateKMatrixAdjointPattyn(void){
+
+        /*Constants*/
+        const int    numdof=NDOF2*NUMVERTICES;
+
+        /*Intermediaries */
+        int        i,j,ig;
+        bool       incomplete_adjoint;
+        double     xyz_list[NUMVERTICES][3];
+        double     Jdet;
+        double     eps1dotdphii,eps1dotdphij;
+        double     eps2dotdphii,eps2dotdphij;
+        double     mu_prime;
+        double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+        double     eps1[3],eps2[3];
+        double     phi[NUMVERTICES];
+        double     dphi[3][NUMVERTICES];
+        GaussPenta *gauss=NULL;
+
+        /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
+        parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
+        ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
+        if(incomplete_adjoint) return Ke;
+
+        /*Retrieve all inputs and parameters*/
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
+        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
+
+        /* Start  looping on the number of gaussian points: */
+        gauss=new GaussPenta(5,5);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
+
+                this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+                matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
+                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
+                eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
+                eps1[2]=epsilon[3];                eps2[2]=epsilon[4];
+
+                for(i=0;i<6;i++){
+                        for(j=0;j<6;j++){
+                                eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
+                                eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
+                                eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
+                                eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
+
+                                Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
+                                Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
+                                Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
+                                Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
+                        }
+                }
+        }
+
+        /*Transform Coordinate System*/
+        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
+
+        /*Clean up and return*/
+        delete gauss;
+        return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateKMatrixAdjointStokes{{{1*/
+ElementMatrix* Penta::CreateKMatrixAdjointStokes(void){
+
+        /*Constants*/
+        const int    numdof=NDOF4*NUMVERTICES;
+
+        /*Intermediaries */
+        int        i,j,ig;
+        bool       incomplete_adjoint;
+        double     xyz_list[NUMVERTICES][3];
+        double     Jdet;
+        double     eps1dotdphii,eps1dotdphij;
+        double     eps2dotdphii,eps2dotdphij;
+        double     eps3dotdphii,eps3dotdphij;
+        double     mu_prime;
+        double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+        double     eps1[3],eps2[3],eps3[3];
+        double     phi[NUMVERTICES];
+        double     dphi[3][NUMVERTICES];
+        GaussPenta *gauss=NULL;
+
+        /*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/
+        parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
+        ElementMatrix* Ke=CreateKMatrixDiagnosticStokes();
+        if(incomplete_adjoint) return Ke;
+
+        /*Retrieve all inputs and parameters*/
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
+        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
+        Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
+
+        /* Start  looping on the number of gaussian points: */
+        gauss=new GaussPenta(5,5);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
+
+                this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+                matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
+                eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
+                eps1[1]=epsilon[2];   eps2[1]=epsilon[1];   eps3[1]=epsilon[4];
+                eps1[2]=epsilon[3];   eps2[2]=epsilon[4];   eps3[2]= -epsilon[0] -epsilon[1];
+
+                for(i=0;i<6;i++){
+                        for(j=0;j<6;j++){
+                                eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
+                                eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
+                                eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
+                                eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
+                                eps3dotdphii=eps3[0]*dphi[0][i]+eps3[1]*dphi[1][i]+eps3[2]*dphi[2][i];
+                                eps3dotdphij=eps3[0]*dphi[0][j]+eps3[1]*dphi[1][j]+eps3[2]*dphi[2][j];
+
+                                Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
+                                Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
+                                Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
+
+                                Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
+                                Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
+                                Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
+
+                                Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
+                                Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
+                                Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
+                        }
+                }
+        }
+
+        /*Transform Coordinate System*/
+        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
+
+        /*Clean up and return*/
+        delete gauss;
+        return Ke;
+}
+/*}}}*/
 /*FUNCTION Penta::CreatePVectorAdjointHoriz{{{1*/
 ElementVector* Penta::CreatePVectorAdjointHoriz(void){
 

proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Penta.h

@@ -215 +215 @@
                 ElementMatrix* CreateKMatrixCouplingMacAyealStokesFriction(void);
                 ElementMatrix* CreateKMatrixCouplingPattynStokes(void);
                 ElementMatrix* CreateKMatrixDiagnosticHoriz(void);
+                ElementMatrix* CreateKMatrixAdjointHoriz(void);
                 ElementVector* CreateDVectorDiagnosticHoriz(void);
                 ElementVector* CreateDVectorDiagnosticStokes(void);
                 ElementMatrix* CreateKMatrixDiagnosticHutter(void);
@@ -274 +275 @@
 
                 #ifdef _HAVE_CONTROL_
                 ElementVector* CreatePVectorAdjointHoriz(void);
+                ElementMatrix* CreateKMatrixAdjointMacAyeal2d(void);
+                ElementMatrix* CreateKMatrixAdjointPattyn(void);
+                ElementMatrix* CreateKMatrixAdjointStokes(void);
                 ElementVector* CreatePVectorAdjointMacAyeal(void);
                 ElementVector* CreatePVectorAdjointPattyn(void);
                 ElementVector* CreatePVectorAdjointStokes(void);