Changeset 11338

Timestamp:

02/06/12 16:59:01 (13 years ago)

Author:

Mathieu Morlighem

Message:

Prepared Stokes Newton's iterations

Location:

issm/trunk-jpl/src/c

Files:

• objects/Elements/Penta.cpp (modified) (2 diffs)
• objects/Elements/Penta.h (modified) (1 diff)
• solvers/solver_newton.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/objects/Elements/Penta.cpp

@@ -788 +788 @@
 #ifdef _HAVE_DIAGNOSTIC_
                 case DiagnosticHorizAnalysisEnum:
-                        Ke=CreateJacobianDiagnosticPattyn();
+                        Ke=CreateJacobianDiagnosticHoriz();
                         break;
 #endif
@@ -7424 +7424 @@
 }
 /*}}}*/
+/*FUNCTION Penta::CreateJacobianDiagnosticHoriz{{{1*/
+ElementMatrix* Penta::CreateJacobianDiagnosticHoriz(void){
+
+        int approximation;
+        inputs->GetInputValue(&approximation,ApproximationEnum);
+
+        switch(approximation){
+                case MacAyealApproximationEnum:
+                        _error_("Not supported yet");
+                case PattynApproximationEnum:
+                        return CreateJacobianDiagnosticPattyn();
+                case StokesApproximationEnum:
+                        return CreateJacobianDiagnosticStokes();
+                default:
+                        _error_("Approximation %s not supported yet",EnumToStringx(approximation));
+        }
+}
+/*}}}*/
 /*FUNCTION Penta::CreateJacobianDiagnosticPattyn{{{1*/
 ElementMatrix* Penta::CreateJacobianDiagnosticPattyn(void){
 
+        /*Constants*/
+        const int    numdof=NDOF2*NUMVERTICES;
+
+        /*Intermediaries */
+        int        i,j,ig;
+        int        approximation;
+        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)*/
+        ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
+
+        /*Retrieve all inputs and parameters*/
+        inputs->GetInputValue(&approximation,ApproximationEnum);
+        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::CreateJacobianDiagnosticStokes{{{1*/
+ElementMatrix* Penta::CreateJacobianDiagnosticStokes(void){
+
+        _error_("Not implemented yet");
         /*Constants*/
         const int    numdof=NDOF2*NUMVERTICES;

issm/trunk-jpl/src/c/objects/Elements/Penta.h

@@ -235 +235 @@
                 ElementMatrix* CreateKMatrixDiagnosticVertVolume(void);
                 ElementMatrix* CreateKMatrixDiagnosticVertSurface(void);
+                ElementMatrix* CreateJacobianDiagnosticHoriz(void);
                 ElementMatrix* CreateJacobianDiagnosticPattyn(void);
+                ElementMatrix* CreateJacobianDiagnosticStokes(void);
                 void           InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
                 void           InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong);

issm/trunk-jpl/src/c/solvers/solver_newton.cpp

@@ -78 +78 @@
 
                 CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kmax);
-                Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); MatFree(&Jff);//VecFree(&pJf);
+                Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); MatFree(&Jff);VecFree(&pJf);
                 VecAXPY(uf,1.,duf);      VecFree(&duf);
                 Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);