Changeset 17370

Timestamp:

03/03/14 09:48:20 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: trying to fix AD

File:

• issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (modified) (7 diffs)

issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

@@ -93 +93 @@
 
         /*Intermediaries */
-        const int  dim = 2; // solve for LSF in horizontal plane only
+        int  dim = 2; // solve for LSF in horizontal plane only
         int i, row, col;
         IssmDouble kappa;
         IssmDouble Jdet, dt, D_scalar;
         IssmDouble h,hx,hy,hz;
-        IssmDouble vel,v[dim],w[dim];
-        IssmDouble calvingrate, c[dim];
+        IssmDouble vel;
+        IssmDouble calvingrate;
         IssmDouble dlsf[dim], norm_dlsf, normal[dim];
         IssmDouble* xyz_list = NULL;
@@ -111 +111 @@
         IssmDouble*    B        = xNew<IssmDouble>(dim*numnodes);
         IssmDouble*    Bprime   = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble     D[dim][dim];
+        IssmDouble*    D        = xNew<IssmDouble>(dim*dim);
+        IssmDouble*    v        = xNew<IssmDouble>(dim);
+        IssmDouble*    w        = xNew<IssmDouble>(dim);
+        IssmDouble*    c        = xNew<IssmDouble>(dim);
 
         /*Retrieve all inputs and parameters*/
@@ -165 +168 @@
                         for(col=0;col<dim;col++)
                                 if(row==col)
-                                        D[row][col]=D_scalar*w[row];
+                                        D[row*dim+col]=D_scalar*w[row];
                                 else
-                                   D[row][col]=0.;
+                                   D[row*dim+col]=0.;
 
                 TripleMultiply(B,dim,numnodes,1,
-                                        &D[0][0],dim,dim,0,
+                                        D,dim,dim,0,
                                         Bprime,dim,numnodes,0,
                                         &Ke->values[0],1);
@@ -177 +180 @@
                 int stabilization=2;
                 vel=0.;
-                for(i=0;i<dim;i++) vel+=pow(v[i],2);
+                for(i=0;i<dim;i++) vel+=v[i]*v[i];
                 vel=sqrt(vel)+1.e-14;
                 switch(stabilization){
@@ -191 +194 @@
                                         for(col=0;col<dim;col++)
                                         if(row==col)
-                                                D[row][col]=D_scalar*kappa;
+                                                D[row*dim+col]=D_scalar*kappa;
                                         else
-                                                D[row][col]=0.;
+                                                D[row*dim+col]=0.;
 
                                 TripleMultiply(Bprime,dim,numnodes,1,
-                                                        &D[0][0],dim,dim,0,
+                                                        D,dim,dim,0,
                                                         Bprime,dim,numnodes,0,
                                                         &Ke->values[0],1);
@@ -206 +209 @@
                                 for(row=0;row<dim;row++) 
                                         for(col=0;col<dim;col++) 
-                                                D[row][col] = D_scalar*h/(2.*vel)*v[row]*v[col];
+                                                D[row*dim+col] = D_scalar*h/(2.*vel)*v[row]*v[col];
 
                                 TripleMultiply(Bprime,dim,numnodes,1,
-                                                        &D[0][0],dim,dim,0,
+                                                        D,dim,dim,0,
                                                         Bprime,dim,numnodes,0,
                                                         &Ke->values[0],1);
@@ -222 +225 @@
         xDelete<IssmDouble>(basis);
         xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(D);
         xDelete<IssmDouble>(Bprime);
+        xDelete<IssmDouble>(v);
+        xDelete<IssmDouble>(w);
+        xDelete<IssmDouble>(c);
         delete gauss;
         return Ke;