From 7d8511fdaf2c2618e4fe962a53c6253b7e1d7e36 Mon Sep 17 00:00:00 2001 From: Roland Kaufmann Date: Thu, 10 Jan 2013 13:28:22 +0100 Subject: [PATCH 1/5] Prepare for processing two-dim. grid Allow a two-dimensional grid to be sent to compute_geometry() without the assertion triggering. Route three-dimensional grids to the existing implementation. --- opm/core/grid/cornerpoint_grid.c | 2 - opm/core/grid/cpgpreprocess/geometry.c | 64 +++++++++++++++++++++----- 2 files changed, 53 insertions(+), 13 deletions(-) diff --git a/opm/core/grid/cornerpoint_grid.c b/opm/core/grid/cornerpoint_grid.c index 2ed774e0..070379d7 100644 --- a/opm/core/grid/cornerpoint_grid.c +++ b/opm/core/grid/cornerpoint_grid.c @@ -138,8 +138,6 @@ void compute_geometry(struct UnstructuredGrid *g) assert (g != NULL); if (g!=NULL) { - assert (g->dimensions == 3); - assert (g->face_centroids != NULL); assert (g->face_normals != NULL); assert (g->face_areas != NULL); diff --git a/opm/core/grid/cpgpreprocess/geometry.c b/opm/core/grid/cpgpreprocess/geometry.c index a49fc874..3cbd626f 100644 --- a/opm/core/grid/cpgpreprocess/geometry.c +++ b/opm/core/grid/cpgpreprocess/geometry.c @@ -5,6 +5,7 @@ #include #include #include "geometry.h" +#include /* ------------------------------------------------------------------ */ static void @@ -25,16 +26,16 @@ norm(const double w[3]) } - /* ------------------------------------------------------------------ */ -void -compute_face_geometry(int ndims, double *coords, int nfaces, - int *nodepos, int *facenodes, double *fnormals, - double *fcentroids, double *fareas) +static void +compute_face_geometry_3d(double *coords, int nfaces, + int *nodepos, int *facenodes, double *fnormals, + double *fcentroids, double *fareas) /* ------------------------------------------------------------------ */ { /* Assume 3D for now */ + const int ndims = 3; int f; double x[3]; double u[3]; @@ -112,15 +113,34 @@ compute_face_geometry(int ndims, double *coords, int nfaces, /* ------------------------------------------------------------------ */ void -compute_cell_geometry(int ndims, double *coords, - int *nodepos, int *facenodes, int *neighbors, - double *fnormals, - double *fcentroids, - int ncells, int *facepos, int *cellfaces, - double *ccentroids, double *cvolumes) +compute_face_geometry(int ndims, double *coords, int nfaces, + int *nodepos, int *facenodes, double *fnormals, + double *fcentroids, double *fareas) /* ------------------------------------------------------------------ */ { + if (ndims == 3) + { + compute_face_geometry_3d(coords, nfaces, nodepos, facenodes, + fnormals, fcentroids, fareas); + } + else + { + assert(0); + } +} + +/* ------------------------------------------------------------------ */ +static void +compute_cell_geometry_3d(double *coords, + int *nodepos, int *facenodes, int *neighbors, + double *fnormals, + double *fcentroids, + int ncells, int *facepos, int *cellfaces, + double *ccentroids, double *cvolumes) +/* ------------------------------------------------------------------ */ +{ + const int ndims = 3; int i,k, f,c; int face,node; double x[3]; @@ -240,3 +260,25 @@ compute_cell_geometry(int ndims, double *coords, cvolumes[c] = volume; } } + +/* ------------------------------------------------------------------ */ +void +compute_cell_geometry(int ndims, double *coords, + int *nodepos, int *facenodes, int *neighbors, + double *fnormals, + double *fcentroids, + int ncells, int *facepos, int *cellfaces, + double *ccentroids, double *cvolumes) +/* ------------------------------------------------------------------ */ +{ + if (ndims == 3) + { + compute_cell_geometry_3d(coords, nodepos, facenodes, + neighbors, fnormals, fcentroids, ncells, + facepos, cellfaces, ccentroids, cvolumes); + } + else + { + assert(0); + } +} From 701526596d76080797494937e06c708ec5e3d319 Mon Sep 17 00:00:00 2001 From: Roland Kaufmann Date: Fri, 18 Jan 2013 14:45:52 +0100 Subject: [PATCH 2/5] Compute properties for edges in 2D --- opm/core/grid/cpgpreprocess/geometry.c | 71 ++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/opm/core/grid/cpgpreprocess/geometry.c b/opm/core/grid/cpgpreprocess/geometry.c index 3cbd626f..1e3eeb92 100644 --- a/opm/core/grid/cpgpreprocess/geometry.c +++ b/opm/core/grid/cpgpreprocess/geometry.c @@ -110,6 +110,71 @@ compute_face_geometry_3d(double *coords, int nfaces, } } +/* ------------------------------------------------------------------ */ +static void +compute_edge_geometry_2d( + /* in */ double *node_coords, + /* in */ int num_edges, + /* in */ int *edge_node_pos, + /* in */ int *edge_nodes, + /* out */ double *edge_normals, + /* out */ double *edge_midpoints, + /* out */ double *edge_lengths) +{ + 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 edge; /* edge index */ + int a_nod, b_nod; /* node indices */ + double a_x, a_y, b_x, b_y; /* node coordinates */ + double v_x, v_y; /* vector elements */ + + /* decompose each edge into a tuple (a,b) between two points and + * compute properties for that face. hopefully the host has enough + * cache pages to keep both input and output at the same time, and + * registers for all the local variables */ + for (edge = 0; edge < num_edges; ++edge) + { + /* an edge in 2D can only have starting and ending point + * check that there are exactly two nodes till the next edge */ + assert (edge_node_pos[edge + 1] - edge_node_pos[edge] == num_dims); + + /* 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]; + + /* extract individual coordinates for the points */ + a_x = node_coords[a_nod * num_dims + x_ofs]; + a_y = node_coords[a_nod * num_dims + y_ofs]; + b_x = node_coords[b_nod * num_dims + x_ofs]; + b_y = node_coords[b_nod * num_dims + y_ofs]; + + /* compute edge center -- average of node coordinates */ + edge_midpoints[edge * num_dims + x_ofs] = (a_x + b_x) * 0.5; + edge_midpoints[edge * num_dims + y_ofs] = (a_y + b_y) * 0.5; + + /* vector from first to last point */ + v_x = b_x - a_x; + v_y = b_y - a_y; + + /* two-dimensional (unary) cross product analog that makes the + * "triple" (dot-cross) product zero, i.e. it's a normal; the + * direction of this vector is such that it will be pointing + * inwards when enumerating nodes clock-wise */ + edge_normals[edge * num_dims + x_ofs] = +v_y; + edge_normals[edge * num_dims + y_ofs] = -v_x; + + /* Euclidian norm in two dimensions is magnitude of edge */ + edge_lengths[edge] = sqrt(v_x*v_x + v_y*v_y); + } +} /* ------------------------------------------------------------------ */ void @@ -123,6 +188,12 @@ compute_face_geometry(int ndims, double *coords, int nfaces, compute_face_geometry_3d(coords, nfaces, nodepos, facenodes, fnormals, fcentroids, fareas); } + else if (ndims == 2) + { + /* two-dimensional interfaces are called 'edges' */ + compute_edge_geometry_2d(coords, nfaces, nodepos, facenodes, + fnormals, fcentroids, fareas); + } else { assert(0); From 4fb25b91efe86b05331511087c60f2168c8e6a8a Mon Sep 17 00:00:00 2001 From: Roland Kaufmann Date: Mon, 21 Jan 2013 10:17:00 +0100 Subject: [PATCH 3/5] 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); From ccf3edb00f2bbfaaafd19d235821118bb876b6fb Mon Sep 17 00:00:00 2001 From: Roland Kaufmann Date: Mon, 21 Jan 2013 22:40:41 +0100 Subject: [PATCH 4/5] Set size properties in grid when allocating structure There is little chance that the client want to set them otherwise since that would make the structure inconsistent, and having all the clients do the work is just needless duplicate labour. --- opm/core/grid.c | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/opm/core/grid.c b/opm/core/grid.c index 12a851b1..6271bc85 100644 --- a/opm/core/grid.c +++ b/opm/core/grid.c @@ -81,6 +81,12 @@ allocate_grid(size_t ndims , G = create_grid_empty(); if (G != NULL) { + /* Grid fields ---------------------------------------- */ + G->dimensions = ndims; + G->number_of_cells = ncells; + G->number_of_faces = nfaces; + G->number_of_nodes = nnodes; + /* Node fields ---------------------------------------- */ nel = nnodes * ndims; G->node_coordinates = malloc(nel * sizeof *G->node_coordinates); From 591cf88ccd0669d98c9e288fe32b409eec94c5e8 Mon Sep 17 00:00:00 2001 From: Roland Kaufmann Date: Mon, 21 Jan 2013 23:06:08 +0100 Subject: [PATCH 5/5] Unit test the geometry computations for 2D grids --- tests/test_geom2d.cpp | 195 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 195 insertions(+) create mode 100644 tests/test_geom2d.cpp diff --git a/tests/test_geom2d.cpp b/tests/test_geom2d.cpp new file mode 100644 index 00000000..a8ffbe2d --- /dev/null +++ b/tests/test_geom2d.cpp @@ -0,0 +1,195 @@ +/* Copyright 2013 Uni Research AS + * This file is licensed under GPL3, see http://www.opm-project.org/ +*/ +#include + +/* --- Boost.Test boilerplate --- */ +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE // Suppress own messages when throw()ing + +#define BOOST_TEST_MODULE CompGeo2DTest +#include +#include + +/* --- our own headers --- */ +#include +#include +#include +#include /* compute_geometry */ + +using namespace std; + +/* Test properties on a grid that looks like this: + * + * /------ + * / | | + * < | | + * \ | | + * \------ + * + */ +struct BackspaceGrid { + UnstructuredGrid* g; + + BackspaceGrid() + { + g = allocate_grid( + 2, /* ndims */ + 2, /* number of cells */ + 6, /* number of edges */ + 6*2, /* six edges with two nodes each */ + 3+4, /* one triangle and one quadrilateral */ + 5); /* five nodes total */ + /* nodes */ + double nodes[] = { + -1., 0., /* node 0 */ + 0., 1., /* node 1 */ + 1., 1., /* node 2 */ + 1., -1., /* node 3 */ + 0., -1., /* node 4 */ + }; + copy(nodes, nodes+sizeof(nodes)/sizeof(nodes[0]), g->node_coordinates); + /* edges */ + int edges[] = { + 0, 1, /* edge 0 */ + 1, 2, /* edge 1 */ + 2, 3, /* edge 2 */ + 3, 4, /* edge 3 */ + 4, 0, /* edge 4 */ + 4, 1, /* edge 5 */ + }; + copy(edges, edges+sizeof(edges)/sizeof(edges[0]), g->face_nodes); + /* starting index in map for each edge */ + int edge_pos[] = { + 0, 2, 4, 6, 8, 10, 12, + }; + copy(edge_pos, edge_pos+sizeof(edge_pos)/sizeof(edge_pos[0]), g->face_nodepos); + /* topology, clock-wise ordering */ + int neighbours[] = { + -1, 0, /* edge 0, between boundary and cell 0 */ + -1, 1, /* edge 1, between boundary and cell 1 */ + -1, 1, /* edge 2, between boundary and cell 1 */ + -1, 1, /* edge 3, between boundary and cell 1 */ + -1, 0, /* edge 4, between boundary and cell 0 */ + 0, 1, /* edge 5, between cell 0 and cell 1 */ + }; + copy(neighbours, neighbours+sizeof(neighbours)/sizeof(neighbours[0]), g->face_cells); + /* cells */ + int cells[] = { + 0, 5, 4, /* cell 0, clockwise */ + 1, 2, 3, 5, /* cell 1, clockwise */ + }; + copy(cells, cells+sizeof(cells)/sizeof(cells[0]), g->cell_faces); + /* starting index in map for each cell */ + int cell_pos[] = { + 0, 3, 7, + }; + copy(cell_pos, cell_pos+sizeof(cell_pos)/sizeof(cell_pos[0]), g->cell_facepos); + /* everything interesting actually happens here, for all tests */ + compute_geometry(g); + } + + ~BackspaceGrid() + { + destroy_grid(g); + } +}; + +BOOST_FIXTURE_TEST_SUITE(CompGeo2D, BackspaceGrid) + +BOOST_AUTO_TEST_CASE(edgeMidpoints) +{ + double midpoints[] = { + -0.5, 0.5, /* edge 0 */ + 0.5, 1., /* edge 1 */ + 1., 0., /* edge 2 */ + 0.5, -1., /* edge 3 */ + -0.5, -0.5, /* edge 4 */ + 0., 0., /* edge 5 */ + }; + BOOST_REQUIRE (sizeof(midpoints)/sizeof(midpoints[0]) == + g->number_of_faces * g->dimensions); + for (int edge = 0; edge < g->number_of_faces; ++edge) + { + for (int dim = 0; dim < g->dimensions; ++dim) + { + BOOST_REQUIRE_CLOSE (g->face_centroids[edge*g->dimensions+dim], + midpoints[edge*g->dimensions+dim], 0.001); + } + } +} + +BOOST_AUTO_TEST_CASE(edgeNormals) +{ + /* inward normal when doing clockwise enumeration */ + double normals[] = { + 1., -1., /* edge 0 */ + 0., -1., /* edge 1 */ + -2., 0., /* edge 2 */ + 0., 1., /* edge 3 */ + 1., 1., /* edge 4 */ + 2., 0., /* edge 5 */ + }; + BOOST_REQUIRE (sizeof(normals)/sizeof(normals[0]) == + g->number_of_faces * g->dimensions); + for (int edge = 0; edge < g->number_of_faces; ++edge) + { + for (int dim = 0; dim < g->dimensions; ++dim) + { + BOOST_REQUIRE_CLOSE (g->face_normals[edge*g->dimensions+dim], + normals[edge*g->dimensions+dim], 0.001); + } + } +} + +BOOST_AUTO_TEST_CASE(edgeLengths) +{ + double lengths[] = { + 1.4142, /* edge 0 */ + 1., /* edge 1 */ + 2., /* edge 2 */ + 1., /* edge 3 */ + 1.4142, /* edge 4 */ + 2., /* edge 5 */ + }; + BOOST_REQUIRE (sizeof(lengths)/sizeof(lengths[0]) == g->number_of_faces); + for (int edge = 0; edge < g->number_of_faces; ++edge) + { + BOOST_REQUIRE_CLOSE (g->face_areas[edge], lengths[edge], 0.001); + } +} + +BOOST_AUTO_TEST_CASE(cellAreas) +{ + double areas[] = { + 1., 2., + }; + BOOST_REQUIRE (sizeof(areas)/sizeof(areas[0]) == g->number_of_cells); + for (int cell = 0; cell < g->number_of_cells; ++cell) + { + BOOST_REQUIRE_CLOSE (g->cell_volumes[cell], areas[cell], 0.001); + } +} + +BOOST_AUTO_TEST_CASE(cellCenters) +{ + double cellCenters[] = { + -1./3., 0.0, /* cell 0 */ + 1./2., 0.0, /* cell 1 */ + }; + BOOST_REQUIRE (sizeof(cellCenters)/sizeof(cellCenters[0]) == + g->number_of_cells * g->dimensions); + for (int cell = 0; cell < g->number_of_cells; ++cell) + { + for (int dim = 0; dim < g->dimensions; ++dim) + { + BOOST_REQUIRE_CLOSE (g->cell_centroids[cell*g->dimensions+dim], + cellCenters[cell*g->dimensions+dim], 0.001); + } + } +} + +BOOST_AUTO_TEST_SUITE_END()