OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
e8ebf0840b
opm-core/opm/core/grid/cart_grid.c
Andreas Lauser 884c5ab027 make config.h the first header to be included in any compile unit 
this is required for consistency amongst the compile units which are
linked into the same library and seems to be forgotten quite
frequently.
2013-04-10 12:56:14 +02:00

709 lines
18 KiB
C
Raw Blame History

/*===========================================================================
//
// File: cart_grid.c
//
// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
//
//==========================================================================*/
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
Copyright 2011 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <stdlib.h>
#include <stdio.h>
#include <opm/core/grid.h>
#include <opm/core/grid/cornerpoint_grid.h>
#include <opm/core/grid/cart_grid.h>
static struct UnstructuredGrid *allocate_cart_grid_3d(int nx, int ny, int nz);
static void fill_cart_topology_3d(struct UnstructuredGrid *G);
static void fill_cart_geometry_3d(struct UnstructuredGrid *G,
const double *x,
const double *y,
const double *z);
static void
fill_layered_geometry_3d(struct UnstructuredGrid *G,
const double *x,
const double *y,
const double *z,
const double *depthz);
struct UnstructuredGrid *
create_grid_cart3d(int nx, int ny, int nz)
{
return create_grid_hexa3d(nx, ny, nz, 1.0, 1.0, 1.0);
}
struct UnstructuredGrid *
create_grid_hexa3d(int nx, int ny, int nz,
double dx, double dy, double dz)
{
int i;
double *x, *y, *z;
struct UnstructuredGrid *G;
x = malloc((nx + 1) * sizeof *x);
y = malloc((ny + 1) * sizeof *y);
z = malloc((nz + 1) * sizeof *z);
if ((x == NULL) || (y == NULL) || (z == NULL)) {
G = NULL;
} else {
for (i = 0; i < nx + 1; i++) { x[i] = i * dx; }
for (i = 0; i < ny + 1; i++) { y[i] = i * dy; }
for (i = 0; i < nz + 1; i++) { z[i] = i * dz; }
G = create_grid_tensor3d(nx, ny, nz, x, y, z,
(const double *) NULL);
}
free(z); free(y); free(x);
return G;
}
/* --------------------------------------------------------------------- */
static struct UnstructuredGrid *allocate_cart_grid_2d(int nx, int ny);
static void fill_cart_topology_2d(struct UnstructuredGrid *G);
static void fill_cart_geometry_2d(struct UnstructuredGrid *G,
const double *x,
const double *y);
struct UnstructuredGrid*
create_grid_cart2d(int nx, int ny, double dx, double dy)
{
int i;
double *x, *y;
struct UnstructuredGrid *G;
x = malloc((nx + 1) * sizeof *x);
y = malloc((ny + 1) * sizeof *y);
if ((x == NULL) || (y == NULL)) {
G = NULL;
} else {
for (i = 0; i < nx + 1; i++) { x[i] = i*dx; }
for (i = 0; i < ny + 1; i++) { y[i] = i*dy; }
G = create_grid_tensor2d(nx, ny, x, y);
}
free(y); free(x);
return G;
}
/* --------------------------------------------------------------------- */
struct UnstructuredGrid *
create_grid_tensor2d(int nx, int ny, const double *x, const double *y)
{
struct UnstructuredGrid *G;
G = allocate_cart_grid_2d(nx, ny);
if (G != NULL)
{
fill_cart_topology_2d(G);
fill_cart_geometry_2d(G, x, y);
}
return G;
}
/* --------------------------------------------------------------------- */
struct UnstructuredGrid *
create_grid_tensor3d(int nx ,
int ny ,
int nz ,
const double *x ,
const double *y ,
const double *z ,
const double *depthz)
{
struct UnstructuredGrid *G;
G = allocate_cart_grid_3d(nx, ny, nz);
if (G != NULL)
{
fill_cart_topology_3d(G);
if (depthz == NULL) {
fill_cart_geometry_3d(G, x, y, z);
}
else {
fill_layered_geometry_3d(G, x, y, z, depthz);
}
}
return G;
}
/* --------------------------------------------------------------------- */
/* Static functions follow: */
/* --------------------------------------------------------------------- */
static struct UnstructuredGrid *
allocate_cart_grid(size_t ndims ,
size_t ncells,
size_t nfaces,
size_t nnodes)
{
size_t nfacenodes, ncellfaces;
nfacenodes = nfaces * (2 * (ndims - 1));
ncellfaces = ncells * (2 * ndims);
return allocate_grid(ndims, ncells, nfaces,
nfacenodes, ncellfaces, nnodes);
}
static struct UnstructuredGrid*
allocate_cart_grid_3d(int nx, int ny, int nz)
{
struct UnstructuredGrid *G;
int Nx, Ny, Nz;
int nxf, nyf, nzf;
int ncells, nfaces, nnodes;
Nx = nx + 1;
Ny = ny + 1;
Nz = nz +1;
nxf = Nx * ny * nz;
nyf = nx * Ny * nz;
nzf = nx * ny * Nz;
ncells = nx * ny * nz ;
nfaces = nxf + nyf + nzf;
nnodes = Nx * Ny * Nz ;
G = allocate_cart_grid(3, ncells, nfaces, nnodes);
if (G != NULL)
{
G->dimensions = 3 ;
G->cartdims[0] = nx;
G->cartdims[1] = ny;
G->cartdims[2] = nz;
G->number_of_cells = ncells;
G->number_of_faces = nfaces;
G->number_of_nodes = nnodes;
}
return G;
}
/* --------------------------------------------------------------------- */
static void
fill_cart_topology_3d(struct UnstructuredGrid *G)
{
int nx, ny, nz;
int Nx, Ny;
int nxf, nyf;
int i,j,k;
int *cfaces, *cfacepos, *fnodes, *fnodepos, *fcells;
nx = G->cartdims[0];
ny = G->cartdims[1];
nz = G->cartdims[2];
Nx = nx+1;
Ny = ny+1;
nxf = Nx*ny*nz;
nyf = nx*Ny*nz;
cfaces = G->cell_faces;
cfacepos = G->cell_facepos;
cfacepos[0] = 0;
for (k=0; k<nz; ++k) {
for (j=0; j<ny; ++j) {
for (i=0; i<nx; ++i) {
*cfaces++ = i+ Nx*(j+ ny* k );
*cfaces++ = i+1+Nx*(j+ ny* k );
*cfaces++ = i+ nx*(j+ Ny* k ) +nxf;
*cfaces++ = i+ nx*(j+1+Ny* k ) +nxf;
*cfaces++ = i+ nx*(j+ ny* k ) +nxf+nyf;
*cfaces++ = i+ nx*(j+ ny*(k+1)) +nxf+nyf;
cfacepos[1] = cfacepos[0]+6;
++cfacepos;
}
}
}
for (k = 0; k < nx * ny * nz; ++k) {
for (i = 0; i < 6; ++i) {
G->cell_facetag[k*6 + i] = i;
}
}
fnodes = G->face_nodes;
fnodepos = G->face_nodepos;
fcells = G->face_cells;
fnodepos[0] = 0;
/* Faces with x-normal */
for (k=0; k<nz; ++k) {
for (j=0; j<ny; ++j) {
for (i=0; i<nx+1; ++i) {
*fnodes++ = i+Nx*(j + Ny * k );
*fnodes++ = i+Nx*(j+1 + Ny * k );
*fnodes++ = i+Nx*(j+1 + Ny *(k+1));
*fnodes++ = i+Nx*(j + Ny *(k+1));
fnodepos[1] = fnodepos[0] + 4;
++fnodepos;
if (i==0) {
*fcells++ = -1;
*fcells++ = i+nx*(j+ny*k);
}
else if (i == nx) {
*fcells++ = i-1+nx*(j+ny*k);
*fcells++ = -1;
}
else {
*fcells++ = i-1 + nx*(j+ny*k);
*fcells++ = i + nx*(j+ny*k);
}
}
}
}
/* Faces with y-normal */
for (k=0; k<nz; ++k) {
for (j=0; j<ny+1; ++j) {
for (i=0; i<nx; ++i) {
*fnodes++ = i+ Nx*(j + Ny * k );
*fnodes++ = i + Nx*(j + Ny *(k+1));
*fnodes++ = i+1 + Nx*(j + Ny *(k+1));
*fnodes++ = i+1 + Nx*(j + Ny * k );
fnodepos[1] = fnodepos[0] + 4;
++fnodepos;
if (j==0) {
*fcells++ = -1;
*fcells++ = i+nx*(j+ny*k);
}
else if (j == ny) {
*fcells++ = i+nx*(j-1+ny*k);
*fcells++ = -1;
}
else {
*fcells++ = i+nx*(j-1+ny*k);
*fcells++ = i+nx*(j+ny*k);
}
}
}
}
/* Faces with z-normal */
for (k=0; k<nz+1; ++k) {
for (j=0; j<ny; ++j) {
for (i=0; i<nx; ++i) {
*fnodes++ = i+ Nx*(j + Ny * k);
*fnodes++ = i+1 + Nx*(j + Ny * k);
*fnodes++ = i+1 + Nx*(j+1 + Ny * k);
*fnodes++ = i+ Nx*(j+1 + Ny * k);
fnodepos[1] = fnodepos[0] + 4;
++fnodepos;
if (k==0) {
*fcells++ = -1;
*fcells++ = i+nx*(j+ny*k);
}
else if (k == nz) {
*fcells++ = i+nx*(j+ny*(k-1));
*fcells++ = -1;
}
else {
*fcells++ = i+nx*(j+ny*(k-1));
*fcells++ = i+nx*(j+ny*k);
}
}
}
}
}
/* --------------------------------------------------------------------- */
static void
fill_cart_geometry_3d(struct UnstructuredGrid *G,
const double *x,
const double *y,
const double *z)
{
int nx, ny, nz;
int i,j,k;
double dx, dy, dz;
double *coord, *ccentroids, *cvolumes;
double *fnormals, *fcentroids, *fareas;
nx = G->cartdims[0];
ny = G->cartdims[1];
nz = G->cartdims[2];
ccentroids = G->cell_centroids;
cvolumes = G->cell_volumes;
for (k=0; k<nz; ++k) {
for (j=0; j<ny; ++j) {
for (i=0; i<nx; ++i) {
*ccentroids++ = (x[i] + x[i + 1]) / 2.0;
*ccentroids++ = (y[j] + y[j + 1]) / 2.0;
*ccentroids++ = (z[k] + z[k + 1]) / 2.0;
dx = x[i + 1] - x[i];
dy = y[j + 1] - y[j];
dz = z[k + 1] - z[k];
*cvolumes++ = dx * dy * dz;
}
}
}
fnormals = G->face_normals;
fcentroids = G->face_centroids;
fareas = G->face_areas;
/* Faces with x-normal */
for (k=0; k<nz; ++k) {
for (j=0; j<ny; ++j) {
for (i=0; i<nx+1; ++i) {
dy = y[j + 1] - y[j];
dz = z[k + 1] - z[k];
*fnormals++ = dy * dz;
*fnormals++ = 0;
*fnormals++ = 0;
*fcentroids++ = x[i];
*fcentroids++ = (y[j] + y[j + 1]) / 2.0;
*fcentroids++ = (z[k] + z[k + 1]) / 2.0;
*fareas++ = dy * dz;
}
}
}
/* Faces with y-normal */
for (k=0; k<nz; ++k) {
for (j=0; j<ny+1; ++j) {
for (i=0; i<nx; ++i) {
dx = x[i + 1] - x[i];
dz = z[k + 1] - z[k];
*fnormals++ = 0;
*fnormals++ = dx * dz;
*fnormals++ = 0;
*fcentroids++ = (x[i] + x[i + 1]) / 2.0;
*fcentroids++ = y[j];
*fcentroids++ = (z[k] + z[k + 1]) / 2.0;
*fareas++ = dx * dz;
}
}
}
/* Faces with z-normal */
for (k=0; k<nz+1; ++k) {
for (j=0; j<ny; ++j) {
for (i=0; i<nx; ++i) {
dx = x[i + 1] - x[i];
dy = y[j + 1] - y[j];
*fnormals++ = 0;
*fnormals++ = 0;
*fnormals++ = dx * dy;
*fcentroids++ = (x[i] + x[i + 1]) / 2.0;
*fcentroids++ = (y[j] + y[j + 1]) / 2.0;
*fcentroids++ = z[k];
*fareas++ = dx * dy;
}
}
}
coord = G->node_coordinates;
for (k=0; k<nz+1; ++k) {
for (j=0; j<ny+1; ++j) {
for (i=0; i<nx+1; ++i) {
*coord++ = x[i];
*coord++ = y[j];
*coord++ = z[k];
}
}
}
}
/* ---------------------------------------------------------------------- */
static void
fill_layered_geometry_3d(struct UnstructuredGrid *G,
const double *x,
const double *y,
const double *z,
const double *depthz)
{
int i , j , k ;
int nx, ny, nz;
const double *depth;
double *coord;
nx = G->cartdims[0]; ny = G->cartdims[1]; nz = G->cartdims[2];
coord = G->node_coordinates;
for (k = 0; k < nz + 1; k++) {
depth = depthz;
for (j = 0; j < ny + 1; j++) {
for (i = 0; i < nx + 1; i++) {
*coord++ = x[i];
*coord++ = y[j];
*coord++ = z[k] + *depth++;
}
}
}
compute_geometry(G);
}
/* --------------------------------------------------------------------- */
static struct UnstructuredGrid*
allocate_cart_grid_2d(int nx, int ny)
{
int nxf, nyf;
int Nx , Ny ;
int ncells, nfaces, nnodes;
struct UnstructuredGrid *G;
Nx = nx + 1;
Ny = ny + 1;
nxf = Nx * ny;
nyf = nx * Ny;
ncells = nx * ny ;
nfaces = nxf + nyf;
nnodes = Nx * Ny ;
G = allocate_cart_grid(2, ncells, nfaces, nnodes);
if (G != NULL)
{
G->dimensions = 2 ;
G->cartdims[0] = nx;
G->cartdims[1] = ny;
G->cartdims[2] = 1 ;
G->number_of_cells = ncells;
G->number_of_faces = nfaces;
G->number_of_nodes = nnodes;
}
return G;
}
/* --------------------------------------------------------------------- */
static void
fill_cart_topology_2d(struct UnstructuredGrid *G)
{
int i,j;
int nx, ny;
int nxf;
int Nx;
int *fnodes, *fnodepos, *fcells, *cfaces, *cfacepos;
cfaces = G->cell_faces;
cfacepos = G->cell_facepos;
nx = G->cartdims[0];
ny = G->cartdims[1];
Nx = nx + 1;
nxf = Nx * ny;
cfacepos[0] = 0;
for (j=0; j<ny; ++j) {
for (i=0; i<nx; ++i) {
*cfaces++ = i+ Nx*j;
*cfaces++ = i+ nx*j +nxf;
*cfaces++ = i+1+Nx*j;
*cfaces++ = i+ nx*(j+1)+nxf;
cfacepos[1] = cfacepos[0]+4;
++cfacepos;
}
}
for (j = 0; j < nx * ny; ++j) {
G->cell_facetag[j*4 + 0] = 0;
G->cell_facetag[j*4 + 1] = 2;
G->cell_facetag[j*4 + 2] = 1;
G->cell_facetag[j*4 + 3] = 3;
}
fnodes = G->face_nodes;
fnodepos = G->face_nodepos;
fcells = G->face_cells;
fnodepos[0] = 0;
/* Faces with x-normal */
for (j=0; j<ny; ++j) {
for (i=0; i<nx+1; ++i) {
*fnodes++ = i+Nx*j;
*fnodes++ = i+Nx*(j+1);
fnodepos[1] = fnodepos[0] + 2;
++fnodepos;
if (i==0) {
*fcells++ = -1;
*fcells++ = i+nx*j;
}
else if (i == nx) {
*fcells++ = i-1+nx*j;
*fcells++ = -1;
}
else {
*fcells++ = i-1 + nx*j;
*fcells++ = i + nx*j;
}
}
}
/* Faces with y-normal */
for (j=0; j<ny+1; ++j) {
for (i=0; i<nx; ++i) {
*fnodes++ = i+1 + Nx*j;
*fnodes++ = i+ Nx*j;
fnodepos[1] = fnodepos[0] + 2;
++fnodepos;
if (j==0) {
*fcells++ = -1;
*fcells++ = i+nx*j;
}
else if (j == ny) {
*fcells++ = i+nx*(j-1);
*fcells++ = -1;
}
else {
*fcells++ = i+nx*(j-1);
*fcells++ = i+nx*j;
}
}
}
}
/* --------------------------------------------------------------------- */
static void
fill_cart_geometry_2d(struct UnstructuredGrid *G,
const double *x,
const double *y)
{
int i,j;
int nx, ny;
double dx, dy;
double *coord, *ccentroids, *cvolumes;
double *fnormals, *fcentroids, *fareas;
nx = G->cartdims[0];
ny = G->cartdims[1];
ccentroids = G->cell_centroids;
cvolumes = G->cell_volumes;
for (j=0; j<ny; ++j) {
for (i=0; i<nx; ++i) {
*ccentroids++ = (x[i] + x[i + 1]) / 2.0;
*ccentroids++ = (y[j] + y[j + 1]) / 2.0;
dx = x[i + 1] - x[i];
dy = y[j + 1] - y[j];
*cvolumes++ = dx * dy;
}
}
fnormals = G->face_normals;
fcentroids = G->face_centroids;
fareas = G->face_areas;
/* Faces with x-normal */
for (j=0; j<ny; ++j) {
for (i=0; i<nx+1; ++i) {
dy = y[j + 1] - y[j];
*fnormals++ = dy;
*fnormals++ = 0;
*fcentroids++ = x[i];
*fcentroids++ = (y[j] + y[j + 1]) / 2.0;
*fareas++ = dy;
}
}
/* Faces with y-normal */
for (j=0; j<ny+1; ++j) {
for (i=0; i<nx; ++i) {
dx = x[i + 1] - x[i];
*fnormals++ = 0;
*fnormals++ = dx;
*fcentroids++ = (x[i] + x[i + 1]) / 2.0;
*fcentroids++ = y[j];
*fareas++ = dx;
}
}
coord = G->node_coordinates;
for (j=0; j<ny+1; ++j) {
for (i=0; i<nx+1; ++i) {
*coord++ = x[i];
*coord++ = y[j];
}
}
}
Reference in New Issue View Git Blame Copy Permalink