[cig-commits] r14503 - seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust

[tan2 at geodynamics.org]

Author: tan2
Date: 2009-03-27 17:19:00 -0700 (Fri, 27 Mar 2009)
New Revision: 14503

Modified:
   seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/read_citcoms_data.c
Log:
Fixed a bug in locating elements

Modified: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/read_citcoms_data.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/read_citcoms_data.c	2009-03-28 00:09:46 UTC (rev 14502)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/read_citcoms_data.c	2009-03-28 00:19:00 UTC (rev 14503)
@@ -723,10 +723,11 @@
 
 
 
-static int is_inside_element(const model_t *model,
-                             double x, double y, double z, int elem)
+static int is_inside_column(const model_t *model,
+                            double x, double y, double z, int elem)
 {
-    /* whether point (x, y, z) is inside element elem */
+    /* Return whether point (x, y, z) is inside the column
+       defined by element elem. The radial range is not checked! */
 
     double corner_x[4], corner_y[4], corner_z[4];
     int i, corners[4];
@@ -816,7 +817,7 @@
     for(i=0; i<4; i++) {
         if(neighbor_elements[i] >= 0 &&
            neighbor_elements[i] < nel &&
-           is_inside_element(model, x, y, z, neighbor_elements[i]))
+           is_inside_column(model, x, y, z, neighbor_elements[i]))
             /* found */
             return neighbor_elements[i];
     }
@@ -837,7 +838,7 @@
     for(nx=0; nx<nox-1; nx++)
         for(ny=0; ny<noy-1; ny++) {
             elem = ijk2elem(nx, ny, nz);
-            if(is_inside_element(model, x, y, z, elem))
+            if(is_inside_column(model, x, y, z, elem))
                 /* found */
                 return elem;
         }
@@ -1234,7 +1235,7 @@
     double x, y, z;
     double xx, yy, zz;
     domain_t *dom;
-    int i;
+    int i, nz;
 
     npoints ++;
 
@@ -1248,7 +1249,9 @@
     if(*hint_domain >= 0) {
         dom = &domains[*hint_domain];
 
-        if(is_inside_element(dom->model, x, y, z, *hint_elem)) {
+        nz = elem % (noz - 1);
+        if(nz == find_layer(dom->model, x, y, z) &&
+           is_inside_column(dom->model, x, y, z, *hint_elem)) {
             /* search is successful */
             hint_works ++;
             get_velocities(dom->model, *hint_elem, x, y, z, *radius,