From 4fb25b91efe86b05331511087c60f2168c8e6a8a Mon Sep 17 00:00:00 2001
From: Roland Kaufmann <roland.kaufmann@uni.no>
Date: Mon, 21 Jan 2013 10:17:00 +0100
Subject: [PATCH] Compute properties for cells in 2D

---
 opm/core/grid/cpgpreprocess/geometry.c | 100 +++++++++++++++++++++++++
 1 file changed, 100 insertions(+)

diff --git a/opm/core/grid/cpgpreprocess/geometry.c b/opm/core/grid/cpgpreprocess/geometry.c
index 1e3eeb92..769d5028 100644
--- a/opm/core/grid/cpgpreprocess/geometry.c
+++ b/opm/core/grid/cpgpreprocess/geometry.c
@@ -332,6 +332,100 @@ compute_cell_geometry_3d(double *coords,
    }
 }
 
+/* ------------------------------------------------------------------ */
+static void
+compute_cell_geometry_2d(
+      /* in  */ double *node_coords,
+      /* in  */ int     *edge_node_pos,
+      /* in  */ int    *edge_nodes,
+      /* in  */ double *edge_midpoints,
+      /* in  */ int     num_cells,
+      /* in  */ int    *cell_edge_pos,
+      /* in  */ int    *cell_edges,
+      /* out */ double *cell_centers,
+      /* out */ double *cell_areas)
+{
+   const int num_dims = 2;
+
+   /* offsets to each of the nodes in a compacted edge */
+   const int a_ofs = 0;
+   const int b_ofs = 1;
+
+   /* offsets to each dimension is a compacted point */
+   const int x_ofs = 0;
+   const int y_ofs = 1;
+
+   int cell;            /* cell index */
+   int num_nodes;       /* number of vertices in current cell */
+   int edge_ndx;        /* relative edge index within cell */
+   int edge;            /* absolute cell index */
+   double center_x;     /* x-coordinate for cell barycenter */
+   double center_y;     /* y-coordinate for cell barycenter */
+   double area;         /* (accumulated) cell area */
+   int a_nod, b_nod;    /* node indices for edge start and end points */
+   double a_x, a_y,
+          b_x, b_y;     /* vectors from center to edge points */
+
+   for (cell = 0; cell < num_cells; ++cell)
+   {
+      /* since the cell is a closed polygon, each point serves as the starting
+       * point of one edge and the ending point of another; thus there is as
+       * many vertices as there are edges */
+      num_nodes = cell_edge_pos[cell + 1] - cell_edge_pos[cell];
+
+      /* to enumerate all vertices of a cell, we would have to expand the
+       * edges and then remove duplicates. however, the centroid of each
+       * edge contains half of the two vertices that are incident on it. if
+       * we instead sum all the face centroids, we get the sum of all the
+       * vertices */
+      center_x = 0.;
+      center_y = 0.;
+      for (edge_ndx = cell_edge_pos[cell];
+           edge_ndx < cell_edge_pos[cell + 1]; ++edge_ndx)
+      {
+         edge = cell_edges[edge_ndx];
+         center_x += edge_midpoints[edge * num_dims + x_ofs];
+         center_y += edge_midpoints[edge * num_dims + y_ofs];
+      }
+      center_x /= (double) num_nodes;
+      center_y /= (double) num_nodes;
+      cell_centers[cell * num_dims + x_ofs] = center_x;
+      cell_centers[cell * num_dims + y_ofs] = center_y;
+
+      /* triangulate the polygon by introducing the cell center and then new
+       * internal edges from this center to the vertices. the total area of
+       * the cell is the sum of area of these sub-triangles */
+      area = 0.;
+      for (edge_ndx = cell_edge_pos[cell];
+           edge_ndx < cell_edge_pos[cell + 1]; ++edge_ndx)
+      {
+         /* indirect lookup of edge index (from array that contains all the
+          * edge indices for a certain cell) */
+         edge = cell_edges[edge_ndx];
+
+         /* get the first and last point on the edge */
+         a_nod = edge_nodes[edge_node_pos[edge] + a_ofs];
+         b_nod = edge_nodes[edge_node_pos[edge] + b_ofs];
+
+         /* vector from center to each of the nodes */
+         a_x = node_coords[a_nod * num_dims + x_ofs] - center_x;
+         a_y = node_coords[a_nod * num_dims + y_ofs] - center_y;
+         b_x = node_coords[b_nod * num_dims + x_ofs] - center_x;
+         b_y = node_coords[b_nod * num_dims + y_ofs] - center_y;
+
+         /* two-dimensional (binary) cross product analog that has length
+          * equal to the parallelogram spanned by the two vectors (but which
+          * is a scalar). the sign tells us the orientation between the nodes
+          * a and b, but we are not interested in that, just the area */
+         area += fabs(a_x * b_y - a_y * b_x);
+      }
+      /* we summed parallelograms which are twice the size of the triangles
+       * that make up the cell; divide out the half for all terms here */
+      area *= 0.5;
+      cell_areas[cell] = area;
+   }
+}
+
 /* ------------------------------------------------------------------ */
 void
 compute_cell_geometry(int ndims, double *coords,
@@ -348,6 +442,12 @@ compute_cell_geometry(int ndims, double *coords,
                                neighbors, fnormals, fcentroids, ncells,
                                facepos, cellfaces, ccentroids, cvolumes);
    }
+   else if (ndims == 2)
+   {
+      compute_cell_geometry_2d(coords, nodepos, facenodes, fcentroids,
+                               ncells, facepos, cellfaces, ccentroids,
+                               cvolumes);
+   }
    else
    {
       assert(0);