Changeset 15453 for issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

Timestamp:

07/07/13 13:58:41 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: preparing Newton for P2 elements

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (9 diffs)

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -3119 +3119 @@
         /*Initialize Element matrix, vectors and Gaussian points*/
         ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal(); //Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)
-        IssmDouble*    dphi   = xNew<IssmDouble>(2*numnodes);
+        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
 
         /*Retrieve all inputs and parameters*/
@@ -3134 +3134 @@
 
                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-                GetNodalFunctionsDerivatives(dphi,&xyz_list[0][0],gauss);
+                GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
 
                 thickness_input->GetInputValue(&thickness, gauss);
@@ -3144 +3144 @@
                 for(i=0;i<numnodes;i++){
                         for(j=0;j<numnodes;j++){
-                                eps1dotdphii=eps1[0]*dphi[0*numnodes+i]+eps1[1]*dphi[1*numnodes+i];
-                                eps1dotdphij=eps1[0]*dphi[0*numnodes+j]+eps1[1]*dphi[1*numnodes+j];
-                                eps2dotdphii=eps2[0]*dphi[0*numnodes+i]+eps2[1]*dphi[1*numnodes+i];
-                                eps2dotdphij=eps2[0]*dphi[0*numnodes+j]+eps2[1]*dphi[1*numnodes+j];
+                                eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i];
+                                eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j];
+                                eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i];
+                                eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j];
 
                                 Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
@@ -3158 +3158 @@
 
         /*Transform Coordinate System*/
-        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
+        TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum);
 
         /*Clean up and return*/
+        xDelete<IssmDouble>(dbasis);
         delete gauss;
-        xDelete<IssmDouble>(dphi);
         return Ke;
 }
@@ -5214 +5214 @@
 
         /*Constants*/
-        const int    numdof=NDOF2*NUMVERTICES;
+        const int numnodes = this->GetNumberOfNodes();
+        const int numdof   = NDOF2*numnodes;
 
         /*Intermediaries */
-        int        i,j;
-        bool       incomplete_adjoint;
-        IssmDouble xyz_list[NUMVERTICES][3];
-        IssmDouble Jdet,thickness;
-        IssmDouble eps1dotdphii,eps1dotdphij;
-        IssmDouble eps2dotdphii,eps2dotdphij;
-        IssmDouble mu_prime;
-        IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];*/
-        IssmDouble eps1[2],eps2[2];
-        IssmDouble dphi[2][NUMVERTICES];
-        GaussTria *gauss=NULL;
+        int         i,j;
+        bool        incomplete_adjoint;
+        IssmDouble  xyz_list[NUMVERTICES][3];
+        IssmDouble  Jdet,thickness;
+        IssmDouble  eps1dotdphii,eps1dotdphij;
+        IssmDouble  eps2dotdphii,eps2dotdphij;
+        IssmDouble  mu_prime;
+        IssmDouble  epsilon[3];/* epsilon=[exx,eyy,exy];*/
+        IssmDouble  eps1[2],eps2[2];
+        GaussTria  *gauss=NULL;
 
         /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/
@@ -5240 +5240 @@
         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
 
+        /*Allocate dbasis*/
+        IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
+
         /* Start  looping on the number of gaussian points: */
         gauss=new GaussTria(2);
@@ -5247 +5250 @@
 
                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-                GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss);
+                GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
 
                 thickness_input->GetInputValue(&thickness, gauss);
@@ -5255 +5258 @@
                 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;
+                for(i=0;i<numnodes;i++){
+                        for(j=0;j<numnodes;j++){
+                                eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i];
+                                eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j];
+                                eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i];
+                                eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j];
+
+                                Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
+                                Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
+                                Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
+                                Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
                         }
                 }
@@ -5271 +5274 @@
 
         /*Transform Coordinate System*/
-        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
+        TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum);
 
         /*Clean up and return*/
         delete gauss;
+        xDelete<IssmDouble>(dbasis);
         //Ke->Transpose();
         return Ke;