OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Compute properties for cells in 2D

Browse Source
This commit is contained in:
Roland Kaufmann 2013-01-21 10:17:00 +01:00
parent 701526596d
commit 4fb25b91ef
1 changed files with 100 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

100
opm/core/grid/cpgpreprocess/geometry.c
View File

@ -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);