[cig-commits] r6568 - in long/3D/Gale/trunk: . src/Gale/plugins/SurfaceProcess

[walter at geodynamics.org]

Author: walter
Date: 2007-04-13 15:47:36 -0700 (Fri, 13 Apr 2007)
New Revision: 6568

Modified:
   long/3D/Gale/trunk/
   long/3D/Gale/trunk/src/Gale/plugins/SurfaceProcess/SurfaceProcess.c
Log:
 r1801 at earth:  boo | 2007-04-13 15:45:56 -0700
 Add in a 3D SurfaceProcess example code



Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
   - 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1800
   + 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1801

Modified: long/3D/Gale/trunk/src/Gale/plugins/SurfaceProcess/SurfaceProcess.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/plugins/SurfaceProcess/SurfaceProcess.c	2007-04-13 22:47:34 UTC (rev 6567)
+++ long/3D/Gale/trunk/src/Gale/plugins/SurfaceProcess/SurfaceProcess.c	2007-04-13 22:47:36 UTC (rev 6568)
@@ -96,6 +96,41 @@
             }
             free(newVal);
           }
+          else if( surface->topo->nDims == 2 ) {
+            int n_i, n, nx, ny;
+            double *newVal;
+            CartesianGenerator *gen;
+
+            n=Mesh_GetLocalSize(surface,MT_VERTEX);
+            gen=((CartesianGenerator *)(surface->generator));
+            nx=gen->range[0]+1;
+            ny=gen->range[1]+1;
+
+            newVal=malloc(sizeof(double)*n);
+            
+            for( n_i = 0; n_i < n; n_i++ ) {
+              newVal[n_i]=0;
+              /* Make sure we are not on the boundary */
+              if( n_i%nx!=0 && (n_i+1)%nx!=0  /* x boundary */
+                  && n_i>nx && n_i<n-nx) {    /* y boundary */
+                double dx = surface->verts[n_i][0]
+                  - surface->verts[n_i - 1][0];
+                double dy = surface->verts[n_i][1]
+                  - surface->verts[n_i - nx][1];
+                
+                newVal[n_i] = 
+                  (heights[n_i - 1] - 2*heights[n_i] + heights[n_i + 1])
+                  / (dx * dx)
+                  + (heights[n_i-nx] - 2*heights[n_i] + heights[n_i+nx])
+                  / (dy * dy);
+              }
+            }
+            for( n_i = 0; n_i < nx; n_i++ ) {
+              heights[n_i] = heights[n_i] + newVal[n_i]*dt;
+            }
+            free(newVal);
+            
+          }
           else {
             abort();
           }