Changeset 16422

Timestamp:

10/16/13 08:20:15 (11 years ago)

Author:

Mathieu Morlighem

Message:

BUG: fixing InputUpdateFromSolutionFS bug, in which we do not transform the solution vector correctly

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• Penta.cpp (modified) (6 diffs)
• Tria.cpp (modified) (1 diff)

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

@@ -10693 +10693 @@
 
         /*Initialize values*/
-        IssmDouble* vvalues  = xNew<IssmDouble>(vnumdof);
-        IssmDouble* pvalues  = xNew<IssmDouble>(pnumdof);
+        IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
         IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
         IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
@@ -10701 +10700 @@
         IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
 
+        /*Prepare coordinate system list*/
+        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
+        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
         /*Get dof list: */
         GetDofListVelocity(&vdoflist,GsetEnum);
@@ -10706 +10710 @@
 
         /*Use the dof list to index into the solution vector: */
-        for(i=0;i<vnumdof;i++) vvalues[i]=solution[vdoflist[i]];
-        for(i=0;i<pnumdof;i++) pvalues[i]=solution[pdoflist[i]];
+        for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
+        for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
 
         /*Transform solution in Cartesian Space*/
-        TransformSolutionCoord(&vvalues[0],nodes,vnumnodes,XYZEnum);
+        TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
 
         /*Ok, we have vx and vy in values, fill in all arrays: */
         for(i=0;i<vnumnodes;i++){
-                vx[i] = vvalues[i*NDOF3+0];
-                vy[i] = vvalues[i*NDOF3+1];
-                vz[i] = vvalues[i*NDOF3+2];
+                vx[i] = values[i*NDOF3+0];
+                vy[i] = values[i*NDOF3+1];
+                vz[i] = values[i*NDOF3+2];
                 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
@@ -10722 +10726 @@
         }
         for(i=0;i<pnumnodes;i++){
-                pressure[i] = pvalues[i];
+                pressure[i] = values[vnumdof+i];
                 if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
         }
@@ -10728 +10732 @@
         /*Recondition pressure and compute vel: */
         this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-        for(i=0;i<pnumnodes;i++) pressure[i]=pressure[i]*FSreconditioning;
-        for(i=0;i<vnumnodes;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+        for(i = 0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
+        for(i = 0;i<vnumnodes;i++) vel[i]      = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
 
         /*Now, we have to move the previous inputs  to old 
@@ -10751 +10755 @@
         xDelete<IssmDouble>(vy);
         xDelete<IssmDouble>(vx);
-        xDelete<IssmDouble>(vvalues);
-        xDelete<IssmDouble>(pvalues);
+        xDelete<IssmDouble>(values);
         xDelete<int>(vdoflist);
         xDelete<int>(pdoflist);
+        xDelete<int>(cs_list);
 }
 /*}}}*/

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

@@ -4417 +4417 @@
         GetDofListPressure(&pdoflist,GsetEnum);
 
+        /*Use the dof list to index into the solution vector: */
         for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
         for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];