Changeset 5571

Timestamp:

08/25/10 10:03:20 (15 years ago)

Author:

Mathieu Morlighem

Message:

Better gradient for balanced thickness

File:

• issm/trunk/src/c/objects/Elements/Tria.cpp (modified) (1 diff)

issm/trunk/src/c/objects/Elements/Tria.cpp

@@ -1342 +1342 @@
         const int    numgrids=3;
         const int    NDOF1=1;
-        const int    numdof=NDOF1*numgrids;
-        double       xyz_list[numgrids][3];
         int          doflist1[numgrids];
 
+        /*Gauss*/
+        double  gauss[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
+
         /* grid data: */
-        double DhDt[numgrids];
-
-        /* gaussian points: */
-        int     num_gauss,ig;
-        double* first_gauss_area_coord  =  NULL;
-        double* second_gauss_area_coord =  NULL;
-        double* third_gauss_area_coord  =  NULL;
-        double* gauss_weights           =  NULL;
-        double  gauss_weight;
-        double  gauss_l1l2l3[3];
+        double lambda;
 
         /*element vector at the gaussian points: */
-        double  grade_g[numgrids]={0.0};
-        double  grade_g_gaussian[numgrids];
-
-        /* Jacobian: */
-        double Jdet;
-
-        /*nodal functions: */
-        double l1l2l3[3];
-        double  lambda;
+        double  gradient_g[numgrids];
 
         /*Inputs*/
         Input* adjoint_input=NULL;
 
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+        /*Retrieve dof list*/
         GetDofList1(&doflist1[0]);
 
-        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
-
         /*Retrieve all inputs we will be needing: */
-        adjoint_input=inputs->GetInput(AdjointEnum);
-
-        /* Start  looping on the number of gaussian points: */
-        for (ig=0; ig<num_gauss; ig++){
-                /*Pick up the gaussian point: */
-                gauss_weight=*(gauss_weights+ig);
-                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
-                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
-                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
-
-                /*Get thickness: */
-                adjoint_input->GetParameterValue(&lambda, gauss_l1l2l3);
-
-                /* Get Jacobian determinant: */
-                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-
-                /* Get nodal functions value at gaussian point:*/
-                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
-
-                /*DhDtuild gradje_g_gaussian vector (actually -dJ/dDhDt): */
-                for (i=0;i<numgrids;i++){
-                        //standard gradient dJ/dki
-                        grade_g_gaussian[i]=-lambda*Jdet*gauss_weight*l1l2l3[i]; 
-                }
-
-                /*Add grade_g_gaussian to grade_g: */
-                for( i=0; i<numgrids;i++) grade_g[i]+=grade_g_gaussian[i];
+        adjoint_input=(Input*)inputs->GetInput(AdjointEnum);
+
+        /* Start  looping on the vertices: */
+        for(i=0; i<numgrids;i++){
+                adjoint_input->GetParameterValue(&lambda,&gauss[i][0]);
+                gradient_g[i]=-lambda;
         }
 
         /*Add grade_g to global vector gradient: */
-        VecSetValues(gradient,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
-
-        xfree((void**)&first_gauss_area_coord);
-        xfree((void**)&second_gauss_area_coord);
-        xfree((void**)&third_gauss_area_coord);
-        xfree((void**)&gauss_weights);
+        VecSetValues(gradient,numgrids,doflist1,(const double*)gradient_g,INSERT_VALUES);
 }
 /*}}}*/