diff --git a/facetopology.c b/facetopology.c index a819fa3a..e9e1ec5f 100644 --- a/facetopology.c +++ b/facetopology.c @@ -13,23 +13,23 @@ //===========================================================================*/ /* -Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. -Copyright 2009, 2010 Statoil ASA. + Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. + Copyright 2009, 2010 Statoil ASA. -This file is part of The Open Reservoir Simulator Project (OpenRS). + This file is part of The Open Reservoir Simulator Project (OpenRS). -OpenRS 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. + OpenRS 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. -OpenRS 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. + OpenRS 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 OpenRS. If not, see . + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ #include @@ -61,132 +61,132 @@ along with OpenRS. If not, see . /* Determine face topology first, then compute intersection. */ /* All intersections that occur are present in the final face geometry.*/ static int *computeFaceTopology(int *a1, - int *a2, - int *b1, - int *b2, - int intersect[4], - int *faces) + int *a2, + int *b1, + int *b2, + int intersect[4], + int *faces) { - int mask[8] = {-1}; - int k; - int *f; - /* Which pillar points should we use? */ - if (a1[1] > b1[1]){ mask[0] = b1[1]; } else { mask[0] = a1[1]; } - if (a2[1] > b2[1]){ mask[2] = b2[1]; } else { mask[2] = a2[1]; } - if (a2[0] > b2[0]){ mask[4] = a2[0]; } else { mask[4] = b2[0]; } - if (a1[0] > b1[0]){ mask[6] = a1[0]; } else { mask[6] = b1[0]; } + int mask[8] = {-1}; + int k; + int *f; + /* Which pillar points should we use? */ + if (a1[1] > b1[1]){ mask[0] = b1[1]; } else { mask[0] = a1[1]; } + if (a2[1] > b2[1]){ mask[2] = b2[1]; } else { mask[2] = a2[1]; } + if (a2[0] > b2[0]){ mask[4] = a2[0]; } else { mask[4] = b2[0]; } + if (a1[0] > b1[0]){ mask[6] = a1[0]; } else { mask[6] = b1[0]; } #if DEBUG - /* Illegal situations */ - if (mask [0] == mask[2] || - mask [0] == mask[4] || - mask [2] == mask[6] || - mask [4] == mask[6]){ - fprintf(stderr, "Illegal Partial pinch!\n"); - } + /* Illegal situations */ + if (mask [0] == mask[2] || + mask [0] == mask[4] || + mask [2] == mask[6] || + mask [4] == mask[6]){ + fprintf(stderr, "Illegal Partial pinch!\n"); + } #endif - /* Partial pinch of face */ - if (mask[0] == mask[6]){ - mask[6] = -1; - } - - - if (mask[2] == mask[4]){ - mask[4] = -1; - } - - - - /* Get shape of face: */ - /* each new intersection will be part of the new face, */ - /* but not all pillar points. This is encoded in mask. */ - - - mask[1] = intersect[3]; /* top-top */ - mask[3] = -1; - mask[5] = intersect[0]; /* bottom-bottom*/ - mask[7] = -1; - - /* bottom-top */ - if (intersect[1] != -1){ - if(a1[0] > b1[1]){ /* intersection[1] left of (any) intersection[0] */ - mask[0] = -1; - mask[6] = -1; - mask[7] = intersect[1]; + /* Partial pinch of face */ + if (mask[0] == mask[6]){ + mask[6] = -1; } - else{ - mask[2] = -1; - mask[4] = -1; - mask[3] = intersect[1]; - } - } - /* top-bottom */ - if (intersect[2] != -1){ - if(a1[1] < b1[0]){ /* intersection[2] left of (any) intersection[3] */ - mask[0] = -1; - mask[6] = -1; - mask[7] = intersect[2]; + + if (mask[2] == mask[4]){ + mask[4] = -1; } - else{ - mask[2] = -1; - mask[4] = -1; - mask[3] = intersect[2]; + + + + /* Get shape of face: */ + /* each new intersection will be part of the new face, */ + /* but not all pillar points. This is encoded in mask. */ + + + mask[1] = intersect[3]; /* top-top */ + mask[3] = -1; + mask[5] = intersect[0]; /* bottom-bottom*/ + mask[7] = -1; + + /* bottom-top */ + if (intersect[1] != -1){ + if(a1[0] > b1[1]){ /* intersection[1] left of (any) intersection[0] */ + mask[0] = -1; + mask[6] = -1; + mask[7] = intersect[1]; + } + else{ + mask[2] = -1; + mask[4] = -1; + mask[3] = intersect[1]; + } } - } - f = faces; - for (k=7; k>=0; --k){ - if(mask[k] != -1){ - *f++ = mask[k]; + + /* top-bottom */ + if (intersect[2] != -1){ + if(a1[1] < b1[0]){ /* intersection[2] left of (any) intersection[3] */ + mask[0] = -1; + mask[6] = -1; + mask[7] = intersect[2]; + } + else{ + mask[2] = -1; + mask[4] = -1; + mask[3] = intersect[2]; + } + } + f = faces; + for (k=7; k>=0; --k){ + if(mask[k] != -1){ + *f++ = mask[k]; + } } - } #if DEBUG>1 - /* Check for repeated nodes:*/ - int i; - fprintf(stderr, "face: "); - for (i=0; i<8; ++i){ - fprintf(stderr, "%d ", mask[i]); - for (k=0; k<8; ++k){ - if (i!=k && mask[i] != -1 && mask[i] == mask[k]){ - fprintf(stderr, "Repeated node in faulted face\n"); - } + /* Check for repeated nodes:*/ + int i; + fprintf(stderr, "face: "); + for (i=0; i<8; ++i){ + fprintf(stderr, "%d ", mask[i]); + for (k=0; k<8; ++k){ + if (i!=k && mask[i] != -1 && mask[i] == mask[k]){ + fprintf(stderr, "Repeated node in faulted face\n"); + } + } } - } - fprintf(stderr, "\n"); + fprintf(stderr, "\n"); #endif - - - return f; - + return f; + + + } /* a) If we assume that the index increase when z increase for each - pillar (but only separately), we can use only the point indices, - since we only need to compare z-values on one pillar at a time. + pillar (but only separately), we can use only the point indices, + since we only need to compare z-values on one pillar at a time. b) We assume input is preprocessed such that no intersections occur - on the first and last lines, for instance by padding the grid - with extra cells. This is convenient in the identification of - (unique) intersections. + on the first and last lines, for instance by padding the grid with + extra cells. This is convenient in the identification of (unique) + intersections. */ #define lineintersection(a1,a2,b1,b2)(((a1>b1)&&(a2b2))) static int faceintersection(int *a1, int *a2, int *b1, int *b2) { - return - max(a1[0],b1[0]) < min(a1[1],b1[1]) || - max(a2[0],b2[0]) < min(a2[1],b2[1]) || - lineintersection(a1[0], a2[0], b1[0], b2[0]) || - lineintersection(a1[1], a2[1], b1[1], b2[1]); + return + max(a1[0],b1[0]) < min(a1[1],b1[1]) || + max(a2[0],b2[0]) < min(a2[1],b2[1]) || + lineintersection(a1[0], a2[0], b1[0], b2[0]) || + lineintersection(a1[1], a2[1], b1[1], b2[1]); } @@ -197,156 +197,158 @@ static int faceintersection(int *a1, int *a2, int *b1, int *b2) /* work should be pointer to 2n ints initialised to zero . */ void findconnections(int n, int *pts[4], - int *intersectionlist, - int *work, - struct processed_grid *out) + int *intersectionlist, + int *work, + struct processed_grid *out) { - /* vectors of point numbers for faces a(b) on pillar 1(2) */ - int *a1 = pts[0]; - int *a2 = pts[1]; - int *b1 = pts[2]; - int *b2 = pts[3]; + /* vectors of point numbers for faces a(b) on pillar 1(2) */ + int *a1 = pts[0]; + int *a2 = pts[1]; + int *b1 = pts[2]; + int *b2 = pts[3]; - /* Intersection record for top line and bottomline of a */ - int *itop = work; - int *ibottom = work + n; - int *f = out->face_nodes + out->face_ptr[out->number_of_faces]; - int *c = out->face_neighbors + 2*out->number_of_faces; + /* Intersection record for top line and bottomline of a */ + int *itop = work; + int *ibottom = work + n; + int *f = out->face_nodes + out->face_ptr[out->number_of_faces]; + int *c = out->face_neighbors + 2*out->number_of_faces; - int k1 = 0; - int k2 = 0; + int k1 = 0; + int k2 = 0; - int i,j=0; - int intersect[4]= {-1}; - int *tmp; - /* for (i=0; i<2*n; work[i++]=-1); */ - - for (i = 0; iface_ptr[++out->number_of_faces] = f - out->face_nodes; + + } + else{ + ; + } + } + } + + /* Non-matching faces */ + else{ + + /* Find new intersection */ + if (lineintersection(a1[i+1],a2[i+1],b1[j+1],b2[j+1])) { + itop[j+1] = out->number_of_nodes++; + + /* store point numbers of intersecting lines */ + *intersectionlist++ = a1[i+1]; + *intersectionlist++ = a2[i+1]; + *intersectionlist++ = b1[j+1]; + *intersectionlist++ = b2[j+1]; + + + }else{ + itop[j+1] = -1; + } + + /* Update intersection record */ + intersect[0] = ibottom[j ]; /* i x j */ + intersect[1] = ibottom[j+1]; /* i x j+1 */ + intersect[2] = itop[j ]; /* i+1 x j */ + intersect[3] = itop[j+1]; /* i+1 x j+1 */ + + + /* Add face to list of faces if no INT_MIN or + * INT_MAX appear in a or b. */ + if (meaningful_face(i,j)){ + + /* + Even indices refer to space between cells, + odd indices refer to cells + */ + int cell_a = i%2 != 0 ? (i-1)/2 : -1; + int cell_b = j%2 != 0 ? (j-1)/2 : -1; + + + + if (cell_a != -1 || cell_b != -1){ + *c++ = cell_a; + *c++ = cell_b; + + f = computeFaceTopology(a1+i, a2+i, b1+j, b2+j, intersect, f); + + out->face_ptr[++out->number_of_faces] = f - out->face_nodes; + } + else{ + ; + } + } + } + } + + /* Update candidates for restart of j for in next i-iteration */ + if (b1[j] < a1[i+1]) { k1 = j; } + if (b2[j] < a2[i+1]) { k2 = j; } + + j = j+1; + } + + + + /* Swap intersection records: top line of a[i,i+1] is bottom + * line of a[i+1,i+2] */ + tmp = itop; itop = ibottom; ibottom = tmp; + + /* Zero out the "new" itop */ + for(j=0;jface_ptr[++out->number_of_faces] = f - out->face_nodes; - - } - else{ - ; - } - } - } - - /* Non-matching faces */ - else{ - - /* Find new intersection */ - if (lineintersection(a1[i+1],a2[i+1],b1[j+1],b2[j+1])) { - itop[j+1] = out->number_of_nodes++; - - /* store point numbers of intersecting lines */ - *intersectionlist++ = a1[i+1]; - *intersectionlist++ = a2[i+1]; - *intersectionlist++ = b1[j+1]; - *intersectionlist++ = b2[j+1]; - - - }else{ - itop[j+1] = -1; - } - - /* Update intersection record */ - intersect[0] = ibottom[j ]; /* i x j */ - intersect[1] = ibottom[j+1]; /* i x j+1 */ - intersect[2] = itop[j ]; /* i+1 x j */ - intersect[3] = itop[j+1]; /* i+1 x j+1 */ - - - /* Add face to list of faces if no INT_MIN or INT_MAX appear in a or b. */ - if (meaningful_face(i,j)){ - - /* - Even indices refer to space between cells, - odd indices refer to cells - */ - int cell_a = i%2 != 0 ? (i-1)/2 : -1; - int cell_b = j%2 != 0 ? (j-1)/2 : -1; - - - - if (cell_a != -1 || cell_b != -1){ - *c++ = cell_a; - *c++ = cell_b; - - f = computeFaceTopology(a1+i, a2+i, b1+j, b2+j, intersect, f); - - out->face_ptr[++out->number_of_faces] = f - out->face_nodes; - } - else{ - ; - } - } - } - } - - /* Update candidates for restart of j for in next i-iteration */ - if (b1[j] < a1[i+1]) { k1 = j; } - if (b2[j] < a2[i+1]) { k2 = j; } - - j = j+1; - } - - - - /* Swap intersection records: top line of a[i,i+1] is bottom line of a[i+1,i+2] */ - tmp = itop; itop = ibottom; ibottom = tmp; - - /* Zero out the "new" itop */ - for(j=0;j. + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ -#ifndef OPENRS_FACETOPOLOGY_HEADER -#define OPENRS_FACETOPOLOGY_HEADER +#ifndef OPM_FACETOPOLOGY_HEADER +#define OPM_FACETOPOLOGY_HEADER void findconnections(int n, int *pts[4], - int *intersectionlist, + int *intersectionlist, int *work, - struct processed_grid *out); + struct processed_grid *out); -#endif /* OPENRS_FACETOPOLOGY_HEADER */ +#endif /* OPM_FACETOPOLOGY_HEADER */ /* Local Variables: */ /* c-basic-offset:4 */ diff --git a/grdecl.h b/grdecl.h index e9a636ae..f9045f32 100644 --- a/grdecl.h +++ b/grdecl.h @@ -1,12 +1,16 @@ -#ifndef GRDECL_H -#define GRDECL_H +#ifndef GRDECL_H_INCLUDED +#define GRDECL_H_INCLUDED struct grdecl{ - int dims[3]; - const double *coord; - const double *zcorn; - const int *actnum; + int dims[3]; + const double *coord; + const double *zcorn; + const int *actnum; }; -#endif +#endif /* GRDECL_H_INCLUDED */ + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ diff --git a/mxgrdecl.c b/mxgrdecl.c index ab356ef7..fdb633e9 100644 --- a/mxgrdecl.c +++ b/mxgrdecl.c @@ -13,23 +13,23 @@ //=======================================================================*/ /* -Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. -Copyright 2009, 2010 Statoil ASA. + Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. + Copyright 2009, 2010 Statoil ASA. -This file is part of The Open Reservoir Simulator Project (OpenRS). + This file is part of The Open Reservoir Simulator Project (OpenRS). -OpenRS 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. + OpenRS 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. -OpenRS 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. + OpenRS 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 OpenRS. If not, see . + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ #include @@ -48,60 +48,64 @@ void mx_init_grdecl(struct grdecl *g, const mxArray *s); /*-------------------------------------------------------*/ void mx_init_grdecl(struct grdecl *g, const mxArray *s) { - int i,n; - mxArray *field; - int numel; - double *tmp; - mxArray *cartdims=NULL, *actnum=NULL, *coord=NULL, *zcorn=NULL; - - if (!mxIsStruct(s) - || !(cartdims = mxGetField(s, 0, "cartDims")) - || !(actnum = mxGetField(s, 0, "ACTNUM")) - || !(coord = mxGetField(s, 0, "COORD")) - || !(zcorn = mxGetField(s, 0, "ZCORN")) - ) - { - char str[]="Input must be a single Matlab struct with fields\n" - "cartDims, ACTNUM, COORD and ZCORN\n"; - mexErrMsgTxt(str); - } - - - numel = mxGetNumberOfElements(cartdims); - if (!mxIsNumeric(cartdims) || numel != 3){ - mexErrMsgTxt("cartDims field must be 3 numbers"); - } - - tmp = mxGetPr(cartdims); - n = 1; - for (i=0; i<3; ++i){ - g->dims[i] = tmp[i]; - n *= tmp[i]; - } + int i,n; + mxArray *field; + int numel; + double *tmp; + mxArray *cartdims=NULL, *actnum=NULL, *coord=NULL, *zcorn=NULL; + + if (!mxIsStruct(s) + || !(cartdims = mxGetField(s, 0, "cartDims")) + || !(actnum = mxGetField(s, 0, "ACTNUM")) + || !(coord = mxGetField(s, 0, "COORD")) + || !(zcorn = mxGetField(s, 0, "ZCORN")) + ) + { + char str[]="Input must be a single Matlab struct with fields\n" + "cartDims, ACTNUM, COORD and ZCORN\n"; + mexErrMsgTxt(str); + } - numel = mxGetNumberOfElements(actnum); - if (mxGetClassID(actnum) != mxINT32_CLASS || - numel != g->dims[0]*g->dims[1]*g->dims[2] ){ - mexErrMsgTxt("ACTNUM field must be nx*ny*nz numbers int32"); - } - g->actnum = mxGetData(actnum); + numel = mxGetNumberOfElements(cartdims); + if (!mxIsNumeric(cartdims) || numel != 3){ + mexErrMsgTxt("cartDims field must be 3 numbers"); + } + + tmp = mxGetPr(cartdims); + n = 1; + for (i=0; i<3; ++i){ + g->dims[i] = tmp[i]; + n *= tmp[i]; + } - - field = mxGetField(s, 0, "COORD"); - numel = mxGetNumberOfElements(coord); - if (mxGetClassID(coord) != mxDOUBLE_CLASS || - numel != 6*(g->dims[0]+1)*(g->dims[1]+1)){ - mexErrMsgTxt("COORD field must have 6*(nx+1)*(ny+1) doubles."); - } - g->coord = mxGetPr(coord); - + numel = mxGetNumberOfElements(actnum); + if (mxGetClassID(actnum) != mxINT32_CLASS || + numel != g->dims[0]*g->dims[1]*g->dims[2] ){ + mexErrMsgTxt("ACTNUM field must be nx*ny*nz numbers int32"); + } + g->actnum = mxGetData(actnum); - numel = mxGetNumberOfElements(zcorn); - if (mxGetClassID(zcorn) != mxDOUBLE_CLASS || - numel != 8*g->dims[0]*g->dims[1]*g->dims[2]){ - mexErrMsgTxt("ZCORN field must have 8*nx*ny*nz doubles."); - } - g->zcorn = mxGetPr(zcorn); + + + field = mxGetField(s, 0, "COORD"); + numel = mxGetNumberOfElements(coord); + if (mxGetClassID(coord) != mxDOUBLE_CLASS || + numel != 6*(g->dims[0]+1)*(g->dims[1]+1)){ + mexErrMsgTxt("COORD field must have 6*(nx+1)*(ny+1) doubles."); + } + g->coord = mxGetPr(coord); + + + numel = mxGetNumberOfElements(zcorn); + if (mxGetClassID(zcorn) != mxDOUBLE_CLASS || + numel != 8*g->dims[0]*g->dims[1]*g->dims[2]){ + mexErrMsgTxt("ZCORN field must have 8*nx*ny*nz doubles."); + } + g->zcorn = mxGetPr(zcorn); } + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ diff --git a/mxgrdecl.h b/mxgrdecl.h index e46c6973..4621c202 100644 --- a/mxgrdecl.h +++ b/mxgrdecl.h @@ -13,30 +13,32 @@ //=======================================================================*/ /* -Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. -Copyright 2009, 2010 Statoil ASA. + Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. + Copyright 2009, 2010 Statoil ASA. -This file is part of The Open Reservoir Simulator Project (OpenRS). + This file is part of The Open Reservoir Simulator Project (OpenRS). -OpenRS 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. + OpenRS 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. -OpenRS 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. + OpenRS 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 OpenRS. If not, see . + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ -#ifndef OPENRS_MXGRDECL_HEADER -#define OPENRS_MXGRDECL_HEADER +#ifndef OPM_MXGRDECL_HEADER +#define OPM_MXGRDECL_HEADER void mx_init_grdecl (struct grdecl *g, const mxArray *s); -#endif /* OPENRS_MXGRDECL_HEADER */ - +#endif /* OPM_MXGRDECL_HEADER */ +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ diff --git a/preprocess.c b/preprocess.c index 9a1d8608..9c8aa69f 100644 --- a/preprocess.c +++ b/preprocess.c @@ -13,23 +13,23 @@ //==========================================================================*/ /* -Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. -Copyright 2009, 2010 Statoil ASA. + Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. + Copyright 2009, 2010 Statoil ASA. -This file is part of The Open Reservoir Simulator Project (OpenRS). + This file is part of The Open Reservoir Simulator Project (OpenRS). -OpenRS 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. + OpenRS 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. -OpenRS 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. + OpenRS 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 OpenRS. If not, see . + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ #include @@ -49,31 +49,30 @@ along with OpenRS. If not, see . void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len); int checkmemory(int nz, struct processed_grid *out, int **intersections); void process_vertical_faces(int direction, - int **intersections, - int *plist, int *work, - struct processed_grid *out); + int **intersections, + int *plist, int *work, + struct processed_grid *out); void process_horizontal_faces(int **intersections, - int *plist, - struct processed_grid *out); + int *plist, + struct processed_grid *out); /*----------------------------------------------------------------- Given a vector with k index running faster than i running faster than j, and Cartesian dimensions , find pointers to the (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of - field. - */ + field. */ static void igetvectors(int dims[3], int i, int j, int *field, int *v[]) { - int im = max(1, i ) - 1; - int ip = min(dims[0], i+1) - 1; - int jm = max(1, j ) - 1; - int jp = min(dims[1], j+1) - 1; + int im = max(1, i ) - 1; + int ip = min(dims[0], i+1) - 1; + int jm = max(1, j ) - 1; + int jp = min(dims[1], j+1) - 1; - v[0] = field + dims[2]*(im + dims[0]* jm); - v[1] = field + dims[2]*(im + dims[0]* jp); - v[2] = field + dims[2]*(ip + dims[0]* jm); - v[3] = field + dims[2]*(ip + dims[0]* jp); + v[0] = field + dims[2]*(im + dims[0]* jm); + v[1] = field + dims[2]*(im + dims[0]* jp); + v[2] = field + dims[2]*(ip + dims[0]* jm); + v[3] = field + dims[2]*(ip + dims[0]* jp); } @@ -84,269 +83,268 @@ static void igetvectors(int dims[3], int i, int j, int *field, int *v[]) /*----------------------------------------------------------------- Special purpose - Convert from k-index to Cartesian index i+nx*(j+ny*k) for every other - element in neighbors. + Convert from k-index to Cartesian index i+nx*(j+ny*k) for every + other element in neighbors. - */ +*/ void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len) { - int k; - if (i<0 || i>=dims[0] || j<0 || j >= dims[1]){ - for(k=0; k=dims[0] || j<0 || j >= dims[1]){ + for(k=0; km; - int n = out->n; + /* Ensure there is enough space */ + int r = (2*nz+2)*(2*nz+2); + int m = out->m; + int n = out->n; - if(out->number_of_faces + r > m){ - m += max(m*0.5, 2*r); - } - if (out->face_ptr[out->number_of_faces] + 6*r > n){ - n += max(n*0.5, 12*r); - } - - if (m != out->m){ - void *p1 = realloc(out->face_neighbors, 2*m * sizeof(*out->face_neighbors)); - void *p2 = realloc(*intersections, 4*m * sizeof(**intersections)); - if (p1 && p2) { - out->face_neighbors = p1; - *intersections = p2; - } else { - return 0; + if(out->number_of_faces + r > m){ + m += max(m*0.5, 2*r); } - } - - if (m != out->m || n != out->n) { - void *p1 = realloc(out->face_nodes, n * sizeof *out->face_nodes); - void *p2 = realloc(out->face_ptr, (m+1) * sizeof *out->face_ptr); - void *p3 = realloc(out->face_tag, m * sizeof *out->face_tag); - - if (p1 && p2 && p3) { - out->face_nodes = p1; - out->face_ptr = p2; - out->face_tag = p3; - } else { - return 0; + if (out->face_ptr[out->number_of_faces] + 6*r > n){ + n += max(n*0.5, 12*r); } - out->m = m; - out->n = n; - } + if (m != out->m){ + void *p1 = realloc(out->face_neighbors, 2*m * sizeof(*out->face_neighbors)); + void *p2 = realloc(*intersections, 4*m * sizeof(**intersections)); + if (p1 && p2) { + out->face_neighbors = p1; + *intersections = p2; + } else { + return 0; + } + } - return 1; + if (m != out->m || n != out->n) { + void *p1 = realloc(out->face_nodes, n * sizeof *out->face_nodes); + void *p2 = realloc(out->face_ptr, (m+1) * sizeof *out->face_ptr); + void *p3 = realloc(out->face_tag, m * sizeof *out->face_tag); + + if (p1 && p2 && p3) { + out->face_nodes = p1; + out->face_ptr = p2; + out->face_tag = p3; + } else { + return 0; + } + + out->m = m; + out->n = n; + } + + return 1; } /*----------------------------------------------------------------- For each vertical face (i.e. i or j constant), - -find point numbers for the corners and - -cell neighbors. - -new points on faults defined by two intgersecting lines. + -find point numbers for the corners and + -cell neighbors. + -new points on faults defined by two intgersecting lines. direction == 0 : constant-i faces. direction == 1 : constant-j faces. - */ +*/ void process_vertical_faces(int direction, - int **intersections, - int *plist, int *work, - struct processed_grid *out) + int **intersections, + int *plist, int *work, + struct processed_grid *out) { - int i,j; - int *cornerpts[4]; - int d[3]; - int f; - enum face_tag tag[] = { LEFT, BACK }; - int *tmp; - int nx = out->dimensions[0]; - int ny = out->dimensions[1]; - int nz = out->dimensions[2]; - int startface; - int num_intersections; - int *ptr; - int len; - for (j=0; jdimensions[0]; + int ny = out->dimensions[1]; + int nz = out->dimensions[2]; + int startface; + int num_intersections; + int *ptr; + int len; + for (j=0; j */ - /* 0 2 2 3 */ - /* rotate clockwise */ - tmp = cornerpts[1]; - cornerpts[1] = cornerpts[0]; - cornerpts[0] = cornerpts[2]; - cornerpts[2] = cornerpts[3]; - cornerpts[3] = tmp; - } + if(direction==1){ + /* 1 3 0 1 */ + /* ---> */ + /* 0 2 2 3 */ + /* rotate clockwise */ + tmp = cornerpts[1]; + cornerpts[1] = cornerpts[0]; + cornerpts[0] = cornerpts[2]; + cornerpts[2] = cornerpts[3]; + cornerpts[3] = tmp; + } - /* int startface = ftab->position; */ - startface = out->number_of_faces; - /* int num_intersections = *npoints - npillarpoints; */ - num_intersections = out->number_of_nodes - - out->number_of_nodes_on_pillars; + /* int startface = ftab->position; */ + startface = out->number_of_faces; + /* int num_intersections = *npoints - npillarpoints; */ + num_intersections = out->number_of_nodes - + out->number_of_nodes_on_pillars; - findconnections(2*nz+2, cornerpts, - *intersections+4*num_intersections, - work, out); + findconnections(2*nz+2, cornerpts, + *intersections+4*num_intersections, + work, out); - ptr = out->face_neighbors + 2*startface; - len = 2*out->number_of_faces - 2*startface; + ptr = out->face_neighbors + 2*startface; + len = 2*out->number_of_faces - 2*startface; - compute_cell_index(out->dimensions, i-1+direction, j-direction, ptr, len); - compute_cell_index(out->dimensions, i, j, ptr+1, len); + compute_cell_index(out->dimensions, i-1+direction, j-direction, ptr, len); + compute_cell_index(out->dimensions, i, j, ptr+1, len); - /* Tag the new faces */ - f = startface; - for (; f < out->number_of_faces; ++f) { - out->face_tag[f] = tag[direction]; - } + /* Tag the new faces */ + f = startface; + for (; f < out->number_of_faces; ++f) { + out->face_tag[f] = tag[direction]; + } + } } - } } static int linearindex(const int dims[3], int i, int j, int k) { - return i+dims[0]*(j+dims[1]*k); + return i+dims[0]*(j+dims[1]*k); } /*----------------------------------------------------------------- For each horizontal face (i.e. k constant), - -find point numbers for the corners and - -cell neighbors. + -find point numbers for the corners and + -cell neighbors. Also define map from logically Cartesian cell index to local cell index 0, ..., #. Exclude cells that are have collapsed coordinates. (This includes cells with ACTNUM==0) - */ +*/ void process_horizontal_faces(int **intersections, - int *plist, - struct processed_grid *out) + int *plist, + struct processed_grid *out) { - int i,j,k; + int i,j,k; - int nx = out->dimensions[0]; - int ny = out->dimensions[1]; - int nz = out->dimensions[2]; + int nx = out->dimensions[0]; + int ny = out->dimensions[1]; + int nz = out->dimensions[2]; - int *cell = out->local_cell_index; - int cellno = 0; - int *f, *n, *c[4]; - int prevcell, thiscell; - int idx; + int *cell = out->local_cell_index; + int cellno = 0; + int *f, *n, *c[4]; + int prevcell, thiscell; + int idx; - /* dimensions of plist */ - int d[3]; - d[0] = 2*nx; - d[1] = 2*ny; - d[2] = 2+2*nz; + /* dimensions of plist */ + int d[3]; + d[0] = 2*nx; + d[1] = 2*ny; + d[2] = 2+2*nz; - for(j=0; jface_nodes + out->face_ptr[out->number_of_faces]; - n = out->face_neighbors + 2*out->number_of_faces; + f = out->face_nodes + out->face_ptr[out->number_of_faces]; + n = out->face_neighbors + 2*out->number_of_faces; - /* Vectors of point numbers */ - igetvectors(d, 2*i+1, 2*j+1, plist, c); + /* Vectors of point numbers */ + igetvectors(d, 2*i+1, 2*j+1, plist, c); - prevcell = -1; + prevcell = -1; - for (k = 1; kdimensions, i,j,(k-1)/2); - cell[idx] = -1; - } - } - else{ + /* If the pinch is a cell: */ + if (k%2){ + idx = linearindex(out->dimensions, i,j,(k-1)/2); + cell[idx] = -1; + } + } + else{ - if (k%2){ - /* Add face */ - *f++ = c[0][k]; - *f++ = c[2][k]; - *f++ = c[3][k]; - *f++ = c[1][k]; + if (k%2){ + /* Add face */ + *f++ = c[0][k]; + *f++ = c[2][k]; + *f++ = c[3][k]; + *f++ = c[1][k]; - out->face_tag[ out->number_of_faces] = TOP; - out->face_ptr[++out->number_of_faces] = f - out->face_nodes; + out->face_tag[ out->number_of_faces] = TOP; + out->face_ptr[++out->number_of_faces] = f - out->face_nodes; - thiscell = linearindex(out->dimensions, i,j,(k-1)/2); - *n++ = prevcell; - *n++ = prevcell = thiscell; + thiscell = linearindex(out->dimensions, i,j,(k-1)/2); + *n++ = prevcell; + *n++ = prevcell = thiscell; - cell[thiscell] = cellno++; + cell[thiscell] = cellno++; - } - else{ - if (prevcell != -1){ - /* Add face */ - *f++ = c[0][k]; - *f++ = c[2][k]; - *f++ = c[3][k]; - *f++ = c[1][k]; + } + else{ + if (prevcell != -1){ + /* Add face */ + *f++ = c[0][k]; + *f++ = c[2][k]; + *f++ = c[3][k]; + *f++ = c[1][k]; - out->face_tag[ out->number_of_faces] = TOP; - out->face_ptr[++out->number_of_faces] = f - out->face_nodes; + out->face_tag[ out->number_of_faces] = TOP; + out->face_ptr[++out->number_of_faces] = f - out->face_nodes; - *n++ = prevcell; - *n++ = prevcell = -1; - } - } - } - } + *n++ = prevcell; + *n++ = prevcell = -1; + } + } + } + } + } } - } - out->number_of_cells = cellno; + out->number_of_cells = cellno; } @@ -363,55 +361,55 @@ void process_horizontal_faces(int **intersections, static void approximate_intersection_pt(int *L, double *c, double *pt) { - double z0 = c[3*L[0]+2]; - double z1 = c[3*L[1]+2]; - double z2 = c[3*L[2]+2]; - double z3 = c[3*L[3]+2]; + double z0 = c[3*L[0]+2]; + double z1 = c[3*L[1]+2]; + double z2 = c[3*L[2]+2]; + double z3 = c[3*L[3]+2]; - double a = (z2-z0)/(z1-z0 - (z3-z2)); - if (isinf(a) || isnan(a)){ - a = 0; - } + double a = (z2-z0)/(z1-z0 - (z3-z2)); + if (isinf(a) || isnan(a)){ + a = 0; + } - pt[0] = c[3*L[0]+0]* (1.0-a) + c[3*L[1]+0]* a; - pt[1] = c[3*L[0]+1]* (1.0-a) + c[3*L[1]+1]* a; - pt[2] = c[3*L[0]+2]* (1.0-a) + c[3*L[1]+2]* a; + pt[0] = c[3*L[0]+0]* (1.0-a) + c[3*L[1]+0]* a; + pt[1] = c[3*L[0]+1]* (1.0-a) + c[3*L[1]+1]* a; + pt[2] = c[3*L[0]+2]* (1.0-a) + c[3*L[1]+2]* a; } /*----------------------------------------------------------------- - Compute x,y and z coordinates for points on each pillar. - Then, append x,y and z coordinates for extra points on faults. - */ + Compute x,y and z coordinates for points on each pillar. Then, + append x,y and z coordinates for extra points on faults. */ static void compute_intersection_coordinates(int *intersections, - struct processed_grid *out) + struct processed_grid *out) { - int n = out->number_of_nodes; - int np = out->number_of_nodes_on_pillars; - int k; - double *pt; - int *itsct = intersections; - /* Make sure the space allocated for nodes match the number of node. */ - void *p = realloc (out->node_coordinates, 3*n*sizeof(double)); - if (p) { - out->node_coordinates = p; - } - else { - fprintf(stderr, "Could not allocate extra space for intersections\n"); - } + int n = out->number_of_nodes; + int np = out->number_of_nodes_on_pillars; + int k; + double *pt; + int *itsct = intersections; + /* Make sure the space allocated for nodes match the number of + * node. */ + void *p = realloc (out->node_coordinates, 3*n*sizeof(double)); + if (p) { + out->node_coordinates = p; + } + else { + fprintf(stderr, "Could not allocate extra space for intersections\n"); + } - /* Append intersections */ - pt = out->node_coordinates + 3*np; + /* Append intersections */ + pt = out->node_coordinates + 3*np; - for (k=np; knode_coordinates, pt); - pt += 3; - itsct += 4; + for (k=np; knode_coordinates, pt); + pt += 3; + itsct += 4; - } + } } @@ -425,9 +423,9 @@ copy_and_permute_actnum(int nx, int ny, int nz, const int *in, int *out) /* Permute actnum such that values of each vertical stack of cells * are adjacent in memory, i.e., - out = [in(0,0,:), in(1,0,:),..., in(nx-1, ny-1,:)] + out = [in(0,0,:), in(1,0,:),..., in(nx-1, ny-1,:)] - in Matlab pseudo-code. + in Matlab pseudo-code. */ for (j=0; jdims[0]; - const int ny = in->dims[1]; - const int nz = in->dims[2]; - const int nc = nx*ny*nz; + const int BIGNUM = 64; + const int nx = in->dims[0]; + const int ny = in->dims[1]; + const int nz = in->dims[2]; + const int nc = nx*ny*nz; - /* internal work arrays */ - int *work; - int *plist; - int *intersections; + /* internal work arrays */ + int *work; + int *plist; + int *intersections; - /* -----------------------------------------------------------------*/ - /* Initialize output structure: - 1) allocate space for grid topology (which may need to be increased) - 2) set Cartesian imensions - */ - out->m = BIGNUM/3; - out->n = BIGNUM; + /* -----------------------------------------------------------------*/ + /* Initialize output structure: + 1) allocate space for grid topology (which may need to be + increased) + 2) set Cartesian imensions + */ + out->m = BIGNUM/3; + out->n = BIGNUM; - out->face_neighbors = malloc( BIGNUM * sizeof *out->face_neighbors); - out->face_nodes = malloc( out->n * sizeof *out->face_nodes); - out->face_ptr = malloc((out->m + 1) * sizeof *out->face_ptr); - out->face_tag = malloc( out->m * sizeof *out->face_tag); - out->face_ptr[0] = 0; + out->face_neighbors = malloc( BIGNUM * sizeof *out->face_neighbors); + out->face_nodes = malloc( out->n * sizeof *out->face_nodes); + out->face_ptr = malloc((out->m + 1) * sizeof *out->face_ptr); + out->face_tag = malloc( out->m * sizeof *out->face_tag); + out->face_ptr[0] = 0; - out->dimensions[0] = nx; - out->dimensions[1] = ny; - out->dimensions[2] = nz; - out->number_of_faces = 0; - out->number_of_nodes = 0; - out->number_of_cells = 0; + out->dimensions[0] = nx; + out->dimensions[1] = ny; + out->dimensions[2] = nz; + out->number_of_faces = 0; + out->number_of_nodes = 0; + out->number_of_cells = 0; - out->node_coordinates = NULL; - out->local_cell_index = malloc(nx*ny*nz * sizeof *out->local_cell_index); + out->node_coordinates = NULL; + out->local_cell_index = malloc(nx*ny*nz * sizeof *out->local_cell_index); - /* Do actual work here:*/ + /* Do actual work here:*/ - /* -----------------------------------------------------------------*/ - /* For each pillar, compare zcorn values for adjacent cells to find - * the unique node z-coordinates specified by the input. While - * here, enumerate unique points and assign point numbers (in plist) - * for each cornerpoint cell. In other words, plist has 8 node - * numbers for each cornerpoint cell.*/ + /* -----------------------------------------------------------------*/ + /* For each pillar, compare zcorn values for adjacent cells to + * find the unique node z-coordinates specified by the input. + * While here, enumerate unique points and assign point numbers + * (in plist) for each cornerpoint cell. In other words, plist has + * 8 node numbers for each cornerpoint cell.*/ - /* initialize grdecl structure "g" that will be processd by - * "finduniquepoints" */ - g.dims[0] = nx; - g.dims[1] = ny; - g.dims[2] = nz; + /* initialize grdecl structure "g" that will be processd by + * "finduniquepoints" */ + g.dims[0] = nx; + g.dims[1] = ny; + g.dims[2] = nz; - actnum = malloc (nc * sizeof *actnum); - g.actnum = copy_and_permute_actnum(nx, ny, nz, in->actnum, actnum); + actnum = malloc (nc * sizeof *actnum); + g.actnum = copy_and_permute_actnum(nx, ny, nz, in->actnum, actnum); - zcorn = malloc (nc * 8 * sizeof *zcorn); - sign = get_zcorn_sign(nx, ny, nz, in->actnum, in->zcorn, &error); - g.zcorn = copy_and_permute_zcorn(nx, ny, nz, in->zcorn, sign, zcorn); + zcorn = malloc (nc * 8 * sizeof *zcorn); + sign = get_zcorn_sign(nx, ny, nz, in->actnum, in->zcorn, &error); + g.zcorn = copy_and_permute_zcorn(nx, ny, nz, in->zcorn, sign, zcorn); - g.coord = in->coord; + g.coord = in->coord; - /* allocate space for cornerpoint numbers plus INT_MIN (INT_MAX) padding */ - plist = malloc( 4*nx*ny*(2*nz+2) * sizeof *plist); + /* allocate space for cornerpoint numbers plus INT_MIN (INT_MAX) + * padding */ + plist = malloc( 4*nx*ny*(2*nz+2) * sizeof *plist); - finduniquepoints(&g, plist, tolerance, out); + finduniquepoints(&g, plist, tolerance, out); - free (zcorn); - free (actnum); + free (zcorn); + free (actnum); - /* -----------------------------------------------------------------*/ - /* Find face topology and face-to-cell connections */ + /* -----------------------------------------------------------------*/ + /* Find face topology and face-to-cell connections */ - /* internal */ - work = malloc(2* (2*nz+2)* sizeof(*work)); - for(i=0; i<4*(nz+1); ++i) { work[i] = -1; } + /* internal */ + work = malloc(2* (2*nz+2)* sizeof(*work)); + for(i=0; i<4*(nz+1); ++i) { work[i] = -1; } - /* internal array to store intersections */ - intersections = malloc(BIGNUM* sizeof(*intersections)); + /* internal array to store intersections */ + intersections = malloc(BIGNUM* sizeof(*intersections)); - process_vertical_faces (0, &intersections, plist, work, out); - process_vertical_faces (1, &intersections, plist, work, out); - process_horizontal_faces ( &intersections, plist, out); + process_vertical_faces (0, &intersections, plist, work, out); + process_vertical_faces (1, &intersections, plist, work, out); + process_horizontal_faces ( &intersections, plist, out); - free (plist); - free (work); + free (plist); + free (work); - /* -----------------------------------------------------------------*/ - /* (re)allocate space for and compute coordinates of nodes that - * arise from intesecting cells (faults) */ - compute_intersection_coordinates(intersections, out); + /* -----------------------------------------------------------------*/ + /* (re)allocate space for and compute coordinates of nodes that + * arise from intesecting cells (faults) */ + compute_intersection_coordinates(intersections, out); - free (intersections); + free (intersections); - /* -----------------------------------------------------------------*/ - /* Enumerate compressed cells: - -make array [0...#cells-1] of global cell numbers - -make [0...nx*ny*nz-1] array of local cell numbers, - lexicographically ordered, used to remap out->face_neighbors - */ - global_cell_index = malloc(out->number_of_cells * - sizeof (*global_cell_index)); - cellnum = 0; - for (i=0; ilocal_cell_index[i]!=-1){ - global_cell_index[cellnum] = i; - out->local_cell_index[i] = cellnum; - cellnum++; + /* -----------------------------------------------------------------*/ + /* Enumerate compressed cells: + -make array [0...#cells-1] of global cell numbers + -make [0...nx*ny*nz-1] array of local cell numbers, + lexicographically ordered, used to remap out->face_neighbors + */ + global_cell_index = malloc(out->number_of_cells * + sizeof (*global_cell_index)); + cellnum = 0; + for (i=0; ilocal_cell_index[i]!=-1){ + global_cell_index[cellnum] = i; + out->local_cell_index[i] = cellnum; + cellnum++; + } } - } - /* Remap out->face_neighbors */ - iptr=out->face_neighbors; - for (i=0; inumber_of_faces*2; ++i, ++iptr){ - if (*iptr != -1){ - *iptr = out->local_cell_index[*iptr]; + /* Remap out->face_neighbors */ + iptr=out->face_neighbors; + for (i=0; inumber_of_faces*2; ++i, ++iptr){ + if (*iptr != -1){ + *iptr = out->local_cell_index[*iptr]; + } } - } - - free(out->local_cell_index); - out->local_cell_index = global_cell_index; - /* if sign==-1 in ZCORN preprocessing, the sign of the z-coordinate - * need to change before we finish */ - if (sign == -1) - { - for (i=2; i<3*out->number_of_nodes; i = i+3) - out->node_coordinates[i]=sign*out->node_coordinates[i]; - } + free(out->local_cell_index); + out->local_cell_index = global_cell_index; + + /* if sign==-1 in ZCORN preprocessing, the sign of the + * z-coordinate need to change before we finish */ + if (sign == -1) + { + for (i=2; i<3*out->number_of_nodes; i = i+3) + out->node_coordinates[i]=sign*out->node_coordinates[i]; + } } /*-------------------------------------------------------*/ void free_processed_grid(struct processed_grid *g) { - if( g ){ - free ( g->face_nodes ); - free ( g->face_ptr ); - free ( g->face_tag ); - free ( g->face_neighbors ); - free ( g->node_coordinates ); - free ( g->local_cell_index ); - } + if( g ){ + free ( g->face_nodes ); + free ( g->face_ptr ); + free ( g->face_tag ); + free ( g->face_neighbors ); + free ( g->node_coordinates ); + free ( g->local_cell_index ); + } } /* Local Variables: */ diff --git a/processgrid.c b/processgrid.c index 03bc0942..0cbfb814 100644 --- a/processgrid.c +++ b/processgrid.c @@ -13,23 +13,23 @@ //===========================================================================*/ /* -Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. -Copyright 2009, 2010, 2011, 2012 Statoil ASA. + Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. + Copyright 2009, 2010, 2011, 2012 Statoil ASA. -This file is part of The Open Reservoir Simulator Project (OpenRS). + This file is part of The Open Reservoir Simulator Project (OpenRS). -OpenRS 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. + OpenRS 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. -OpenRS 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. + OpenRS 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 OpenRS. If not, see . + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ /* Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. */ @@ -51,30 +51,30 @@ static mxArray * allocate_nodes(size_t nnodes) /* ---------------------------------------------------------------------- */ { - size_t nflds; - const char *fields[] = { "num", "coords" }; + size_t nflds; + const char *fields[] = { "num", "coords" }; - mxArray *nodes, *num, *coords; + mxArray *nodes, *num, *coords; - nflds = sizeof(fields) / sizeof(fields[0]); + nflds = sizeof(fields) / sizeof(fields[0]); - nodes = mxCreateStructMatrix(1, 1, nflds, fields); + nodes = mxCreateStructMatrix(1, 1, nflds, fields); - num = mxCreateDoubleScalar(nnodes); - coords = mxCreateDoubleMatrix(nnodes, 3, mxREAL); + num = mxCreateDoubleScalar(nnodes); + coords = mxCreateDoubleMatrix(nnodes, 3, mxREAL); - if ((nodes != NULL) && (num != NULL) && (coords != NULL)) { - mxSetField(nodes, 0, "num" , num); - mxSetField(nodes, 0, "coords", coords); - } else { - if (coords != NULL) { mxDestroyArray(coords); } - if (num != NULL) { mxDestroyArray(num); } - if (nodes != NULL) { mxDestroyArray(nodes); } + if ((nodes != NULL) && (num != NULL) && (coords != NULL)) { + mxSetField(nodes, 0, "num" , num); + mxSetField(nodes, 0, "coords", coords); + } else { + if (coords != NULL) { mxDestroyArray(coords); } + if (num != NULL) { mxDestroyArray(num); } + if (nodes != NULL) { mxDestroyArray(nodes); } - nodes = NULL; - } + nodes = NULL; + } - return nodes; + return nodes; } @@ -83,17 +83,17 @@ static void fill_nodes(mxArray *nodes, struct processed_grid *grid) /* ---------------------------------------------------------------------- */ { - size_t i, j, nnodes; - double *coords; + size_t i, j, nnodes; + double *coords; - nnodes = grid->number_of_nodes; - coords = mxGetPr(mxGetField(nodes, 0, "coords")); + nnodes = grid->number_of_nodes; + coords = mxGetPr(mxGetField(nodes, 0, "coords")); - for (j = 0; j < 3; ++j) { - for (i = 0; i < nnodes; ++i) { - *coords++ = grid->node_coordinates[3*i + j]; + for (j = 0; j < 3; ++j) { + for (i = 0; i < nnodes; ++i) { + *coords++ = grid->node_coordinates[3*i + j]; + } } - } } @@ -102,40 +102,40 @@ static mxArray * allocate_faces(size_t nf, size_t nfacenodes) /* ---------------------------------------------------------------------- */ { - size_t nflds; - const char *fields[] = { "num", "neighbors", "nodePos", "nodes", "tag" }; + size_t nflds; + const char *fields[] = { "num", "neighbors", "nodePos", "nodes", "tag" }; - mxArray *faces, *num, *neighbors, *nodePos, *nodes, *tag; + mxArray *faces, *num, *neighbors, *nodePos, *nodes, *tag; - nflds = sizeof(fields) / sizeof(fields[0]); + nflds = sizeof(fields) / sizeof(fields[0]); - faces = mxCreateStructMatrix(1, 1, nflds, fields); + faces = mxCreateStructMatrix(1, 1, nflds, fields); - num = mxCreateDoubleScalar (nf); - neighbors = mxCreateNumericMatrix(nf , 2, mxINT32_CLASS, mxREAL); - nodePos = mxCreateNumericMatrix(nf + 1 , 1, mxINT32_CLASS, mxREAL); - nodes = mxCreateNumericMatrix(nfacenodes, 1, mxINT32_CLASS, mxREAL); - tag = mxCreateNumericMatrix(nf , 1, mxINT32_CLASS, mxREAL); + num = mxCreateDoubleScalar (nf); + neighbors = mxCreateNumericMatrix(nf , 2, mxINT32_CLASS, mxREAL); + nodePos = mxCreateNumericMatrix(nf + 1 , 1, mxINT32_CLASS, mxREAL); + nodes = mxCreateNumericMatrix(nfacenodes, 1, mxINT32_CLASS, mxREAL); + tag = mxCreateNumericMatrix(nf , 1, mxINT32_CLASS, mxREAL); - if ((faces != NULL) && (num != NULL) && (neighbors != NULL) && - (nodePos != NULL) && (nodes != NULL) && (tag != NULL)) { - mxSetField(faces, 0, "num" , num); - mxSetField(faces, 0, "neighbors", neighbors); - mxSetField(faces, 0, "nodePos" , nodePos); - mxSetField(faces, 0, "nodes" , nodes); - mxSetField(faces, 0, "tag" , tag); - } else { - if (tag != NULL) { mxDestroyArray(tag); } - if (nodes != NULL) { mxDestroyArray(nodes); } - if (nodePos != NULL) { mxDestroyArray(nodePos); } - if (neighbors != NULL) { mxDestroyArray(neighbors); } - if (num != NULL) { mxDestroyArray(num); } - if (faces != NULL) { mxDestroyArray(faces); } + if ((faces != NULL) && (num != NULL) && (neighbors != NULL) && + (nodePos != NULL) && (nodes != NULL) && (tag != NULL)) { + mxSetField(faces, 0, "num" , num); + mxSetField(faces, 0, "neighbors", neighbors); + mxSetField(faces, 0, "nodePos" , nodePos); + mxSetField(faces, 0, "nodes" , nodes); + mxSetField(faces, 0, "tag" , tag); + } else { + if (tag != NULL) { mxDestroyArray(tag); } + if (nodes != NULL) { mxDestroyArray(nodes); } + if (nodePos != NULL) { mxDestroyArray(nodePos); } + if (neighbors != NULL) { mxDestroyArray(neighbors); } + if (num != NULL) { mxDestroyArray(num); } + if (faces != NULL) { mxDestroyArray(faces); } - faces = NULL; - } + faces = NULL; + } - return faces; + return faces; } @@ -144,33 +144,33 @@ static void fill_faces(mxArray *faces, struct processed_grid *grid) /* ---------------------------------------------------------------------- */ { - size_t i, f, nf, nfn; + size_t i, f, nf, nfn; - int *pi; + int *pi; - nf = grid->number_of_faces; + nf = grid->number_of_faces; - /* Fill faces.neighbors */ - pi = mxGetData(mxGetField(faces, 0, "neighbors")); - for (i = 0; i < 2; i++) { - for (f = 0; f < nf; f++) { - /* Add one for one-based indexing in M */ - *pi++ = grid->face_neighbors[2*f + i] + 1; + /* Fill faces.neighbors */ + pi = mxGetData(mxGetField(faces, 0, "neighbors")); + for (i = 0; i < 2; i++) { + for (f = 0; f < nf; f++) { + /* Add one for one-based indexing in M */ + *pi++ = grid->face_neighbors[2*f + i] + 1; + } } - } - /* Fill faces.nodePos */ - pi = mxGetData(mxGetField(faces, 0, "nodePos")); - for (i = 0; i <= nf; i++) { pi[i] = grid->face_ptr[i] + 1; } + /* Fill faces.nodePos */ + pi = mxGetData(mxGetField(faces, 0, "nodePos")); + for (i = 0; i <= nf; i++) { pi[i] = grid->face_ptr[i] + 1; } - /* Fill faces.nodes */ - pi = mxGetData(mxGetField(faces, 0, "nodes")); - nfn = grid->face_ptr[nf]; /* Total number of face nodes */ - for (i = 0; i < nfn; i++) { pi[i] = grid->face_nodes[i] + 1; } + /* Fill faces.nodes */ + pi = mxGetData(mxGetField(faces, 0, "nodes")); + nfn = grid->face_ptr[nf]; /* Total number of face nodes */ + for (i = 0; i < nfn; i++) { pi[i] = grid->face_nodes[i] + 1; } - /* Fill faces.tag */ - pi = mxGetData(mxGetField(faces, 0, "tag")); - for (f = 0; f < nf; f++) { pi[f] = grid->face_tag[f] + 1; } + /* Fill faces.tag */ + pi = mxGetData(mxGetField(faces, 0, "tag")); + for (f = 0; f < nf; f++) { pi[f] = grid->face_tag[f] + 1; } } @@ -179,18 +179,18 @@ static size_t count_halffaces(size_t nf, const int *neighbors) /* ---------------------------------------------------------------------- */ { - int c1, c2; - size_t nhf, f; + int c1, c2; + size_t nhf, f; - for (f = nhf = 0; f < nf; f++) { - c1 = neighbors[2*f + 0]; - c2 = neighbors[2*f + 1]; + for (f = nhf = 0; f < nf; f++) { + c1 = neighbors[2*f + 0]; + c2 = neighbors[2*f + 1]; - nhf += c1 >= 0; - nhf += c2 >= 0; - } + nhf += c1 >= 0; + nhf += c2 >= 0; + } - return nhf; + return nhf; } @@ -199,37 +199,37 @@ static mxArray * allocate_cells(size_t nc, size_t ncf) /* ---------------------------------------------------------------------- */ { - size_t nflds; - const char *fields[] = { "num", "facePos", "faces", "indexMap" }; + size_t nflds; + const char *fields[] = { "num", "facePos", "faces", "indexMap" }; - mxArray *cells, *num, *facePos, *faces, *indexMap; + mxArray *cells, *num, *facePos, *faces, *indexMap; - nflds = sizeof(fields) / sizeof(fields[0]); + nflds = sizeof(fields) / sizeof(fields[0]); - cells = mxCreateStructMatrix(1, 1, nflds, fields); + cells = mxCreateStructMatrix(1, 1, nflds, fields); - num = mxCreateDoubleScalar (nc); - facePos = mxCreateNumericMatrix(nc + 1, 1, mxINT32_CLASS, mxREAL); - faces = mxCreateNumericMatrix(ncf , 2, mxINT32_CLASS, mxREAL); - indexMap = mxCreateNumericMatrix(nc , 1, mxINT32_CLASS, mxREAL); + num = mxCreateDoubleScalar (nc); + facePos = mxCreateNumericMatrix(nc + 1, 1, mxINT32_CLASS, mxREAL); + faces = mxCreateNumericMatrix(ncf , 2, mxINT32_CLASS, mxREAL); + indexMap = mxCreateNumericMatrix(nc , 1, mxINT32_CLASS, mxREAL); - if ((cells != NULL) && (num != NULL) && (facePos != NULL) && - (faces != NULL) && (indexMap != NULL)) { - mxSetField(cells, 0, "num" , num ); - mxSetField(cells, 0, "facePos" , facePos ); - mxSetField(cells, 0, "faces" , faces ); - mxSetField(cells, 0, "indexMap", indexMap); - } else { - if (indexMap != NULL) { mxDestroyArray(indexMap); } - if (faces != NULL) { mxDestroyArray(faces); } - if (facePos != NULL) { mxDestroyArray(facePos); } - if (num != NULL) { mxDestroyArray(num); } - if (cells != NULL) { mxDestroyArray(cells); } + if ((cells != NULL) && (num != NULL) && (facePos != NULL) && + (faces != NULL) && (indexMap != NULL)) { + mxSetField(cells, 0, "num" , num ); + mxSetField(cells, 0, "facePos" , facePos ); + mxSetField(cells, 0, "faces" , faces ); + mxSetField(cells, 0, "indexMap", indexMap); + } else { + if (indexMap != NULL) { mxDestroyArray(indexMap); } + if (faces != NULL) { mxDestroyArray(faces); } + if (facePos != NULL) { mxDestroyArray(facePos); } + if (num != NULL) { mxDestroyArray(num); } + if (cells != NULL) { mxDestroyArray(cells); } - cells = NULL; - } + cells = NULL; + } - return cells; + return cells; } @@ -238,69 +238,69 @@ static void fill_cells(mxArray *cells, struct processed_grid *grid) /* ---------------------------------------------------------------------- */ { - size_t c, nc, f, nf, i; + size_t c, nc, f, nf, i; - int c1, c2, cf_tag, nhf; - int *pi1, *pi2; + int c1, c2, cf_tag, nhf; + int *pi1, *pi2; - nc = grid->number_of_cells; - nf = grid->number_of_faces; + nc = grid->number_of_cells; + nf = grid->number_of_faces; - /* Simultaneously fill cells.facePos and cells.faces by transposing the - * neighbours mapping. */ - pi1 = mxGetData(mxGetField(cells, 0, "facePos")); - pi2 = mxGetData(mxGetField(cells, 0, "faces" )); - for (i = 0; i < nc + 1; i++) { pi1[i] = 0; } + /* Simultaneously fill cells.facePos and cells.faces by transposing the + * neighbours mapping. */ + pi1 = mxGetData(mxGetField(cells, 0, "facePos")); + pi2 = mxGetData(mxGetField(cells, 0, "faces" )); + for (i = 0; i < nc + 1; i++) { pi1[i] = 0; } - /* 1) Count connections (i.e., faces per cell). */ - for (f = 0; f < nf; f++) { - c1 = grid->face_neighbors[2*f + 0]; - c2 = grid->face_neighbors[2*f + 1]; + /* 1) Count connections (i.e., faces per cell). */ + for (f = 0; f < nf; f++) { + c1 = grid->face_neighbors[2*f + 0]; + c2 = grid->face_neighbors[2*f + 1]; - if (c1 >= 0) { pi1[c1 + 1] += 1; } - if (c2 >= 0) { pi1[c2 + 1] += 1; } - } - - /* 2) Define start pointers (really, position *end* pointers at start). */ - for (c = 1; c <= nc; c++) { - pi1[0] += pi1[c]; - pi1[c] = pi1[0] - pi1[c]; - } - - /* 3) Fill connection structure whilst advancing end pointers. */ - nhf = pi1[0]; - pi1[0] = 0; - - mxAssert (((size_t) nhf) == mxGetM(mxGetField(cells, 0, "faces")), - "Number of half faces (SIZE(cells.faces,1)) incorrectly " - "determined earlier."); - - for (f = 0; f < nf; f++) { - cf_tag = 2*grid->face_tag[f] + 1; /* [1, 3, 5] */ - c1 = grid->face_neighbors[2*f + 0]; - c2 = grid->face_neighbors[2*f + 1]; - - if (c1 >= 0) { - pi2[ pi1[ c1 + 1 ] + 0*nhf ] = f + 1; - pi2[ pi1[ c1 + 1 ] + 1*nhf ] = cf_tag + 1; /* out */ - - pi1[ c1 + 1 ] += 1; + if (c1 >= 0) { pi1[c1 + 1] += 1; } + if (c2 >= 0) { pi1[c2 + 1] += 1; } } - if (c2 >= 0) { - pi2[ pi1[ c2 + 1 ] + 0*nhf ] = f + 1; - pi2[ pi1[ c2 + 1 ] + 1*nhf ] = cf_tag + 0; /* in */ - pi1[ c2 + 1 ] += 1; + /* 2) Define start pointers (really, position *end* pointers at start). */ + for (c = 1; c <= nc; c++) { + pi1[0] += pi1[c]; + pi1[c] = pi1[0] - pi1[c]; } - } - /* Finally, adjust pointer array for one-based indexing in M. */ - for (i = 0; i < nc + 1; i++) { pi1[i] += 1; } + /* 3) Fill connection structure whilst advancing end pointers. */ + nhf = pi1[0]; + pi1[0] = 0; - /* Fill cells.indexMap. Note that despite the name, 'local_cell_index' - * really *is* the (zero-based) indexMap of the 'processed_grid'. */ - pi1 = mxGetData(mxGetField(cells, 0, "indexMap")); - for (c = 0; c < nc; c++) { pi1[c] = grid->local_cell_index[c] + 1; } + mxAssert (((size_t) nhf) == mxGetM(mxGetField(cells, 0, "faces")), + "Number of half faces (SIZE(cells.faces,1)) incorrectly " + "determined earlier."); + + for (f = 0; f < nf; f++) { + cf_tag = 2*grid->face_tag[f] + 1; /* [1, 3, 5] */ + c1 = grid->face_neighbors[2*f + 0]; + c2 = grid->face_neighbors[2*f + 1]; + + if (c1 >= 0) { + pi2[ pi1[ c1 + 1 ] + 0*nhf ] = f + 1; + pi2[ pi1[ c1 + 1 ] + 1*nhf ] = cf_tag + 1; /* out */ + + pi1[ c1 + 1 ] += 1; + } + if (c2 >= 0) { + pi2[ pi1[ c2 + 1 ] + 0*nhf ] = f + 1; + pi2[ pi1[ c2 + 1 ] + 1*nhf ] = cf_tag + 0; /* in */ + + pi1[ c2 + 1 ] += 1; + } + } + + /* Finally, adjust pointer array for one-based indexing in M. */ + for (i = 0; i < nc + 1; i++) { pi1[i] += 1; } + + /* Fill cells.indexMap. Note that despite the name, 'local_cell_index' + * really *is* the (zero-based) indexMap of the 'processed_grid'. */ + pi1 = mxGetData(mxGetField(cells, 0, "indexMap")); + for (c = 0; c < nc; c++) { pi1[c] = grid->local_cell_index[c] + 1; } } @@ -309,52 +309,52 @@ static mxArray * allocate_grid(struct processed_grid *grid, const char *func) /* ---------------------------------------------------------------------- */ { - size_t nflds, nhf; - const char *fields[] = { "nodes", "faces", "cells", - "type", "cartDims", "griddim" }; + size_t nflds, nhf; + const char *fields[] = { "nodes", "faces", "cells", + "type", "cartDims", "griddim" }; - mxArray *G, *nodes, *faces, *cells; - mxArray *type, *typestr, *cartDims, *griddim; + mxArray *G, *nodes, *faces, *cells; + mxArray *type, *typestr, *cartDims, *griddim; - nflds = sizeof(fields) / sizeof(fields[0]); - nhf = count_halffaces(grid->number_of_faces, grid->face_neighbors); + nflds = sizeof(fields) / sizeof(fields[0]); + nhf = count_halffaces(grid->number_of_faces, grid->face_neighbors); - G = mxCreateStructMatrix(1, 1, nflds, fields); + G = mxCreateStructMatrix(1, 1, nflds, fields); - nodes = allocate_nodes(grid->number_of_nodes); - faces = allocate_faces(grid->number_of_faces, - grid->face_ptr[ grid->number_of_faces ]); - cells = allocate_cells(grid->number_of_cells, nhf); - type = mxCreateCellMatrix(1, 1); - typestr = mxCreateString(func); - cartDims = mxCreateDoubleMatrix(1, 3, mxREAL); - griddim = mxCreateDoubleScalar(3); + nodes = allocate_nodes(grid->number_of_nodes); + faces = allocate_faces(grid->number_of_faces, + grid->face_ptr[ grid->number_of_faces ]); + cells = allocate_cells(grid->number_of_cells, nhf); + type = mxCreateCellMatrix(1, 1); + typestr = mxCreateString(func); + cartDims = mxCreateDoubleMatrix(1, 3, mxREAL); + griddim = mxCreateDoubleScalar(3); - if ((G != NULL) && (nodes != NULL) && (faces != NULL) && - (cells != NULL) && (type != NULL) && (typestr != NULL) && - (cartDims != NULL) && (griddim != NULL)) { - mxSetCell(type, 0, typestr); + if ((G != NULL) && (nodes != NULL) && (faces != NULL) && + (cells != NULL) && (type != NULL) && (typestr != NULL) && + (cartDims != NULL) && (griddim != NULL)) { + mxSetCell(type, 0, typestr); - mxSetField(G, 0, "nodes" , nodes ); - mxSetField(G, 0, "faces" , faces ); - mxSetField(G, 0, "cells" , cells ); - mxSetField(G, 0, "type" , type ); - mxSetField(G, 0, "cartDims", cartDims); - mxSetField(G, 0, "griddim" , griddim ); - } else { - if (griddim != NULL) { mxDestroyArray(griddim); } - if (cartDims != NULL) { mxDestroyArray(cartDims); } - if (typestr != NULL) { mxDestroyArray(typestr); } - if (type != NULL) { mxDestroyArray(type); } - if (cells != NULL) { mxDestroyArray(cells); } - if (faces != NULL) { mxDestroyArray(faces); } - if (nodes != NULL) { mxDestroyArray(nodes); } - if (G != NULL) { mxDestroyArray(G); } + mxSetField(G, 0, "nodes" , nodes ); + mxSetField(G, 0, "faces" , faces ); + mxSetField(G, 0, "cells" , cells ); + mxSetField(G, 0, "type" , type ); + mxSetField(G, 0, "cartDims", cartDims); + mxSetField(G, 0, "griddim" , griddim ); + } else { + if (griddim != NULL) { mxDestroyArray(griddim); } + if (cartDims != NULL) { mxDestroyArray(cartDims); } + if (typestr != NULL) { mxDestroyArray(typestr); } + if (type != NULL) { mxDestroyArray(type); } + if (cells != NULL) { mxDestroyArray(cells); } + if (faces != NULL) { mxDestroyArray(faces); } + if (nodes != NULL) { mxDestroyArray(nodes); } + if (G != NULL) { mxDestroyArray(G); } - G = NULL; - } + G = NULL; + } - return G; + return G; } @@ -363,16 +363,16 @@ static void fill_grid(mxArray *G, struct processed_grid *grid) /* ---------------------------------------------------------------------- */ { - double *pr; + double *pr; - pr = mxGetPr(mxGetField(G, 0, "cartDims")); - pr[0] = grid->dimensions[0]; - pr[1] = grid->dimensions[1]; - pr[2] = grid->dimensions[2]; + pr = mxGetPr(mxGetField(G, 0, "cartDims")); + pr[0] = grid->dimensions[0]; + pr[1] = grid->dimensions[1]; + pr[2] = grid->dimensions[2]; - fill_nodes(mxGetField(G, 0, "nodes"), grid); - fill_faces(mxGetField(G, 0, "faces"), grid); - fill_cells(mxGetField(G, 0, "cells"), grid); + fill_nodes(mxGetField(G, 0, "nodes"), grid); + fill_faces(mxGetField(G, 0, "faces"), grid); + fill_cells(mxGetField(G, 0, "cells"), grid); } @@ -381,23 +381,23 @@ static int args_ok(int nlhs, int nrhs, const mxArray *prhs[]) /* ---------------------------------------------------------------------- */ { - int ok; + int ok; - ok = (nlhs == 1) && ((nrhs == 1) || (nrhs == 2)); + ok = (nlhs == 1) && ((nrhs == 1) || (nrhs == 2)); - ok = ok && !mxIsEmpty(prhs[0]); - ok = ok && mxIsStruct(prhs[0]); + ok = ok && !mxIsEmpty(prhs[0]); + ok = ok && mxIsStruct(prhs[0]); - ok = ok && (mxGetFieldNumber(prhs[0], "cartDims") >= 0); - ok = ok && (mxGetFieldNumber(prhs[0], "COORD" ) >= 0); - ok = ok && (mxGetFieldNumber(prhs[0], "ZCORN" ) >= 0); - ok = ok && (mxGetFieldNumber(prhs[0], "ACTNUM" ) >= 0); + ok = ok && (mxGetFieldNumber(prhs[0], "cartDims") >= 0); + ok = ok && (mxGetFieldNumber(prhs[0], "COORD" ) >= 0); + ok = ok && (mxGetFieldNumber(prhs[0], "ZCORN" ) >= 0); + ok = ok && (mxGetFieldNumber(prhs[0], "ACTNUM" ) >= 0); - if (ok && (nrhs == 2)) { - ok = mxIsDouble(prhs[1]) && (mxGetNumberOfElements(prhs[1]) == 1); - } + if (ok && (nrhs == 2)) { + ok = mxIsDouble(prhs[1]) && (mxGetNumberOfElements(prhs[1]) == 1); + } - return ok; + return ok; } @@ -406,56 +406,60 @@ static double define_tolerance(int nrhs, const mxArray *prhs[]) /* ---------------------------------------------------------------------- */ { - double tol; + double tol; - tol = 0.0; + tol = 0.0; - if (nrhs == 2) { - tol = mxGetScalar(prhs[1]); - } + if (nrhs == 2) { + tol = mxGetScalar(prhs[1]); + } - return tol; + return tol; } /* G = processgrid(grdecl) G = processgrid(grdecl, tolerance) - */ +*/ /* ---------------------------------------------------------------------- */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* ---------------------------------------------------------------------- */ { - double tolerance; - char errmsg[1023 + 1]; - struct grdecl grdecl; - struct processed_grid g; + double tolerance; + char errmsg[1023 + 1]; + struct grdecl grdecl; + struct processed_grid g; - if (args_ok(nlhs, nrhs, prhs)) { - mx_init_grdecl(&grdecl, prhs[0]); - tolerance = define_tolerance(nrhs, prhs); + if (args_ok(nlhs, nrhs, prhs)) { + mx_init_grdecl(&grdecl, prhs[0]); + tolerance = define_tolerance(nrhs, prhs); - process_grdecl(&grdecl, tolerance, &g); + process_grdecl(&grdecl, tolerance, &g); - plhs[0] = allocate_grid(&g, mexFunctionName()); + plhs[0] = allocate_grid(&g, mexFunctionName()); - if (plhs[0] != NULL) { - fill_grid(plhs[0], &g); + if (plhs[0] != NULL) { + fill_grid(plhs[0], &g); + } else { + /* Failed to create grid structure. Return empty. */ + plhs[0] = mxCreateDoubleMatrix(0, 0, mxREAL); + } + + free_processed_grid(&g); } else { - /* Failed to create grid structure. Return empty. */ - plhs[0] = mxCreateDoubleMatrix(0, 0, mxREAL); + sprintf(errmsg, + "Calling sequence is\n\t" + "G = %s(grdecl)\t%%or\n\t" + "G = %s(grdecl, tolerance)\n" + "The 'grdecl' must be a valid structure with fields\n" + "\t'cartDims', 'COORD', 'ZCORN', and 'ACTNUM'", + mexFunctionName(), mexFunctionName()); + mexErrMsgTxt(errmsg); } - - free_processed_grid(&g); - } else { - sprintf(errmsg, - "Calling sequence is\n\t" - "G = %s(grdecl)\t%%or\n\t" - "G = %s(grdecl, tolerance)\n" - "The 'grdecl' must be a valid structure with fields\n" - "\t'cartDims', 'COORD', 'ZCORN', and 'ACTNUM'", - mexFunctionName(), mexFunctionName()); - mexErrMsgTxt(errmsg); - } } + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ diff --git a/uniquepoints.c b/uniquepoints.c index bf589c26..628c179e 100644 --- a/uniquepoints.c +++ b/uniquepoints.c @@ -20,7 +20,7 @@ OpenRS 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 + the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OpenRS is distributed in the hope that it will be useful, @@ -50,309 +50,306 @@ /*----------------------------------------------------------------- - Compare function passed to qsort -*/ + Compare function passed to qsortx */ static int compare(const void *a, const void *b) { - if (*(double*)a < *(double*) b) return -1; - else return 1; + if (*(double*)a < *(double*) b) return -1; + else return 1; } /*----------------------------------------------------------------- - Creat sorted list of z-values in zcorn with actnum==1 -*/ -static int createSortedList(double *list, int n, int m, - const double *z[], const int *a[]) + Creat sorted list of z-values in zcorn with actnum==1x */ +static int createSortedList(double *list, int n, int m, + const double *z[], const int *a[]) { - int i,j; - double *ptr = list; - for (i=0; i apart in of increasing - doubles. -*/ + Remove points less than apart in of increasing + doubles. */ static int uniquify(int n, double *list, double tolerance) { - int i; - int pos; - double val; - - assert (!(tolerance < 0.0)); + int i; + int pos; + double val; - if (n<1) return 0; - pos = 0; - val = list[pos++];/* Keep first value */ - - for (i=1; i - away from the last point (list[n-1]), we remove this - second-to-last point as it cannot be distinguished from "final" - point. - */ - list[pos-1] = list[n-1]; + list[pos-2] + tolerance < list[n-1]. - return pos; + If, however, the second to last point is less than + away from the last point (list[n-1]), we remove this + second-to-last point as it cannot be distinguished from "final" + point. + */ + list[pos-1] = list[n-1]; + + return pos; } /*----------------------------------------------------------------- - Along single pillar: -*/ -static int assignPointNumbers(int begin, - int end, - const double *zlist, - int n, - const double *zcorn, - const int *actnum, - int *plist, - double tolerance) + Along single pillar: */ +static int assignPointNumbers(int begin, + int end, + const double *zlist, + int n, + const double *zcorn, + const int *actnum, + int *plist, + double tolerance) { - /* n - number of cells */ - /* zlist - list of len unique z-values */ - /* start - number of unique z-values processed before. */ + /* n - number of cells */ + /* zlist - list of len unique z-values */ + /* start - number of unique z-values processed before. */ - int i, k; - /* All points should now be within tolerance of a listed point. */ + int i, k; + /* All points should now be within tolerance of a listed point. */ - const double *z = zcorn; - const int *a = actnum; - int *p = plist; + const double *z = zcorn; + const int *a = actnum; + int *p = plist; - k = begin; - *p++ = INT_MIN; /* Padding to ease processing of faults */ - for (i=0; i with k index running faster than i running - faster than j, and Cartesian dimensions , find pointers to the - (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of - field. - */ -static void igetvectors(const int dims[3], int i, int j, - const int *field, const int *v[]) + Given a vector with k index running faster than i running + faster than j, and Cartesian dimensions , find pointers to the + (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of + field. */ +static void igetvectors(const int dims[3], int i, int j, + const int *field, const int *v[]) { - - int im = max(1, i ) - 1; - int ip = min(dims[0], i+1) - 1; - int jm = max(1, j ) - 1; - int jp = min(dims[1], j+1) - 1; - - v[0] = field + dims[2]*(im + dims[0]* jm); - v[1] = field + dims[2]*(im + dims[0]* jp); - v[2] = field + dims[2]*(ip + dims[0]* jm); - v[3] = field + dims[2]*(ip + dims[0]* jp); + + int im = max(1, i ) - 1; + int ip = min(dims[0], i+1) - 1; + int jm = max(1, j ) - 1; + int jp = min(dims[1], j+1) - 1; + + v[0] = field + dims[2]*(im + dims[0]* jm); + v[1] = field + dims[2]*(im + dims[0]* jp); + v[2] = field + dims[2]*(ip + dims[0]* jm); + v[3] = field + dims[2]*(ip + dims[0]* jp); } /*----------------------------------------------------------------- - Given a vector with k index running faster than i running - faster than j, and Cartesian dimensions , find pointers to the - (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of - field. - */ -static void dgetvectors(const int dims[3], int i, int j, - const double *field, const double *v[]) + Given a vector with k index running faster than i running + faster than j, and Cartesian dimensions , find pointers to the + (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of + field. */ +static void dgetvectors(const int dims[3], int i, int j, + const double *field, const double *v[]) { - - int im = max(1, i ) - 1; - int ip = min(dims[0], i+1) - 1; - int jm = max(1, j ) - 1; - int jp = min(dims[1], j+1) - 1; - - v[0] = field + dims[2]*(im + dims[0]* jm); - v[1] = field + dims[2]*(im + dims[0]* jp); - v[2] = field + dims[2]*(ip + dims[0]* jm); - v[3] = field + dims[2]*(ip + dims[0]* jp); + + int im = max(1, i ) - 1; + int ip = min(dims[0], i+1) - 1; + int jm = max(1, j ) - 1; + int jp = min(dims[1], j+1) - 1; + + v[0] = field + dims[2]*(im + dims[0]* jm); + v[1] = field + dims[2]*(im + dims[0]* jp); + v[2] = field + dims[2]*(ip + dims[0]* jm); + v[3] = field + dims[2]*(ip + dims[0]* jp); } /*----------------------------------------------------------------- Given a z coordinate, find x and y coordinates on line defined by coord. Coord points to a vector of 6 doubles [x0,y0,z0,x1,y1,z1]. - */ + */ static void interpolate_pillar(const double *coord, double *pt) { - double a = (pt[2]-coord[2])/(coord[5]-coord[2]); - if (isinf(a) || isnan(a)){ - a = 0; - } + double a = (pt[2]-coord[2])/(coord[5]-coord[2]); + if (isinf(a) || isnan(a)){ + a = 0; + } #if 0 - pt[0] = coord[0] + a*(coord[3]-coord[0]); - pt[1] = coord[1] + a*(coord[4]-coord[1]); + pt[0] = coord[0] + a*(coord[3]-coord[0]); + pt[1] = coord[1] + a*(coord[4]-coord[1]); #else - pt[0] = (1.0 - a)*coord[0] + a*coord[3]; - pt[1] = (1.0 - a)*coord[1] + a*coord[4]; + pt[0] = (1.0 - a)*coord[0] + a*coord[3]; + pt[1] = (1.0 - a)*coord[1] + a*coord[4]; #endif } /*----------------------------------------------------------------- - Assign point numbers p such that "zlist(p)==zcorn". - Assume that coordinate number is arranged in a - sequence such that the natural index is (k,i,j) -*/ + Assign point numbers p such that "zlist(p)==zcorn". Assume that + coordinate number is arranged in a sequence such that the natural + index is (k,i,j) */ int finduniquepoints(const struct grdecl *g, - /* return values: */ - int *plist, /* list of point numbers on each pillar*/ - double tolerance, - struct processed_grid *out) - + /* return values: */ + int *plist, /* list of point numbers on + * each pillar*/ + double tolerance, + struct processed_grid *out) + { - - const int nx = out->dimensions[0]; - const int ny = out->dimensions[1]; - const int nz = out->dimensions[2]; - const int nc = g->dims[0]*g->dims[1]*g->dims[2]; + + const int nx = out->dimensions[0]; + const int ny = out->dimensions[1]; + const int nz = out->dimensions[2]; + const int nc = g->dims[0]*g->dims[1]*g->dims[2]; - /* ztab->data may need extra space temporarily due to simple boundary treatement */ - int npillarpoints = 8*(nx+1)*(ny+1)*nz; - int npillars = (nx+1)*(ny+1); + /* zlist may need extra space temporarily due to simple boundary + * treatement */ + int npillarpoints = 8*(nx+1)*(ny+1)*nz; + int npillars = (nx+1)*(ny+1); - double *zlist = malloc(npillarpoints*sizeof *zlist); - int *zptr = malloc((npillars+1)*sizeof *zptr); + double *zlist = malloc(npillarpoints*sizeof *zlist); + int *zptr = malloc((npillars+1)*sizeof *zptr); - int i,j,k; + int i,j,k; - int d1[3]; - int len = 0; - double *zout = zlist; - int pos = 0; - double *coord = (double*)g->coord; - double *pt; - const double *z[4]; - const int *a[4]; - int *p; - int pix, cix; - int zix; + int d1[3]; + int len = 0; + double *zout = zlist; + int pos = 0; + double *coord = (double*)g->coord; + double *pt; + const double *z[4]; + const int *a[4]; + int *p; + int pix, cix; + int zix; - d1[0] = 2*g->dims[0]; - d1[1] = 2*g->dims[1]; - d1[2] = 2*g->dims[2]; + d1[0] = 2*g->dims[0]; + d1[1] = 2*g->dims[1]; + d1[2] = 2*g->dims[2]; - out->node_coordinates = malloc (3*8*nc*sizeof(*out->node_coordinates)); + out->node_coordinates = malloc (3*8*nc*sizeof(*out->node_coordinates)); - zptr[pos++] = zout - zlist; + zptr[pos++] = zout - zlist; - pt = out->node_coordinates; + pt = out->node_coordinates; - /* Loop over pillars, find unique points on each pillar */ - for (j=0; j < g->dims[1]+1; ++j){ - for (i=0; i < g->dims[0]+1; ++i){ - - /* Get positioned pointers for actnum and zcorn data */ - igetvectors(g->dims, i, j, g->actnum, a); - dgetvectors(d1, 2*i, 2*j, g->zcorn, z); + /* Loop over pillars, find unique points on each pillar */ + for (j=0; j < g->dims[1]+1; ++j){ + for (i=0; i < g->dims[0]+1; ++i){ - len = createSortedList( zout, d1[2], 4, z, a); - len = uniquify (len, zout, tolerance); + /* Get positioned pointers for actnum and zcorn data */ + igetvectors(g->dims, i, j, g->actnum, a); + dgetvectors(d1, 2*i, 2*j, g->zcorn, z); - /* Assign unique points */ - for (k=0; knumber_of_nodes_on_pillars = zptr[pos-1]; - out->number_of_nodes = zptr[pos-1]; + out->number_of_nodes_on_pillars = zptr[pos-1]; + out->number_of_nodes = zptr[pos-1]; - /* Loop over all vertical sets of zcorn values, assign point numbers */ - p = plist; - for (j=0; j < 2*g->dims[1]; ++j){ - for (i=0; i < 2*g->dims[0]; ++i){ - - /* pillar index */ - pix = (i+1)/2 + (g->dims[0]+1)*((j+1)/2); - - /* cell column position */ - cix = g->dims[2]*((i/2) + (j/2)*g->dims[0]); + /* Loop over all vertical sets of zcorn values, assign point + * numbers */ + p = plist; + for (j=0; j < 2*g->dims[1]; ++j){ + for (i=0; i < 2*g->dims[0]; ++i){ - /* zcorn column position */ - zix = 2*g->dims[2]*(i+2*g->dims[0]*j); - - if (!assignPointNumbers(zptr[pix], zptr[pix+1], zlist, - 2*g->dims[2], - g->zcorn + zix, g->actnum + cix, - p, tolerance)){ - fprintf(stderr, "Something went wrong in assignPointNumbers"); - return 0; - } + /* pillar index */ + pix = (i+1)/2 + (g->dims[0]+1)*((j+1)/2); - p += 2 + 2*g->dims[2]; + /* cell column position */ + cix = g->dims[2]*((i/2) + (j/2)*g->dims[0]); + + /* zcorn column position */ + zix = 2*g->dims[2]*(i+2*g->dims[0]*j); + + if (!assignPointNumbers(zptr[pix], zptr[pix+1], zlist, + 2*g->dims[2], + g->zcorn + zix, g->actnum + cix, + p, tolerance)){ + fprintf(stderr, "Something went wrong in assignPointNumbers"); + return 0; + } + + p += 2 + 2*g->dims[2]; + } } - } - free(zptr); - free(zlist); + free(zptr); + free(zlist); - return 1; + return 1; } /* Local Variables: */ /* c-basic-offset:4 */ /* End: */ - diff --git a/uniquepoints.h b/uniquepoints.h index ac3574bf..b9d698ac 100644 --- a/uniquepoints.h +++ b/uniquepoints.h @@ -13,34 +13,34 @@ //==========================================================================*/ /* -Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. -Copyright 2009, 2010 Statoil ASA. + Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. + Copyright 2009, 2010 Statoil ASA. -This file is part of The Open Reservoir Simulator Project (OpenRS). + This file is part of The Open Reservoir Simulator Project (OpenRS). -OpenRS 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. + OpenRS 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. -OpenRS 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. + OpenRS 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 OpenRS. If not, see . + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ -#ifndef OPENRS_UNIQUEPOINTS_HEADER -#define OPENRS_UNIQUEPOINTS_HEADER +#ifndef OPM_UNIQUEPOINTS_HEADER +#define OPM_UNIQUEPOINTS_HEADER int finduniquepoints(const struct grdecl *g, /* input */ - int *p, /* for each z0 in zcorn, z0 = z[p0] */ - double t, /* tolerance*/ - struct processed_grid *out); + int *p, /* for each z0 in zcorn, z0 = z[p0] */ + double t, /* tolerance*/ + struct processed_grid *out); -#endif /* OPENRS_UNIQUEPOINTS_HEADER */ +#endif /* OPM_UNIQUEPOINTS_HEADER */ /* Local Variables: */ /* c-basic-offset:4 */