[cig-commits] r13926 - in long/3D/Gale/trunk: . src/Underworld/plugins/EulerDeform

[walter at geodynamics.org]

Author: walter
Date: 2009-01-22 15:28:12 -0800 (Thu, 22 Jan 2009)
New Revision: 13926

Modified:
   long/3D/Gale/trunk/
   long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c
Log:
 r2467 at dante:  boo | 2009-01-22 15:26:37 -0800
 Generalize the corner interpolation



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

Modified: long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c	2009-01-22 23:28:10 UTC (rev 13925)
+++ long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c	2009-01-22 23:28:12 UTC (rev 13926)
@@ -671,12 +671,16 @@
 	return True;
 }
 
+Bool _EulerDeform_LineInterp(const double** crds,const double* pnt,
+                             unsigned fromDim,unsigned toDim, 
+			     double* val);
+Bool _EulerDeform_QuadYInterp(double** crds,const double* pnt,double* val);
+
 /* Remesh the left or right top corners in 2D or 3D */
 
 void EulerDeform_Remesh_Corner(Mesh *mesh, int corner, int inside,
                                double side_coord) {
   IJK		ijk;
-  unsigned n_corner, n_interior, n, n_in, max_z, z;
   Grid *grid;
   grid =
     *(Grid**)ExtensionManager_Get( mesh->info, mesh, 
@@ -684,12 +688,9 @@
                                                                "vertexGrid" ) );
   ijk[0]=corner;
   ijk[1]=grid->sizes[1]-1;
-  max_z=1;
-  if(Mesh_GetDimSize(mesh)==3)
-    max_z=grid->sizes[2];
-  
-  for(ijk[2]=0;ijk[2]<max_z;++ijk[2])
+  if(Mesh_GetDimSize(mesh)==2)
     {
+      unsigned n_corner, n_interior, n, n_in;
       n_corner=RegularMeshUtils_Node_3DTo1D(mesh,ijk);
       ijk[0]=inside;
       n_interior=RegularMeshUtils_Node_3DTo1D(mesh,ijk);
@@ -698,25 +699,72 @@
          && Mesh_GlobalToDomain(mesh,MT_VERTEX,
                                 n_interior,&n_in))
         {
-          double x0,y0,x1,y1;
-          double delta;
-          
-          x0=mesh->verts[n][0];
-          y0=mesh->verts[n][1];
-          x1=mesh->verts[n_in][0];
-          y1=mesh->verts[n_in][1];
-          
-          if((x0-side_coord)*(corner-inside)<0)
+          double *crds[2];
+
+          crds[0]=mesh->verts[n];
+          crds[1]=mesh->verts[n_in];
+          if(!_EulerDeform_LineInterp((const double**)crds,&side_coord,0,1,
+                                       &(mesh->verts[n][1])))
             {
-              printf("The side is moving in the wrong direction: %g %g\n",
-                     x0,side_coord);
+              printf("The side is moving in the wrong direction.\n");
               abort();
             }
           mesh->verts[n][0]=side_coord;
-          delta=(x1-side_coord)/(x1-x0);
-          mesh->verts[n][1]=y0*delta + y1*(1-delta);
         }
     }
+  else /* 3D */
+    {
+      for(ijk[2]=0; ijk[2]<grid->sizes[2]; ++ijk[2])
+        {
+          unsigned n_corner, n_interior, n, n_in;
+          ijk[0]=corner;
+          n_corner=RegularMeshUtils_Node_3DTo1D(mesh,ijk);
+          ijk[0]=inside;
+          n_interior=RegularMeshUtils_Node_3DTo1D(mesh,ijk);
+          if(Mesh_GlobalToDomain(mesh,MT_VERTEX,n_corner,&n)
+             && Mesh_GlobalToDomain(mesh,MT_VERTEX,n_interior,&n_in))
+            {
+              IJK ijk_off;
+              unsigned n_temp_off, n_temp_off_in, n_off, n_off_in;
+              double *crds[4];
+              double temp, pnt[2];
+              int offset=1;
+
+              crds[0]=mesh->verts[n];
+              crds[1]=mesh->verts[n_in];
+
+              ijk_off[0]=ijk[0];
+              ijk_off[1]=ijk[1];
+
+              if(crds[0][2]<crds[1][2])
+                offset=-1;
+
+              ijk_off[2]=ijk[2]+offset;
+              n_temp_off=RegularMeshUtils_Node_3DTo1D(mesh,ijk_off);
+              n_temp_off_in=RegularMeshUtils_Node_3DTo1D(mesh,ijk_off);
+              if(!Mesh_GlobalToDomain(mesh,MT_VERTEX,n_temp_off,&n_off_in)
+                 || !Mesh_GlobalToDomain(mesh,MT_VERTEX,n_temp_off_in,&n_off_in))
+                {
+                  printf("Can not map the neighbor of this coordinate to the local grid.\n %d %d %d %d",
+                         ijk[0],ijk[1],ijk[2], ijk_off[2]);
+                  abort();
+                }
+              crds[2]=mesh->verts[n_off_in];
+              crds[3]=mesh->verts[n_off];
+
+              pnt[0]=side_coord;
+              pnt[2]=crds[0][2];
+              if(!_EulerDeform_QuadYInterp(crds,pnt,&temp))
+                {
+                  printf("Unable to interpolate to this point\n %g %g %d %d %d\n",
+                         pnt[0],pnt[1],ijk[0],ijk[1],ijk[2]);
+                  abort();
+                }
+              mesh->verts[n][0]=side_coord;
+              mesh->verts[n][1]=temp;
+            }
+        }
+    }
 }
 
 void EulerDeform_Remesh( TimeIntegratee* crdAdvector, EulerDeform_Context* edCtx ) {