diff --git a/Makefile.am b/Makefile.am index 583e67ff..22bc09ec 100644 --- a/Makefile.am +++ b/Makefile.am @@ -42,6 +42,7 @@ $(ERT_LDFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) # # Please try to keep the list sorted. +<<<<<<< HEAD libopmcore_la_SOURCES = \ opm/core/GridManager.cpp \ opm/core/eclipse/EclipseGridInspector.cpp \ @@ -265,6 +266,227 @@ opm/core/wells/InjectionSpecification.hpp \ opm/core/wells/ProductionSpecification.hpp \ opm/core/wells/WellCollection.hpp \ opm/core/wells/WellsGroup.hpp \ +======= +libopmcore_la_SOURCES = \ +opm/core/GridManager.cpp \ +opm/core/eclipse/EclipseGridInspector.cpp \ +opm/core/eclipse/EclipseGridParser.cpp \ +opm/core/fluid/BlackoilPropertiesBasic.cpp \ +opm/core/fluid/BlackoilPropertiesFromDeck.cpp \ +opm/core/fluid/IncompPropertiesBasic.cpp \ +opm/core/fluid/IncompPropertiesFromDeck.cpp \ +opm/core/fluid/PvtPropertiesBasic.cpp \ +opm/core/fluid/PvtPropertiesIncompFromDeck.cpp \ +opm/core/fluid/RockBasic.cpp \ +opm/core/fluid/RockCompressibility.cpp \ +opm/core/fluid/RockFromDeck.cpp \ +opm/core/fluid/SaturationPropsBasic.cpp \ +opm/core/fluid/SaturationPropsFromDeck.cpp \ +opm/core/fluid/blackoil/BlackoilPvtProperties.cpp \ +opm/core/fluid/blackoil/SinglePvtDead.cpp \ +opm/core/fluid/blackoil/SinglePvtInterface.cpp \ +opm/core/fluid/blackoil/SinglePvtLiveGas.cpp \ +opm/core/fluid/blackoil/SinglePvtLiveOil.cpp \ +opm/core/grid.c \ +opm/core/grid/cart_grid.c \ +opm/core/grid/cornerpoint_grid.c \ +opm/core/grid/cpgpreprocess/facetopology.c \ +opm/core/grid/cpgpreprocess/geometry.c \ +opm/core/grid/cpgpreprocess/preprocess.c \ +opm/core/grid/cpgpreprocess/uniquepoints.c \ +opm/core/linalg/LinearSolverFactory.cpp \ +opm/core/linalg/LinearSolverInterface.cpp \ +opm/core/linalg/sparse_sys.c \ +opm/core/newwells.c \ +opm/core/pressure/CompressibleTpfa.cpp \ +opm/core/pressure/FlowBCManager.cpp \ +opm/core/pressure/IncompTpfa.cpp \ +opm/core/pressure/cfsh.c \ +opm/core/pressure/flow_bc.c \ +opm/core/pressure/fsh.c \ +opm/core/pressure/fsh_common_impl.c \ +opm/core/pressure/ifsh.c \ +opm/core/pressure/mimetic/hybsys.c \ +opm/core/pressure/mimetic/hybsys_global.c \ +opm/core/pressure/mimetic/mimetic.c \ +opm/core/pressure/msmfem/coarse_conn.c \ +opm/core/pressure/msmfem/coarse_sys.c \ +opm/core/pressure/msmfem/dfs.c \ +opm/core/pressure/msmfem/hash_set.c \ +opm/core/pressure/msmfem/ifsh_ms.c \ +opm/core/pressure/msmfem/partition.c \ +opm/core/pressure/tpfa/cfs_tpfa.c \ +opm/core/pressure/tpfa/cfs_tpfa_residual.c \ +opm/core/pressure/tpfa/compr_bc.c \ +opm/core/pressure/tpfa/compr_quant.c \ +opm/core/pressure/tpfa/compr_quant_general.c \ +opm/core/pressure/tpfa/compr_source.c \ +opm/core/pressure/tpfa/ifs_tpfa.c \ +opm/core/pressure/tpfa/trans_tpfa.c \ +opm/core/pressure/well.c \ +opm/core/simulator/SimulatorReport.cpp \ +opm/core/simulator/SimulatorTimer.cpp \ +opm/core/simulator/SimulatorTwophase.cpp \ +opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp \ +opm/core/transport/reorder/TransportModelInterface.cpp \ +opm/core/transport/reorder/TransportModelTwophase.cpp \ +opm/core/transport/reorder/nlsolvers.c \ +opm/core/transport/reorder/reordersequence.cpp \ +opm/core/transport/reorder/tarjan.c \ +opm/core/transport/spu_explicit.c \ +opm/core/transport/spu_implicit.c \ +opm/core/transport/transport_source.c \ +opm/core/utility/MonotCubicInterpolator.cpp \ +opm/core/utility/StopWatch.cpp \ +opm/core/utility/miscUtilities.cpp \ +opm/core/utility/miscUtilitiesBlackoil.cpp \ +opm/core/utility/parameters/Parameter.cpp \ +opm/core/utility/parameters/ParameterGroup.cpp \ +opm/core/utility/parameters/ParameterTools.cpp \ +opm/core/utility/parameters/ParameterXML.cpp \ +opm/core/utility/parameters/tinyxml/tinystr.cpp \ +opm/core/utility/parameters/tinyxml/tinyxml.cpp \ +opm/core/utility/parameters/tinyxml/tinyxmlerror.cpp \ +opm/core/utility/parameters/tinyxml/tinyxmlparser.cpp \ +opm/core/utility/writeVtkData.cpp \ +opm/core/vag_format/vag.cpp \ +opm/core/wells/InjectionSpecification.cpp \ +opm/core/wells/ProductionSpecification.cpp \ +opm/core/wells/WellCollection.cpp \ +opm/core/wells/WellsGroup.cpp \ +opm/core/wells/WellsManager.cpp + +nobase_include_HEADERS = \ +opm/core/GridAdapter.hpp \ +opm/core/GridManager.hpp \ +opm/core/eclipse/CornerpointChopper.hpp \ +opm/core/eclipse/EclipseGridInspector.hpp \ +opm/core/eclipse/EclipseGridParser.hpp \ +opm/core/eclipse/EclipseGridParserHelpers.hpp \ +opm/core/eclipse/EclipseUnits.hpp \ +opm/core/eclipse/SpecialEclipseFields.hpp \ +opm/core/fluid/BlackoilPropertiesBasic.hpp \ +opm/core/fluid/BlackoilPropertiesFromDeck.hpp \ +opm/core/fluid/BlackoilPropertiesInterface.hpp \ +opm/core/fluid/IncompPropertiesBasic.hpp \ +opm/core/fluid/IncompPropertiesFromDeck.hpp \ +opm/core/fluid/IncompPropertiesInterface.hpp \ +opm/core/fluid/PvtPropertiesBasic.hpp \ +opm/core/fluid/PvtPropertiesIncompFromDeck.hpp \ +opm/core/fluid/RockBasic.hpp \ +opm/core/fluid/RockCompressibility.hpp \ +opm/core/fluid/RockFromDeck.hpp \ +opm/core/fluid/SaturationPropsBasic.hpp \ +opm/core/fluid/SaturationPropsFromDeck.hpp \ +opm/core/fluid/SimpleFluid2p.hpp \ +opm/core/fluid/blackoil/BlackoilPhases.hpp \ +opm/core/fluid/blackoil/BlackoilPvtProperties.hpp \ +opm/core/fluid/blackoil/SinglePvtConstCompr.hpp \ +opm/core/fluid/blackoil/SinglePvtDead.hpp \ +opm/core/fluid/blackoil/SinglePvtInterface.hpp \ +opm/core/fluid/blackoil/SinglePvtLiveGas.hpp \ +opm/core/fluid/blackoil/SinglePvtLiveOil.hpp \ +opm/core/fluid/blackoil/phaseUsageFromDeck.hpp \ +opm/core/grid.h \ +opm/core/grid/cart_grid.h \ +opm/core/grid/cornerpoint_grid.h \ +opm/core/grid/cpgpreprocess/facetopology.h \ +opm/core/grid/cpgpreprocess/geometry.h \ +opm/core/grid/cpgpreprocess/grdecl.h \ +opm/core/grid/cpgpreprocess/preprocess.h \ +opm/core/grid/cpgpreprocess/uniquepoints.h \ +opm/core/linalg/LinearSolverFactory.hpp \ +opm/core/linalg/LinearSolverInterface.hpp \ +opm/core/linalg/blas_lapack.h \ +opm/core/linalg/sparse_sys.h \ +opm/core/newwells.h \ +opm/core/pressure/CompressibleTpfa.hpp \ +opm/core/pressure/FlowBCManager.hpp \ +opm/core/pressure/HybridPressureSolver.hpp \ +opm/core/pressure/IncompTpfa.hpp \ +opm/core/pressure/TPFACompressiblePressureSolver.hpp \ +opm/core/pressure/TPFAPressureSolver.hpp \ +opm/core/pressure/flow_bc.h \ +opm/core/pressure/fsh.h \ +opm/core/pressure/fsh_common_impl.h \ +opm/core/pressure/mimetic/hybsys.h \ +opm/core/pressure/mimetic/hybsys_global.h \ +opm/core/pressure/mimetic/mimetic.h \ +opm/core/pressure/msmfem/coarse_conn.h \ +opm/core/pressure/msmfem/coarse_sys.h \ +opm/core/pressure/msmfem/dfs.h \ +opm/core/pressure/msmfem/hash_set.h \ +opm/core/pressure/msmfem/ifsh_ms.h \ +opm/core/pressure/msmfem/partition.h \ +opm/core/pressure/tpfa/cfs_tpfa.h \ +opm/core/pressure/tpfa/cfs_tpfa_residual.h \ +opm/core/pressure/tpfa/compr_bc.h \ +opm/core/pressure/tpfa/compr_quant.h \ +opm/core/pressure/tpfa/compr_quant_general.h \ +opm/core/pressure/tpfa/compr_source.h \ +opm/core/pressure/tpfa/ifs_tpfa.h \ +opm/core/pressure/tpfa/trans_tpfa.h \ +opm/core/simulator/BlackoilState.hpp \ +opm/core/simulator/SimulatorReport.hpp \ +opm/core/simulator/SimulatorTimer.hpp \ +opm/core/simulator/SimulatorTwophase.hpp \ +opm/core/simulator/TwophaseState.hpp \ +opm/core/simulator/WellState.hpp \ +opm/core/transport/CSRMatrixBlockAssembler.hpp \ +opm/core/transport/CSRMatrixUmfpackSolver.hpp \ +opm/core/transport/GravityColumnSolver.hpp \ +opm/core/transport/GravityColumnSolver_impl.hpp \ +opm/core/transport/ImplicitAssembly.hpp \ +opm/core/transport/ImplicitTransport.hpp \ +opm/core/transport/JacobianSystem.hpp \ +opm/core/transport/NormSupport.hpp \ +opm/core/transport/SimpleFluid2pWrapper.hpp \ +opm/core/transport/SinglePointUpwindTwoPhase.hpp \ +opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp \ +opm/core/transport/reorder/TransportModelInterface.hpp \ +opm/core/transport/reorder/TransportModelTwophase.hpp \ +opm/core/transport/reorder/nlsolvers.h \ +opm/core/transport/reorder/reordersequence.h \ +opm/core/transport/reorder/tarjan.h \ +opm/core/transport/spu_explicit.h \ +opm/core/transport/spu_implicit.h \ +opm/core/transport/transport_source.h \ +opm/core/utility/Average.hpp \ +opm/core/utility/ColumnExtract.hpp \ +opm/core/utility/ErrorMacros.hpp \ +opm/core/utility/Factory.hpp \ +opm/core/utility/MonotCubicInterpolator.hpp \ +opm/core/utility/RootFinders.hpp \ +opm/core/utility/SparseTable.hpp \ +opm/core/utility/SparseVector.hpp \ +opm/core/utility/StopWatch.hpp \ +opm/core/utility/UniformTableLinear.hpp \ +opm/core/utility/Units.hpp \ +opm/core/utility/buildUniformMonotoneTable.hpp \ +opm/core/utility/initState.hpp \ +opm/core/utility/initState_impl.hpp \ +opm/core/utility/linInt.hpp \ +opm/core/utility/linearInterpolation.hpp \ +opm/core/utility/miscUtilities.hpp \ +opm/core/utility/miscUtilitiesBlackoil.hpp \ +opm/core/utility/parameters/Parameter.hpp \ +opm/core/utility/parameters/ParameterGroup.hpp \ +opm/core/utility/parameters/ParameterGroup_impl.hpp \ +opm/core/utility/parameters/ParameterMapItem.hpp \ +opm/core/utility/parameters/ParameterRequirement.hpp \ +opm/core/utility/parameters/ParameterStrings.hpp \ +opm/core/utility/parameters/ParameterTools.hpp \ +opm/core/utility/parameters/ParameterXML.hpp \ +opm/core/utility/parameters/tinyxml/tinystr.h \ +opm/core/utility/parameters/tinyxml/tinyxml.h \ +opm/core/utility/writeVtkData.hpp \ +opm/core/vag_format/vag.hpp \ +opm/core/well.h \ +opm/core/wells/InjectionSpecification.hpp \ +opm/core/wells/ProductionSpecification.hpp \ +opm/core/wells/WellCollection.hpp \ +opm/core/wells/WellsGroup.hpp \ +>>>>>>> master opm/core/wells/WellsManager.hpp # ---------------------------------------------------------------------- diff --git a/opm/core/grid/cpgpreprocess/facetopology.c b/opm/core/grid/cpgpreprocess/facetopology.c index 6898e2b4..15c71c9a 100644 --- a/opm/core/grid/cpgpreprocess/facetopology.c +++ b/opm/core/grid/cpgpreprocess/facetopology.c @@ -13,39 +13,38 @@ //===========================================================================*/ /* -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 -#include -#include #include -#include #include +#include +#include +#include +#include #include "preprocess.h" -#include "sparsetable.h" #include "facetopology.h" /* No checking of input arguments in this code! */ -#define min(i,j) ((i)<(j) ? (i) : (j)) -#define max(i,j) ((i)>(j) ? (i) : (j)) +#define MIN(i,j) ((i)<(j) ? (i) : (j)) +#define MAX(i,j) ((i)>(j) ? (i) : (j)) #define DEBUG 1 @@ -59,296 +58,309 @@ along with OpenRS. If not, see . /*------------------------------------------------------*/ -/* Determine face geometry first, then compute intersections. */ +/* 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) +static int * +computeFaceTopology(const int *a1, const int *a2, + const int *b1, const int *b2, + int intersect[4], int *faces) { - int mask[8] = {-1}; - int k, *f; + int mask[8]; + 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]; } + for (k = 0; k < 8; k++) { mask[k] = -1; } + + /* 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. +/* 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. - b) We assume no intersections occur on the first and last lines. - This is convenient in the identification of (unique) intersections. + 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. */ -#define lineintersection(a1,a2,b1,b2)(((a1>b1)&&(a2b2))) -static int faceintersection(int *a1, int *a2, int *b1, int *b2) +#define LINE_INTERSECTION(a1, a2, b1, b2) \ + (((a1 > b1) && (a2 < b2)) || \ + ((a1 < b1) && (a2 > b2))) + +static int +faceintersection(const int *a1, const int *a2, + const int *b1, const 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]) || + LINE_INTERSECTION(a1[0], a2[0], b1[0], b2[0]) || + LINE_INTERSECTION(a1[1], a2[1], b1[1], b2[1]); } -#define meaningful_face !((a1[i] ==INT_MIN) && (b1[j] ==INT_MIN)) && \ - !((a1[i+1]==INT_MAX) && (b1[j+1]==INT_MAX)) +#define MEANINGFUL_FACE(i, j) \ + (! ((a1[i] == INT_MIN) && (b1[j] == INT_MIN)) && \ + ! ((a1[i+1] == INT_MAX) && (b1[j+1] == INT_MAX))) + /* 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; - int *tmp; + /* 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}; - /* 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 (LINE_INTERSECTION(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; - while(jface_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; } - /* --------------------------------------------------------- */ - /* face a(i,i+1) and face b(j,j+1) have nonzero intersection */ - /* --------------------------------------------------------- */ - if (faceintersection(a1+i, a2+i, b1+j, b2+j)){ - - - /* Completely matching faces */ - if (a1[i]==b1[j] && a1[i+1]==b1[j+1] && - a2[i]==b2[j] && a2[i+1]==b2[j+1]){ - - /* Add face to list of faces if not any first or last points are involved. */ - if (meaningful_face){ - - int cell_a = i%2 ? (i-1)/2 : -1; - int cell_b = j%2 ? (j-1)/2 : -1; - - if (cell_a != -1 || cell_b != -1){ - *c++ = cell_a; - *c++ = cell_b; - - /* face */ - *f++ = a1[i]; - *f++ = a2[i]; - /* avoid duplicating nodes in pinched faces */ - if (a2[i+1] != a2[i]) *f++ = a2[i+1]; - if (a1[i+1] != a1[i]) *f++ = a1[i+1]; - - out->face_ptr[++out->number_of_faces] = f - out->face_nodes; - - } - else{ - ; - /* - fprintf(stderr, - "Warning. For some reason we get face connecting void spaces\n"); - */ - } - - } - } - - /* 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){ - - /* - Even indices refer to space between cells, - odd indices refer to cells - */ - int cell_a = i%2 ? (i-1)/2 : -1; - int cell_b = j%2 ? (j-1)/2 : -1; + j = j+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); + /* 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; - out->face_ptr[++out->number_of_faces] = f - out->face_nodes; - } - else{ - ; - /* - fprintf(stderr, - "Warning. For some reason we get face connecting void spaces\n"); - */ - } - } - } - } + /* 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/opm/core/grid/cpgpreprocess/geometry.c b/opm/core/grid/cpgpreprocess/geometry.c index 43abeed9..a49fc874 100644 --- a/opm/core/grid/cpgpreprocess/geometry.c +++ b/opm/core/grid/cpgpreprocess/geometry.c @@ -131,6 +131,12 @@ compute_cell_geometry(int ndims, double *coords, double ccell[3]; double cface[3] = {0}; const double twothirds = 0.666666666666666666666666666667; + + int ndigits; + + ndigits = ((int) (log(ncells) / log(10.0))) + 1; + + for (c=0; c0)){ - fprintf(stderr, "Internal error in compute_cell_geometry: negative volume\n"); - } /* face centroid of triangle */ for (i=0; i 0.0)) { + fprintf(stderr, + "Internal error in mex_compute_geometry(%*d): " + "negative volume\n", ndigits, c); + } + for (i=0; i. + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ -#include -#include -#include #include +#include +#include +#include +#include + #include - -#include - - +#include "grdecl.h" +#include "mxgrdecl.h" /* Get COORD, ZCORN, ACTNUM and DIMS from mxArray. */ /*-------------------------------------------------------*/ void mx_init_grdecl(struct grdecl *g, const mxArray *s) { - int i,n; - mxArray *field; - int numel; + int i,n; + size_t numel; + mxArray *cartdims=NULL, *actnum=NULL, *coord=NULL, *zcorn=NULL; - 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"); - } - - double *tmp = mxGetPr(cartdims); - n = 1; - for (i=0; i<3; ++i){ - g->dims[i] = tmp[i]; - n *= tmp[i]; - } + if (!mxIsStruct(s) + || !(cartdims = mxGetField(s, 0, "cartDims")) + || !(coord = mxGetField(s, 0, "COORD")) + || !(zcorn = mxGetField(s, 0, "ZCORN")) + ) + { + char str[]="Input must be a single MATLAB struct with fields\n" + "'cartDims', 'COORD' and 'ZCORN'. ACTNUM may be included.\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"); + } + + if (mxIsDouble(cartdims)) { + double *tmp = mxGetPr(cartdims); + for (i = 0; i < 3; ++i) { + g->dims[i] = (int) tmp[i]; + } + } + else if (mxIsInt32(cartdims)) { + int *tmp = mxGetData(cartdims); + memcpy(g->dims, tmp, 3 * sizeof *g->dims); + } + + n = g->dims[0]; + for (i = 1; i < 3; i++) { n *= g->dims[ 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); - + if ((actnum = mxGetField(s, 0, "ACTNUM")) != NULL) { + numel = mxGetNumberOfElements(actnum); + if ((! mxIsInt32(actnum)) || (numel != (size_t)(n))) { + mexErrMsgTxt("ACTNUM field must be nx*ny*nz numbers int32"); + } + g->actnum = mxGetData(actnum); + } + else { + g->actnum = NULL; + } - 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); + + numel = mxGetNumberOfElements(coord); + if ((! mxIsDouble(coord)) || + numel != (size_t)(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 ((! mxIsDouble(zcorn)) || + numel != (size_t)(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/opm/core/grid/cpgpreprocess/mxgrdecl.h b/opm/core/grid/cpgpreprocess/mxgrdecl.h index 54fa2267..4621c202 100644 --- a/opm/core/grid/cpgpreprocess/mxgrdecl.h +++ b/opm/core/grid/cpgpreprocess/mxgrdecl.h @@ -1,4 +1,4 @@ -//=========================================================================== +/*========================================================================= // // File: mxgrdecl.h // @@ -10,33 +10,35 @@ // // $Revision$ // -//=========================================================================== +//=======================================================================*/ /* -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/opm/core/grid/cpgpreprocess/preprocess.c b/opm/core/grid/cpgpreprocess/preprocess.c index f5048a8f..bcd39c15 100644 --- a/opm/core/grid/cpgpreprocess/preprocess.c +++ b/opm/core/grid/cpgpreprocess/preprocess.c @@ -13,59 +13,89 @@ //==========================================================================*/ /* -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 . */ -#include -#include -#include #include -#include #include +#include +#include +#include +#include + #include "preprocess.h" #include "uniquepoints.h" #include "facetopology.h" -#define min(i,j) ((i)<(j) ? (i) : (j)) -#define max(i,j) ((i)>(j) ? (i) : (j)) +#define MIN(i,j) ((i)<(j) ? (i) : (j)) +#define MAX(i,j) ((i)>(j) ? (i) : (j)) +static void +compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len); + +static int +checkmemory(int nz, struct processed_grid *out, int **intersections); + +static void +process_vertical_faces(int direction, + int **intersections, + int *plist, int *work, + struct processed_grid *out); + +static void +process_horizontal_faces(int **intersections, + int *plist, + struct processed_grid *out); + +static int +linearindex(const int dims[3], int i, int j, int k) +{ + assert (0 <= i); + assert (0 <= j); + assert (0 <= k); + + assert (i < dims[0]); + assert (j < dims[1]); + assert (k < dims[2]); + + return i + dims[0]*(j + dims[1]*k); +} /*----------------------------------------------------------------- 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(int dims[3], int i, int j, int *field, int *v[]) + 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); } @@ -76,263 +106,295 @@ 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. - */ -static void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len) +*/ +static void +compute_cell_index(const int dims[3], int i, int j, + int *neighbors, int len) { + int k; - int k; - if (i<0 || i>=dims[0] || j<0 || j >= dims[1]){ - for(k=0; k= dims[0])) || /* 'i' outside [0, dims[0]) */ + ((j < 0) || (j >= dims[1]))) { /* 'j' outside [0, dims[1]) */ + + for (k = 0; k < len; k += 2) { + neighbors[k] = -1; /* Neighbour is outside domain */ + } } - }else{ - for(k=0; km; - int n = out->n; + /* Ensure there is enough space to manage the (pathological) case + * of every single cell on one side of a fault connecting to all + * cells on the other side of the fault (i.e., an all-to-all cell + * connectivity pairing). */ + r = (2*nz + 2) * (2*nz + 2); + m = out->m; + 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 / 2, 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 / 2, 12 * r); } - out->m = m; - out->n = n; - } + ok = m == out->m; + if (! ok) { + void *p1, *p2, *p3, *p4; - return 1; + p1 = realloc(*intersections , 4*m * sizeof **intersections); + p2 = realloc(out->face_neighbors, 2*m * sizeof *out->face_neighbors); + p3 = realloc(out->face_ptr , (m+1) * sizeof *out->face_ptr); + p4 = realloc(out->face_tag , 1*m * sizeof *out->face_tag); + + if (p1 != NULL) { *intersections = p1; } + if (p2 != NULL) { out->face_neighbors = p2; } + if (p3 != NULL) { out->face_ptr = p3; } + if (p4 != NULL) { out->face_tag = p4; } + + ok = (p1 != NULL) && (p2 != NULL) && (p3 != NULL) && (p4 != NULL); + + if (ok) { out->m = m; } + } + + if (ok && (n != out->n)) { + void *p1; + + p1 = realloc(out->face_nodes, n * sizeof *out->face_nodes); + + ok = p1 != NULL; + + if (ok) { + out->face_nodes = p1; + out->n = n; + } + } + + return ok; } /*----------------------------------------------------------------- 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. - */ -static void process_vertical_faces(int direction, - int **intersections, - int *plist, int *work, - struct processed_grid *out) +*/ +static void +process_vertical_faces(int direction, + int **intersections, + int *plist, int *work, + struct processed_grid *out) { - int i,j; - int d[3]; - int *cornerpts[4]; + 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; - enum face_tag tag[] = { LEFT, BACK }; + assert ((direction == 0) || (direction == 1)); - int nx = out->dimensions[0]; - int ny = out->dimensions[1]; - int nz = out->dimensions[2]; + d[0] = 2 * (nx + 0); + d[1] = 2 * (ny + 0); + d[2] = 2 * (nz + 1); - int startface, num_intersections, f, len, *ptr; + for (j = 0; j < ny + direction; ++j) { + for (i = 0; i < nx + (1 - direction); ++i) { - 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 */ - int *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; + /* Establish new connections (faces) along pillar pair. */ + findconnections(2*nz + 2, cornerpts, + *intersections + 4*num_intersections, + work, out); - findconnections(2*nz+2, cornerpts, - *intersections+4*num_intersections, - work, out); + /* Start of ->face_neighbors[] for this set of connections. */ + ptr = out->face_neighbors + 2*startface; - ptr = out->face_neighbors + 2*startface; - len = 2*out->number_of_faces - 2*startface; + /* Total number of cells (both sides) connected by this + * set of connections (faces). */ + 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); + /* Derive inter-cell connectivity (i.e. ->face_neighbors) + * of global (uncompressed) cells for this set of + * connections (faces). */ + 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); } /*----------------------------------------------------------------- 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) - */ -static void process_horizontal_faces(int **intersections, - int *plist, - struct processed_grid *out) +*/ +static void +process_horizontal_faces(int **intersections, + 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 *cell = out->local_cell_index; + int cellno = 0; + int *f, *n, *c[4]; + int prevcell, thiscell; + int idx; - int *f, *n, *c[4], prevcell, thiscell; - - /* 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; - - - /* Vectors of point numbers */ - igetvectors(d, 2*i+1, 2*j+1, plist, c); - - prevcell = -1; - - - for (k = 1; kdimensions, 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]; - - 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; - - 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]; - - out->face_tag[ out->number_of_faces] = TOP; - out->face_ptr[++out->number_of_faces] = f - out->face_nodes; - - *n++ = prevcell; - *n++ = prevcell = -1; + if (! checkmemory(nz, out, intersections)) { + fprintf(stderr, + "Could not allocate enough space in " + "process_horizontal_faces()\n"); + exit(1); + } + + + 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); + + prevcell = -1; + + + for (k = 1; kdimensions, 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]; + + 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; + + 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]; + + out->face_tag[ out->number_of_faces] = TOP; + out->face_ptr[++out->number_of_faces] = f - out->face_nodes; + + *n++ = prevcell; + *n++ = prevcell = -1; + } + } + } } - } } - } } - } - out->number_of_cells = cellno; + out->number_of_cells = cellno; } @@ -348,251 +410,538 @@ static void process_horizontal_faces(int **intersections, */ static void approximate_intersection_pt(int *L, double *c, double *pt) { + double a; + double z0, z1, z2, z3; + double b1, b2; + double x1, y1; + double x2, y2; + double z; - 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]; + /* no intersection on pillars expected here! */ + assert (L[0] != L[2]); + assert (L[1] != L[3]); - double a = (z2-z0)/(z1-z0 - (z3-z2)); - if (isinf(a) || isnan(a)){ - a = 0; - } + z0 = c[3*L[0] + 2]; + z1 = c[3*L[1] + 2]; + z2 = c[3*L[2] + 2]; + z3 = c[3*L[3] + 2]; + + /* find parameter a where lines L0L1 and L2L3 have same + * z-coordinate */ + if (fabs((z1 - z0) - (z3 - z2)) > 0.0) { + + a = (z2 - z0) / ((z1 - z0) - (z3 - z2)); + + } else { + + a = 0; + + } + + /* the corresponding z-coordinate is */ + z = z0*(1.0 - a) + z1*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; + /* find point (x1, y1, z) on pillar 1 */ + b1 = (z2 - z) / (z2 - z0); + b2 = (z - z0) / (z2 - z0); + x1 = c[3*L[0] + 0]*b1 + c[3*L[2] + 0]*b2; + y1 = c[3*L[0] + 1]*b1 + c[3*L[2] + 1]*b2; + /* find point (x2, y2, z) on pillar 2 */ + b1 = (z - z3) / (z1 - z3); + b2 = (z1 - z) / (z1 - z3); + x2 = c[3*L[1] + 0]*b1 + c[3*L[3] + 0]*b2; + y2 = c[3*L[1] + 1]*b1 + c[3*L[3] + 1]*b2; + + /* horizontal lines are by definition ON the bilinear surface + spanned by L0, L1, L2 and L3. find point (x, y, z) on + horizontal line between point (x1, y1, z) and (x2, y2, z).*/ + pt[0] = x1*(1.0 - a) + x2*a; + pt[1] = y1*(1.0 - a) + y2*a; + pt[2] = z; } /*----------------------------------------------------------------- - 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) { - int n = out->number_of_nodes; - int np = out->number_of_nodes_on_pillars; - - int k, *itsct; - double *pt; - - /* 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; - itsct = intersections; + /* 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; - } + } +} + + +/* ------------------------------------------------------------------ */ +static int* +copy_and_permute_actnum(int nx, int ny, int nz, const int *in, int *out) +/* ------------------------------------------------------------------ */ +{ + int i,j,k; + int *ptr = 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,:)] + * + * in MATLAB pseudo-code. + */ + if (in != NULL) { + for (j = 0; j < ny; ++j) { + for (i = 0; i < nx; ++i) { + for (k = 0; k < nz; ++k) { + *ptr++ = in[i + nx*(j + ny*k)]; + } + } + } + } + else { + /* No explicit ACTNUM. Assume all cells active. */ + for (i = 0; i < nx * ny * nz; i++) { + out[ i ] = 1; + } + } + + return out; +} + +/* ------------------------------------------------------------------ */ +static double* +copy_and_permute_zcorn(int nx, int ny, int nz, const double *in, + double sign, double *out) +/* ------------------------------------------------------------------ */ +{ + int i,j,k; + double *ptr = out; + /* Permute zcorn such that values of each vertical stack of cells + * are adjacent in memory, i.e., + + out = [in(0,0,:), in(1,0,:),..., in(2*nx-1, 2*ny-1,:)] + + in Matlab pseudo-code. + */ + for (j=0; j<2*ny; ++j){ + for (i=0; i<2*nx; ++i){ + for (k=0; k<2*nz; ++k){ + *ptr++ = sign * in[i+2*nx*(j+2*ny*k)]; + } + } + } + return out; +} + +/* ------------------------------------------------------------------ */ +static int +get_zcorn_sign(int nx, int ny, int nz, const int *actnum, + const double *zcorn, int *error) +/* ------------------------------------------------------------------ */ +{ + /* Ensure that zcorn (i.e., depth) is strictly nondecreasing in + the k-direction. This is required by the processign algorithm. + + 1) if z(i,j,k) <= z(i,j,k+1) for all (i,j,k), return 1.0 + + 2) if -z(i,j,k) <=-z(i,j,k+1) for all (i,j,k), return -1.0 + + 3) if (1) and (2) fails, return -1.0, and set *error = 1. + + */ + int sign; + int i, j, k; + int c1, c2; + double z1, z2; + + for (sign = 1; sign>-2; sign = sign - 2) + { + *error = 0; + + for (j=0; j<2*ny; ++j){ + for (i=0; i<2*nx; ++i){ + for (k=0; k<2*nz-1; ++k){ + z1 = sign*zcorn[i+2*nx*(j+2*ny*(k))]; + z2 = sign*zcorn[i+2*nx*(j+2*ny*(k+1))]; + + c1 = i/2 + nx*(j/2 + ny*k/2); + c2 = i/2 + nx*(j/2 + ny*(k+1)/2); + + if (((actnum == NULL) || + (actnum[c1] && actnum[c2])) + && (z2 < z1)) { + + fprintf(stderr, "\nZCORN should be strictly " + "nondecreasing along pillars!\n"); + *error = 1; + goto end; + } + } + } + } + + end: + if (!*error){ + break; + } + } + + if (*error){ + fprintf(stderr, "Attempt to reverse sign in ZCORN failed.\n" + "Grid definition may be broken\n"); + } + + return sign; +} + + +static void +ind2sub(const size_t nx, + const size_t ny, + const size_t nz, + size_t c , + size_t *i, size_t *j, size_t *k) +{ + assert (c < (nx * ny * nz)); + +#if !defined(NDEBUG) + (void) nz; +#endif + + *i = c % nx; c /= nx; + *j = c % ny; + *k = c / ny; +} + + +/* ---------------------------------------------------------------------- */ +static double +vert_size(const struct grdecl *in, + const size_t c , + const size_t off[8]) +/* ---------------------------------------------------------------------- */ +{ + size_t i, j, k, start; + double dz; + const double *zcorn; + + ind2sub(in->dims[0], in->dims[1], in->dims[2], c, &i, &j, &k); + + zcorn = in->zcorn; + start = (k * off[4]) + (j * off[2]) + (i * 2); + + for (k = 0, dz = 0.0; (! (fabs(dz) > 0)) && (k < 4); k++) { + dz = zcorn[start + off[k + 4]] - zcorn[start + off[k]]; + } + + return dz; +} + + +/* ---------------------------------------------------------------------- */ +static int +is_lefthanded(const struct grdecl *in) +/* ---------------------------------------------------------------------- */ +{ + int active, searching; + size_t nx, ny, nz, c; + size_t origin, imax, jmax; + size_t off[8]; + double dx[2], dy[2], dz, triple; + const double *pt_coord; + + nx = in->dims[0]; + ny = in->dims[1]; + nz = in->dims[2]; + + off[0] = 0; + off[1] = off[0] + 1; + off[2] = off[0] + (2 * nx); + off[3] = off[2] + 1; + off[4] = off[0] + ((2 * nx) * (2 * ny)); + off[5] = off[4] + 1; + off[6] = off[4] + (2 * nx); + off[7] = off[6] + 1; + + pt_coord = in->coord; + + origin = 0; + imax = (nx + 0) * 1 * (2 * 3); + jmax = (nx + 1) * (ny + 0) * (2 * 3); + + dx[0] = pt_coord[imax + 0] - pt_coord[origin + 0]; + dy[0] = pt_coord[imax + 1] - pt_coord[origin + 1]; + + dx[1] = pt_coord[jmax + 0] - pt_coord[origin + 0]; + dy[1] = pt_coord[jmax + 1] - pt_coord[origin + 1]; + + c = 0; dz = 0.0; + do { + active = (in->actnum == NULL) || (in->actnum[c] != 0); + + if (active) { + dz = vert_size(in, c, off); + } + + searching = ! (active && (fabs(dz) > 0.0)); + + c += 1; + } while (searching && (c < (nx * ny * nz))); + + assert (! searching); /* active && (fabs(dz) > 0) */ + + /* Compute vector triple product to distinguish left-handed (<0) + * from right-handed (>0) coordinate systems. */ + triple = dz * (dx[0]*dy[1] - dx[1]*dy[0]); + + assert (fabs(triple) > 0.0); + + return triple < 0.0; +} + + +/* ---------------------------------------------------------------------- */ +static void +reverse_face_nodes(struct processed_grid *out) +/* ---------------------------------------------------------------------- */ +{ + int f, t, *i, *j; + + for (f = 0; f < out->number_of_faces; f++) { + i = out->face_nodes + (out->face_ptr[f + 0] + 0); + j = out->face_nodes + (out->face_ptr[f + 1] - 1); + + assert (i <= j); + + while (i < j) { + t = *i; + *i = *j; + *j = t; + + i += 1; + j -= 1; + } + } } /*----------------------------------------------------------------- Public interface - */ +*/ void process_grdecl(const struct grdecl *in, double tolerance, struct processed_grid *out) { - int i,j,k; + struct grdecl g; - int nx = in->dims[0]; - int ny = in->dims[1]; - int nz = in->dims[2]; + size_t i; + int sign, error, left_handed; + int cellnum; - int nc = nx*ny*nz; - struct grdecl g; + int *actnum, *iptr; + int *global_cell_index; - int *actnum, *iptr, sign, error; - double *zcorn , *dptr; + double *zcorn; - int *intersections, *work, *plist, *ptr, *global_cell_index; + const size_t BIGNUM = 64; + const int nx = in->dims[0]; + const int ny = in->dims[1]; + const int nz = in->dims[2]; + const size_t nc = ((size_t) nx) * ((size_t) ny) * ((size_t) nz); - const int BIGNUM = 64; + /* internal work arrays */ + int *work; + int *plist; + int *intersections; - g.dims[0] = nx; - g.dims[1] = ny; - g.dims[2] = nz; - actnum = malloc (nc * sizeof *actnum); - zcorn = malloc (nc * 8 * sizeof *zcorn); - /* Permute actnum */ - iptr = actnum; - if (in->actnum != NULL) { - for (j=0; jactnum[i+nx*(j+ny*k)]; + + + /* -----------------------------------------------------------------*/ + /* Initialize output structure: + 1) allocate space for grid topology (which may need to be + increased) + 2) set Cartesian imensions + */ + out->m = (int) (BIGNUM / 3); + out->n = (int) 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->dimensions[0] = in->dims[0]; + out->dimensions[1] = in->dims[1]; + out->dimensions[2] = in->dims[2]; + out->number_of_faces = 0; + out->number_of_nodes = 0; + out->number_of_cells = 0; + + out->node_coordinates = NULL; + out->local_cell_index = malloc(nc * sizeof *out->local_cell_index); + + + + /* 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.*/ + + /* initialize grdecl structure "g" that will be processd by + * "finduniquepoints" */ + g.dims[0] = in->dims[0]; + g.dims[1] = in->dims[1]; + g.dims[2] = in->dims[2]; + + 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); + + g.coord = in->coord; + + + /* allocate space for cornerpoint numbers plus INT_MIN (INT_MAX) + * padding */ + plist = malloc(8 * (nc + ((size_t)nx)*((size_t)ny)) * sizeof *plist); + + finduniquepoints(&g, plist, tolerance, out); + + free (zcorn); + free (actnum); + + /* Determine if coordinate system is left handed or not. */ + left_handed = is_lefthanded(in); + if (left_handed) { + /* Reflect Y coordinates about XZ plane to create right-handed + * coordinate system whilst processing intersections. */ + for (i = 1; i < ((size_t) 3) * out->number_of_nodes; i += 3) { + out->node_coordinates[i] = -out->node_coordinates[i]; } - } } - } else { - for (i=0; i-2; sign = sign - 2){ - error = 0; - /* Ensure that zcorn is strictly nondecreasing in k-direction */ - for (j=0; j<2*ny; ++j){ - for (i=0; i<2*nx; ++i){ - for (k=0; k<2*nz-1; ++k){ - double z1 = sign*in->zcorn[i+2*nx*(j+2*ny*(k))]; - double z2 = sign*in->zcorn[i+2*nx*(j+2*ny*(k+1))]; + process_vertical_faces (0, &intersections, plist, work, out); + process_vertical_faces (1, &intersections, plist, work, out); + process_horizontal_faces ( &intersections, plist, out); - int c1 = i/2 + nx*(j/2 + ny*(k/2)); - int c2 = i/2 + nx*(j/2 + ny*((k+1)/2)); + free (plist); + free (work); - if (((in->actnum == NULL) || (in->actnum[c1] && in->actnum[c2])) && (z2 < z1)){ - fprintf(stderr, "\nZCORN should be strictly nondecreasing along pillars!\n"); + /* -----------------------------------------------------------------*/ + /* (re)allocate space for and compute coordinates of nodes that + * arise from intesecting cells (faults) */ + compute_intersection_coordinates(intersections, out); - error = 1; - goto end; - } + 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(nc * sizeof *global_cell_index); + cellnum = 0; + for (i = 0; i < nc; ++i) { + if (out->local_cell_index[i] != -1) { + global_cell_index[cellnum] = (int) i; + out->local_cell_index[i] = cellnum; + cellnum++; } - } } - end: - if (!error){ - break; + /* Remap out->face_neighbors */ + iptr = out->face_neighbors; + for (i = 0; i < ((size_t) 2) * out->number_of_faces; ++i, ++iptr) { + if (*iptr != -1){ + *iptr = out->local_cell_index[*iptr]; + } } - } - if (error){ - fprintf(stderr, "Attempt to reverse sign in ZCORN failed.\n" - "Grid definition may be broken\n"); - } - /* Permute zcorn */ - dptr = zcorn; - for (j=0; j<2*ny; ++j){ - for (i=0; i<2*nx; ++i){ - for (k=0; k<2*nz; ++k){ - *dptr++ = sign*in->zcorn[i+2*nx*(j+2*ny*k)]; - } + free(out->local_cell_index); + out->local_cell_index = global_cell_index; + + /* Reflect Y coordinate back to original position if left-handed + * coordinate system was detected and handled earlier. */ + if (left_handed) { + for (i = 1; i < ((size_t) 3) * out->number_of_nodes; i += 3) { + out->node_coordinates[i] = -out->node_coordinates[i]; + } } - } - - - - g.zcorn = zcorn; - g.coord = in->coord; - - - /* Code below assumes k index runs fastests, ie. that dimensions of - table is permuted to (dims[2], dims[0], dims[1]) */ - - 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->dimensions[0] = g.dims[0]; - out->dimensions[1] = g.dims[1]; - out->dimensions[2] = g.dims[2]; - 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); - - - /* internal */ - intersections = malloc(BIGNUM* sizeof(*intersections)); - - /* internal */ - work = malloc(2* (2*nz+2)* sizeof(*work)); - for(i=0; i<4*(nz+1); ++i) work[i] = -1; - - /* Allocate space for cornerpoint numbers plus INT_MIN (INT_MAX) padding */ - plist = malloc( 4*nx*ny*(2*nz+2) * sizeof *plist); - - - /* Do actual work here:*/ - - finduniquepoints(&g, plist, tolerance, out); - - free (zcorn); - free (actnum); - - 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); - - compute_intersection_coordinates(intersections, out); - - free (intersections); - - - /* Convert to local cell numbers in face_neighbors */ - ptr=out->face_neighbors;; - for (i=0; inumber_of_faces*2; ++i, ++ptr){ - if (*ptr != -1){ - *ptr = out->local_cell_index[*ptr]; + /* 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 < ((size_t) 3) * out->number_of_nodes; i += 3) + out->node_coordinates[i] *= sign; } - } - /* Invert global-to-local map */ - global_cell_index = malloc(out->number_of_cells * - sizeof (*global_cell_index)); - for (i=0; ilocal_cell_index[i]!=-1){ - global_cell_index[out->local_cell_index[i]] = i; + /* If an odd number of coordinate reflections were applied, the + * processing routines--especially facetopology()--will produce + * node orderings that lead to normals pointing from 2 to 1. + * Reverse nodes to reestablish expected normal direction (and + * positive cell volumes). */ + if (left_handed ^ (sign == -1)) { + reverse_face_nodes(out); } - } - free(out->local_cell_index); - out->local_cell_index = global_cell_index; - - /* HACK: undo sign reversal in ZCORN prepprocessing */ - 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/opm/core/grid/cpgpreprocess/preprocess.h b/opm/core/grid/cpgpreprocess/preprocess.h index e0453cc4..9592ce92 100644 --- a/opm/core/grid/cpgpreprocess/preprocess.h +++ b/opm/core/grid/cpgpreprocess/preprocess.h @@ -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, @@ -53,26 +53,26 @@ extern "C" { /* Output structure holding grid topology */ struct processed_grid{ - int m; - int n; + int m; /** Upper bound on "number_of_faces" */ + int n; /** Upper bound on "number_of_nodes" */ int dimensions[3]; /* Cartesian dimension */ int number_of_faces; int *face_nodes; /* Nodes numbers of each face sequentially. */ int *face_ptr; /* Start position for each face in face_nodes. */ int *face_neighbors; /* Global cell numbers. 2 ints per face sequentially */ enum face_tag *face_tag; - - int number_of_nodes; - int number_of_nodes_on_pillars; + + int number_of_nodes; + int number_of_nodes_on_pillars; /** Total number of unique cell vertices that lie on pillars. */ double *node_coordinates; /* 3 doubles per node, sequentially */ - + int number_of_cells; /* number of active cells */ int *local_cell_index; /* Global to local map */ }; - void process_grdecl (const struct grdecl *g, - double tol, + void process_grdecl (const struct grdecl *g, + double tol, struct processed_grid *out); void free_processed_grid(struct processed_grid *g); diff --git a/opm/core/grid/cpgpreprocess/processgrid.c b/opm/core/grid/cpgpreprocess/processgrid.c index f15ee344..9ef848ef 100644 --- a/opm/core/grid/cpgpreprocess/processgrid.c +++ b/opm/core/grid/cpgpreprocess/processgrid.c @@ -13,255 +13,453 @@ //===========================================================================*/ /* -Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. -Copyright 2009, 2010 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. */ /* Mex gateway by Jostein R. Natvig, SINTEF ICT. */ -#include -#include -#include #include +#include +#include +#include + #include #include "preprocess.h" #include "mxgrdecl.h" -void fill_grid(mxArray **out, struct processed_grid *grid) +/* ---------------------------------------------------------------------- */ +static mxArray * +allocate_nodes(size_t nnodes) +/* ---------------------------------------------------------------------- */ { - const char *names[] = {"nodes", "faces", "cells", "cartDims", "type"}; - mxArray *G = mxCreateStructMatrix(1,1,sizeof names / sizeof names[0],names); + size_t nflds; + const char *fields[] = { "num", "coords" }; - int i,j; - double *ptr; + mxArray *nodes, *num, *coords; + nflds = sizeof(fields) / sizeof(fields[0]); - /* nodes */ - const char *n2[] = {"num", "coords"}; - mxArray *nodes = mxCreateStructMatrix(1,1,sizeof n2 / sizeof n2[0],n2); - mxSetField(nodes, 0, "num", mxCreateDoubleScalar(grid->number_of_nodes)); + nodes = mxCreateStructMatrix(1, 1, nflds, fields); - mxArray *nodecoords = mxCreateDoubleMatrix(grid->number_of_nodes, 3, mxREAL); - ptr = mxGetPr(nodecoords); + num = mxCreateDoubleScalar(nnodes); + coords = mxCreateDoubleMatrix(nnodes, 3, mxREAL); - for (j=0;j<3;++j) - { - for(i=0; inumber_of_nodes; ++i) - { - ptr[i+grid->number_of_nodes*j] = grid->node_coordinates[3*i+j]; + 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; } - } - mxSetField(nodes, 0, "coords", nodecoords); - mxSetField(G, 0, "nodes", nodes); - - /* faces */ - const char *n3[] = {"num", "neighbors", "nodes", "numNodes", "nodePos", "tag"}; - mxArray *faces = mxCreateStructMatrix(1,1,sizeof n3 / sizeof n3[0], n3); - - - mxSetField(faces, 0, "num", mxCreateDoubleScalar(grid->number_of_faces)); - - mxArray *faceneighbors = mxCreateNumericMatrix(grid->number_of_faces, 2, - mxINT32_CLASS, mxREAL); - { - int *iptr = mxGetData(faceneighbors); - for(j=0; j<2; ++j) - { - for (i=0; inumber_of_faces; ++i) - { - int ix = grid->face_neighbors[2*i+j]; - if (ix == -1) - { - iptr[i+grid->number_of_faces*j] = 0; - } - else - { - iptr[i+grid->number_of_faces*j] = ix+1; - } - } - } - } - mxSetField(faces, 0, "neighbors", faceneighbors); - mxArray *numnodes = mxCreateNumericMatrix(grid->number_of_faces, 1, - mxINT32_CLASS, mxREAL); - mxArray *nodepos = mxCreateNumericMatrix(grid->number_of_faces+1, 1, - mxINT32_CLASS, mxREAL); - { - int *iptr1 = mxGetData(numnodes); - int *iptr2 = mxGetData(nodepos); - iptr2[0] = 1; - for (i=0; inumber_of_faces; ++i) - { - iptr1[i] = grid->face_ptr[i+1] - grid->face_ptr[i]; - iptr2[i+1] = iptr2[i] + iptr1[i]; - } - } - mxSetField(faces, 0, "numNodes", numnodes); - mxSetField(faces, 0, "nodePos", nodepos); - - mxArray *tags = mxCreateNumericMatrix(grid->number_of_faces, 1, - mxINT32_CLASS, mxREAL); - { - int *iptr = mxGetData(tags); - for (i = 0; i < grid->number_of_faces; ++i) - { - iptr[i] = grid->face_tag[i] + 1; - } - } - mxSetField(faces, 0, "tag", tags); - - - const char *n4[] = {"num", "faces", "facePos", "indexMap"}; - mxArray *cells = mxCreateStructMatrix(1,1,sizeof n4 / sizeof n4[0], n4); - - mxSetField(cells, 0, "num", mxCreateDoubleScalar(grid->number_of_cells)); - - mxArray *map = mxCreateNumericMatrix(grid->number_of_cells, 1, - mxINT32_CLASS, mxREAL); - { - int *iptr = mxGetData(map); - for(i=0; inumber_of_cells; ++i) - { - iptr[i] = grid->local_cell_index[i]+1; - } - } - mxSetField(cells, 0, "indexMap", map); - - mxArray *facepos = mxCreateNumericMatrix(grid->number_of_cells+1, 1, - mxINT32_CLASS, mxREAL); - { - int *iptr = mxGetData(facepos); - for(i=0; inumber_of_cells+1; ++i) - { - iptr[i] = 0; - } - for (i=0; i<2*grid->number_of_faces; ++i) - { - int c=grid->face_neighbors[i]; - if(c != -1) - { - iptr[c+1]++; - } - } - iptr[0] = 1; - for(i=0; inumber_of_cells; ++i) - { - iptr[i+1] += iptr[i]; - } - } - mxSetField(cells, 0, "facePos", facepos); - - - - int *counter = calloc(grid->number_of_cells, sizeof(*counter)); - int num_half_faces = 0; - { - int *iptr = mxGetData(facepos); - for(i=0; inumber_of_cells; ++i) - { - counter[i] = num_half_faces; - num_half_faces += iptr[i+1]-iptr[i]; - } - } - - mxArray *cellfaces = mxCreateNumericMatrix(num_half_faces, 1, - mxINT32_CLASS, mxREAL); - { - int *iptr = mxGetData(cellfaces); - for (i=0; inumber_of_faces; ++i) - { - int c1 = grid->face_neighbors[2*i]; - int c2 = grid->face_neighbors[2*i+1]; - if(c1 != -1) iptr[counter[c1]++] = i+1; - if(c2 != -1) iptr[counter[c2]++] = i+1; - } - } - mxSetField(cells, 0, "faces", cellfaces); - - mxSetField(G, 0, "cells", cells); - - int n = grid->face_ptr[grid->number_of_faces]; - - mxArray *facenodes = mxCreateNumericMatrix(n, 1, - mxINT32_CLASS, mxREAL); - { - int *iptr = mxGetData(facenodes); - for (i=0; iface_nodes[i]+1; - } - } - mxSetField(faces, 0, "nodes", facenodes); - - mxSetField(G, 0, "faces", faces); - - - free(counter); - - /* cartDims */ - mxArray *cartDims = mxCreateDoubleMatrix(1, 3, mxREAL); - ptr = mxGetPr(cartDims); - ptr[0] = grid->dimensions[0]; - ptr[1] = grid->dimensions[1]; - ptr[2] = grid->dimensions[2]; - - mxSetField(G, 0, "cartDims", cartDims); - - /* type */ - mxArray *type = mxCreateCellMatrix(1, 1); - mxSetCell(type, 0, mxCreateString("processgrid")); - - mxSetField(G, 0, "type", type); - - out[0] = G; + return nodes; } -/* Gateway routine for Matlab mex function. */ -/*-------------------------------------------------------*/ -void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +/* ---------------------------------------------------------------------- */ +static void +fill_nodes(mxArray *nodes, struct processed_grid *grid) +/* ---------------------------------------------------------------------- */ { - /* Set up data passed from Matlab */ - struct grdecl g; - struct processed_grid out; - double tolerance = 0.0; + size_t i, j, nnodes; + double *coords; - mx_init_grdecl(&g, prhs[0]); - if (nrhs == 2) - { - tolerance = mxGetScalar (prhs[1]); - } + nnodes = grid->number_of_nodes; + coords = mxGetPr(mxGetField(nodes, 0, "coords")); - process_grdecl(&g, tolerance, &out); - - - if (plhs >0) - { - /* write to matlab struct */ - fill_grid(plhs, &out); - } - - - /* Free whatever was allocated in initGrdecl. */ - free_processed_grid(&out); + for (j = 0; j < 3; ++j) { + for (i = 0; i < nnodes; ++i) { + *coords++ = grid->node_coordinates[3*i + j]; + } + } } + + +/* ---------------------------------------------------------------------- */ +static mxArray * +allocate_faces(size_t nf, size_t nfacenodes) +/* ---------------------------------------------------------------------- */ +{ + size_t nflds; + const char *fields[] = { "num", "neighbors", "nodePos", "nodes", "tag" }; + + mxArray *faces, *num, *neighbors, *nodePos, *nodes, *tag; + + nflds = sizeof(fields) / sizeof(fields[0]); + + 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); + + 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; + } + + return faces; +} + + +/* ---------------------------------------------------------------------- */ +static void +fill_faces(mxArray *faces, struct processed_grid *grid) +/* ---------------------------------------------------------------------- */ +{ + size_t i, f, nf, nfn; + + int *pi; + + 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.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.tag */ + pi = mxGetData(mxGetField(faces, 0, "tag")); + for (f = 0; f < nf; f++) { pi[f] = grid->face_tag[f] + 1; } +} + + +/* ---------------------------------------------------------------------- */ +static size_t +count_halffaces(size_t nf, const int *neighbors) +/* ---------------------------------------------------------------------- */ +{ + int c1, c2; + size_t nhf, f; + + for (f = nhf = 0; f < nf; f++) { + c1 = neighbors[2*f + 0]; + c2 = neighbors[2*f + 1]; + + nhf += c1 >= 0; + nhf += c2 >= 0; + } + + return nhf; +} + + +/* ---------------------------------------------------------------------- */ +static mxArray * +allocate_cells(size_t nc, size_t ncf) +/* ---------------------------------------------------------------------- */ +{ + size_t nflds; + const char *fields[] = { "num", "facePos", "faces", "indexMap" }; + + mxArray *cells, *num, *facePos, *faces, *indexMap; + + nflds = sizeof(fields) / sizeof(fields[0]); + + 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); + + 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; + } + + return cells; +} + + +/* ---------------------------------------------------------------------- */ +static void +fill_cells(mxArray *cells, struct processed_grid *grid) +/* ---------------------------------------------------------------------- */ +{ + size_t c, nc, f, nf, i; + + int c1, c2, cf_tag, nhf; + int *pi1, *pi2; + + 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; } + + /* 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 (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; } +} + + +/* ---------------------------------------------------------------------- */ +static mxArray * +allocate_grid(struct processed_grid *grid, const char *func) +/* ---------------------------------------------------------------------- */ +{ + size_t nflds, nhf; + const char *fields[] = { "nodes", "faces", "cells", + "type", "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); + + 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); + + 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); } + + G = NULL; + } + + return G; +} + + +/* ---------------------------------------------------------------------- */ +static void +fill_grid(mxArray *G, struct processed_grid *grid) +/* ---------------------------------------------------------------------- */ +{ + double *pr; + + 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); +} + + +/* ---------------------------------------------------------------------- */ +static int +args_ok(int nlhs, int nrhs, const mxArray *prhs[]) +/* ---------------------------------------------------------------------- */ +{ + int ok; + + ok = (nlhs == 1) && ((nrhs == 1) || (nrhs == 2)); + + 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); + + if (ok && (nrhs == 2)) { + ok = mxIsDouble(prhs[1]) && (mxGetNumberOfElements(prhs[1]) == 1); + } + + return ok; +} + + +/* ---------------------------------------------------------------------- */ +static double +define_tolerance(int nrhs, const mxArray *prhs[]) +/* ---------------------------------------------------------------------- */ +{ + double tol; + + tol = 0.0; + + if (nrhs == 2) { + tol = mxGetScalar(prhs[1]); + } + + 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; + + if (args_ok(nlhs, nrhs, prhs)) { + mx_init_grdecl(&grdecl, prhs[0]); + tolerance = define_tolerance(nrhs, prhs); + + process_grdecl(&grdecl, tolerance, &g); + + plhs[0] = allocate_grid(&g, mexFunctionName()); + + 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 { + 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'", + mexFunctionName(), mexFunctionName()); + mexErrMsgTxt(errmsg); + } +} + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ diff --git a/opm/core/grid/cpgpreprocess/processgrid.m b/opm/core/grid/cpgpreprocess/processgrid.m index af3a1bcf..372cdbdc 100644 --- a/opm/core/grid/cpgpreprocess/processgrid.m +++ b/opm/core/grid/cpgpreprocess/processgrid.m @@ -1,12 +1,14 @@ -function varargout = processgrid(varargin) +function G = processgrid(varargin) %Compute grid topology and geometry from pillar grid description. % % SYNOPSIS: -% G = processGRDECL(grdecl) +% G = processgrid(grdecl) +% G = processgrid(grdecl,ztol) % % PARAMETERS: % grdecl - Raw pillar grid structure, as defined by function % 'readGRDECL', with fields COORDS, ZCORN and, possibly, ACTNUM. +% ztol - tolerance for unique points % % RETURNS: % G - Valid grid definition containing connectivity, cell @@ -32,10 +34,33 @@ function varargout = processgrid(varargin) % $Date$ % $Revision$ -% Build MEX edition of same. -% -buildmex CFLAGS='$CFLAGS -Wall -fPIC' processgrid.c preprocess.c ... - uniquepoints.c facetopology.c sparsetable.c mxgrdecl.c + G = processgrid_mex(varargin{:}); + G.griddim = 3; + G = splitDisconnectedGrid(G, false); +end -% Call MEX edition. -[varargout{1:nargout}] = processgrid(varargin{:}); +function G = splitDisconnectedGrid(G, verbose) + % Check if grid is connected + ix = all(G.faces.neighbors~=0, 2); + I = [G.faces.neighbors(ix,1);G.faces.neighbors(ix,2)]; + J = [G.faces.neighbors(ix,2);G.faces.neighbors(ix,1)]; + N = double(max(G.faces.neighbors(:))); + A = sparse(double(I),double(J),1,N,N)+speye(N); + clear ix I J + [a,b,c,d]=dmperm(A); %#ok + clear A b d + if numel(c) > 2, + dispif(verbose, '\nGrid has %d disconnected components\n', ... + numel(c)- 1); + % Partition grid into connected subgrids + for i = 1:numel(c) - 1, + g(i) = extractSubgrid(G, a(c(i):c(i+1)-1)); %#ok + sz(i) = g(i).cells.num; %#ok + g(i).cartDims = G.cartDims; %#ok + end + + % Return largest (in number of cells) grid first + [i,i] = sort(-sz); %#ok + G = g(i); + end +end diff --git a/opm/core/grid/cpgpreprocess/readvector.cpp b/opm/core/grid/cpgpreprocess/readvector.cpp deleted file mode 100644 index 7b5adc5a..00000000 --- a/opm/core/grid/cpgpreprocess/readvector.cpp +++ /dev/null @@ -1,282 +0,0 @@ -/*=========================================================================== -// -// Author: Jostein R. Natvig -// -//==========================================================================*/ - - -/* - Copyright 2011 SINTEF ICT, Applied Mathematics. -*/ - - -#include -#include -#include -#include -#include - -#include - -#include "readvector.hpp" - -struct ParserState { - int error; -}; - -template -static void parse_next_token (FILE *fp, struct ParserState *s, - T *v, int *count); -static char *get_token_string (FILE *fp, struct ParserState *s, char *buf); -static void skip_line (FILE *fp, struct ParserState *s); -static char *convert_string (int *value, char *str); -static char *convert_string (double *value, char *str); -static void split_token_string(char *str, char *rstr, char *vstr); -static int get_repeat_count (char *rstr, struct ParserState *s); -template -static void get_value (char *vstr, struct ParserState *s, T *v); - -/* ------------------------------------------------------------------ */ - -static const int bufsz = 1024; - -/* ------------------------------------------------------------------ */ - -template -static void readvector(FILE *fp, struct ParserState *s, std::vector& v) -{ - T value; - int count, i; - - count = 1; - s->error = 0; - while (count > 0){ - - parse_next_token(fp, s, &value, &count); - for (i=0; ierror) { - assert(0); - } -} - - -/* ------------------------------------------------------------------ */ -template -static void parse_next_token(FILE *fp, struct ParserState *s, T *v, int *r) -{ - char str[bufsz]; - char rstr[bufsz]; - char vstr[bufsz]; - - str[0]='\0'; - while (strlen(str)==0) - { - get_token_string(fp, s, str); - } - - if (str[0] == '/') - { - /* end-of-vector marker */ - *r = 0; - } - else - { - split_token_string(str, rstr, vstr); - *r = get_repeat_count(rstr, s); - get_value(vstr, s, v); - - if (s->error) - { - /* signal abort to caller */ - *r = 0; - } - } -} - -/* ------------------------------------------------------------------ */ - -/* split string 'rrrr*vvvv' in strings 'rrrr' and 'vvvv' */ -static void split_token_string(char *str, char *rstr, char *vstr) -{ - char *ptr; - - ptr=strchr(str, '*'); - if ((ptr != NULL) && (ptr!=str)) { - while(str!=ptr) - { - *rstr++=*str++; - } - *rstr='\0'; - - str++; - } - else - { - rstr[0]='\0'; - } - - while((*vstr++=*str++) ) - { - continue ; - } - *vstr='\0'; -} - -/* ------------------------------------------------------------------ */ - -static int get_repeat_count(char *rstr, struct ParserState *s) -{ - int r; - char *ptr; - - if (strlen(rstr)>0) - { - ptr=convert_string(&r, rstr); - - if ((ptr==rstr) || (strlen(ptr)>0) || (r<1)) - { - s->error = 1; - } - } - else - { - r = 1; - } - - return r; -} - -/* ------------------------------------------------------------------ */ -template -static void get_value(char *vstr, struct ParserState *s, T *v) -{ - char *ptr; - - if (strlen(vstr)>0) - { - /* Convert FORTRAN style floats 3.13D-6 to IEEE */ - for(ptr=vstr; *ptr; ++ptr) - { - if (*ptr=='D') - { - *ptr='E'; - } - } - - ptr = convert_string(v, vstr); - - if ((ptr==vstr) || (strlen(ptr)>0)) - { - s->error = 1; - } - } - else - { - s->error = 1; - } -} - - -/* ------------------------------------------------------------------ */ - -static char *get_token_string(FILE *fp, struct ParserState *s, char *str) -{ - char *ptr = str; - int c; - - /* Skip leading blanks */ - while((c = fgetc(fp)) != EOF && isspace(c)) - { - ; - } - - *ptr++ = c; - - /* Catch end marker */ - if (c == '/') { - skip_line(fp, s); - *ptr++ = '\0'; - return str; - } - - while((c = fgetc(fp)) != EOF && !isspace(c)) - { - *ptr++ = c; - - /* Break and skip rest of line if '--' if encountered */ - if (*(ptr-1) == '-' && ptr-str>1 && *(ptr-2) == '-'){ - ptr = ptr - 2; - skip_line(fp, s); - break; - } - - /* If end marker is encontered, push character back onto stream. */ - if (c=='/') { - ungetc(c, fp); - ptr--; - break; - } - - assert(ptr-str < bufsz); - } - - *ptr='\0'; - return str; -} - -/* ------------------------------------------------------------------ */ - -static void skip_line(FILE *fp, struct ParserState *s) -{ - static_cast(s); - int c; - while((c = fgetc(fp))!=EOF && c != '\n') { - ; - } -} - -/* ------------------------------------------------------------------ */ - -// int version -static char *convert_string(int *value, char *str) -{ - char *q; - *value = strtol(str, &q, 10); - return q; -} - -/* ------------------------------------------------------------------ */ - -// double version -static char *convert_string(double *value, char *str) -{ - char *q; - *value = strtod(str, &q); - return q; -} - -/* ------------------------------------------------------------------ */ -template -void read_vector_from_file(const std::string& fn, std::vector& v) -{ - FILE *fp = fopen(fn.c_str(), "r"); - struct ParserState s = { 0 }; - readvector(fp, &s, v); - fclose(fp); -} - - -void read_vector_from_file(const std::string& fn, std::vector& v) -{ - read_vector_from_file(fn, v); -} - -void read_vector_from_file(const std::string& fn, std::vector& v) -{ - read_vector_from_file(fn, v); -} diff --git a/opm/core/grid/cpgpreprocess/readvector.hpp b/opm/core/grid/cpgpreprocess/readvector.hpp deleted file mode 100644 index 04cdb977..00000000 --- a/opm/core/grid/cpgpreprocess/readvector.hpp +++ /dev/null @@ -1,24 +0,0 @@ -/*=========================================================================== -// -// File: readvector.hpp -// -// Created: 2011-11-30 09:35:14+0100 -// -// Author: Jostein R. Natvig -// -//==========================================================================*/ - - -/* - Copyright 2011 SINTEF ICT, Applied Mathematics. -*/ - -#ifndef READVECTOR_HPP_HEADER -#define READVECTOR_HPP_HEADER - -#include - -void read_vector_from_file(const std::string&, std::vector& v); -void read_vector_from_file(const std::string&, std::vector& v); - -#endif /* READVECTOR_HPP_HEADER */ diff --git a/opm/core/grid/cpgpreprocess/sparsetable.c b/opm/core/grid/cpgpreprocess/sparsetable.c deleted file mode 100644 index 9ef36bab..00000000 --- a/opm/core/grid/cpgpreprocess/sparsetable.c +++ /dev/null @@ -1,100 +0,0 @@ -/*=========================================================================== -// -// File: sparsetable.c -// -// Created: Fri Jun 19 08:48:04 2009 -// -// Author: Jostein R. Natvig -// -// $Date$ -// -// $Revision$ -// -//==========================================================================*/ - -/* -Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. -Copyright 2009, 2010 Statoil ASA. - -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 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 . -*/ - -#include -#include -#include -#include -#include - - -#include "sparsetable.h" - -void free_sparse_table (sparse_table_t *tab) -{ - if(tab->ptr) free(tab->ptr); - if(tab->data) free(tab->data); - - free(tab); -} - -sparse_table_t *malloc_sparse_table(int m, int n, int datasz) -{ - size_t alloc_sz; - sparse_table_t *tab = malloc(sizeof *tab); - tab->m = m; - tab->n = n; - tab->position = 0; - if (!(tab->ptr = malloc((m+1) * sizeof (*tab->ptr)))){ - fprintf(stderr, "Could not allocate space for sparse ptr\n"); - free_sparse_table(tab); - return NULL; - } - - alloc_sz = datasz; - alloc_sz *= n; - if(!(tab->data = malloc(alloc_sz))){ - fprintf(stderr, "Could not allocate space for sparse data(%d)\n", n); - free_sparse_table(tab); - return NULL; - } - - return tab; -} - -sparse_table_t *realloc_sparse_table(sparse_table_t *tab, int m, int n, int datasz) -{ - void *p = realloc(tab->ptr, (m+1) * sizeof (*tab->ptr)); - if (p){ - tab->ptr = p; - }else{ - fprintf(stderr, "Could not reallocate space for sparse ptr\n"); - free_sparse_table(tab); - return NULL; - } - - p = realloc(tab->data, n * datasz); - if(p){ - tab->data = p; - }else{ - fprintf(stderr, "Could not reallocate space for sparse data(%d)\n", n); - free_sparse_table(tab); - return NULL; - } - - tab->m = m; - tab->n = n; - - return tab; -} diff --git a/opm/core/grid/cpgpreprocess/sparsetable.h b/opm/core/grid/cpgpreprocess/sparsetable.h deleted file mode 100644 index f28df362..00000000 --- a/opm/core/grid/cpgpreprocess/sparsetable.h +++ /dev/null @@ -1,53 +0,0 @@ -/*=========================================================================== -// -// File: sparsetable.h -// -// Created: Fri Jun 19 08:47:45 2009 -// -// Author: Jostein R. Natvig -// -// $Date$ -// -// $Revision$ -// -//===========================================================================*/ - -/* -Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. -Copyright 2009, 2010 Statoil ASA. - -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 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 . -*/ - -#ifndef OPENRS_SPARSETABLE_HEADER -#define OPENRS_SPARSETABLE_HEADER - -typedef struct{ - int m; /* number of rows */ - int *ptr; /* row pointer of size m+1 */ - int position; /* first position in ptr that is not filled. */ - - int n; /* size of data */ - void *data; /* sparse table data */ - - -} sparse_table_t; - -void free_sparse_table (sparse_table_t *tab); -sparse_table_t *malloc_sparse_table (int m, int n, int datasz); -sparse_table_t *realloc_sparse_table (sparse_table_t *tab, int m, int n, int datasz); - -#endif /* OPENRS_SPARSETABLE_HEADER */ diff --git a/opm/core/grid/cpgpreprocess/uniquepoints.c b/opm/core/grid/cpgpreprocess/uniquepoints.c index 1513e858..4e6779d9 100644 --- a/opm/core/grid/cpgpreprocess/uniquepoints.c +++ b/opm/core/grid/cpgpreprocess/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, @@ -33,326 +33,350 @@ */ #include -#include -#include -#include -#include #include +#include +#include #include +#include +#include -#include "sparsetable.h" #include "preprocess.h" #include "uniquepoints.h" -#define min(i,j) ((i)<(j) ? (i) : (j)) -#define max(i,j) ((i)>(j) ? (i) : (j)) -#define overlap(a1,a2,b1,b2) max(a1,b1) < min(a2,b2) - +#define MIN(i,j) (((i) < (j)) ? (i) : (j)) +#define MAX(i,j) (((i) > (j)) ? (i) : (j)) /*----------------------------------------------------------------- - Compare function passed to qsort -*/ + Compare function passed to qsortx */ static int compare(const void *a, const void *b) { - if (*(const double*)a < *(const double*) b) return -1; - else return 1; + const double a0 = *(const double*) a; + const double b0 = *(const double*) b; + + /* { -1, a < b + * compare(a,b) = { 0, a = b + * { 1, a > b */ + return (a0 > b0) - (a0 < b0); } /*----------------------------------------------------------------- - 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, pos = 0; - double val; + int i; + int pos; + double val; - assert (!(tolerance < 0.0)); + assert (!(tolerance < 0.0)); - if (n<1) return 0; + if (n<1) return 0; + pos = 0; + val = list[pos++];/* Keep first value */ - 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[]) +/* ---------------------------------------------------------------------- */ +static void +vector_positions(const int dims[3] , + const int i , + const int j , + size_t start[4]) +/* ---------------------------------------------------------------------- */ { - - 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); + size_t im, ip, jm, jp; + + im = MAX(1, i ) - 1; + jm = MAX(1, j ) - 1; + ip = MIN(dims[0], i+1) - 1; + jp = MIN(dims[1], j+1) - 1; + + start[ 0 ] = dims[2] * (im + dims[0]*jm); + start[ 1 ] = dims[2] * (im + dims[0]*jp); + start[ 2 ] = dims[2] * (ip + dims[0]*jm); + start[ 3 ] = 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 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); + size_t p, start[4]; + + vector_positions(dims, i, j, start); + + for (p = 0; p < 4; p++) { + v[p] = field + start[p]; + } +} + + +/*----------------------------------------------------------------- + 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[]) +{ + size_t p, start[4]; + + vector_positions(dims, i, j, start); + + for (p = 0; p < 4; p++) { + v[p] = field + start[p]; + } } /*----------------------------------------------------------------- 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; + + if (fabs(coord[5] - coord[2]) > 0) { + + a = (pt[2] - coord[2]) / (coord[5] - coord[2]); + + } else { + + 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) + { - int nx = out->dimensions[0]; - int ny = out->dimensions[1]; - int nz = out->dimensions[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); - sparse_table_t *ztab = malloc_sparse_table(npillars, - npillarpoints, - sizeof(double)); + 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]; + + + /* 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); - int nc = g->dims[0]*g->dims[1]*g->dims[2]; - double *zlist = ztab->data; /* casting void* to double* */ - int *zptr = ztab->ptr; + int i,j,k; - int i,j,k; + int d1[3]; + int len = 0; + double *zout = zlist; + int pos = 0; + 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; + const double *coord = g->coord; - const double *coord = g->coord; - double *pt; - int *p; + 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){ + /* 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){ - const int *a[4]; - const double *z[4]; - - /* 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); + /* 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); - len = createSortedList( zout, d1[2], 4, z, a); - len = uniquify (len, zout, tolerance); + len = createSortedList( zout, d1[2], 4, z, a); + len = uniquify (len, zout, tolerance); - /* 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 */ - int pix = (i+1)/2 + (g->dims[0]+1)*((j+1)/2); - - /* cell column position */ - int 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 */ - int zix = 2*g->dims[2]*(i+2*g->dims[0]*j); - - const int *a = g->actnum + cix; - const double *z = g->zcorn + zix; + /* pillar index */ + pix = (i+1)/2 + (g->dims[0]+1)*((j+1)/2); - if (!assignPointNumbers(zptr[pix], zptr[pix+1], zlist, - 2*g->dims[2], z, a, p, tolerance)){ - fprintf(stderr, "Something went wrong in assignPointNumbers"); - return 0; - } + /* cell column position */ + cix = g->dims[2]*((i/2) + (j/2)*g->dims[0]); - p += 2 + 2*g->dims[2]; + /* 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_sparse_table(ztab); - - - return 1; + return 1; } /* Local Variables: */ /* c-basic-offset:4 */ /* End: */ - diff --git a/opm/core/grid/cpgpreprocess/uniquepoints.h b/opm/core/grid/cpgpreprocess/uniquepoints.h index ac3574bf..b9d698ac 100644 --- a/opm/core/grid/cpgpreprocess/uniquepoints.h +++ b/opm/core/grid/cpgpreprocess/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 */ diff --git a/opm/core/linalg/sparse_sys.h b/opm/core/linalg/sparse_sys.h index 1485c2f4..ce34bc51 100644 --- a/opm/core/linalg/sparse_sys.h +++ b/opm/core/linalg/sparse_sys.h @@ -37,13 +37,13 @@ extern "C" { */ struct CSRMatrix { - size_t m; /** Number of rows */ - size_t nnz; /** Number of structurally non-zero elements */ + size_t m; /**< Number of rows */ + size_t nnz; /**< Number of structurally non-zero elements */ - int *ia; /** Row pointers */ - int *ja; /** Column indices */ + int *ia; /**< Row pointers */ + int *ja; /**< Column indices */ - double *sa; /** Matrix elements */ + double *sa; /**< Matrix elements */ }; diff --git a/tests/Makefile.am b/tests/Makefile.am index ae4ee69d..1ec5dd26 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -52,8 +52,6 @@ test_lapack_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) test_readpolymer_SOURCES = test_readpolymer.cpp test_read_vag_SOURCES = test_read_vag.cpp -test_readvector_SOURCES = test_readvector.cpp - test_sf2p_SOURCES = test_sf2p.cpp test_writeVtkData_SOURCES = test_writeVtkData.cpp @@ -77,4 +75,4 @@ endif if HAVE_ERT noinst_PROGRAMS += test_ert test_ert_SOURCES = test_ert.cpp -endif \ No newline at end of file +endif diff --git a/tests/test_readvector.cpp b/tests/test_readvector.cpp deleted file mode 100644 index 56c58015..00000000 --- a/tests/test_readvector.cpp +++ /dev/null @@ -1,67 +0,0 @@ -#include -#include -#include - -#include -#include -#include - -#include - -static struct UnstructuredGrid* -read_grid(const std::string& dir) -{ - std::string fn; - fn = dir + '/' + "zcorn.txt"; - - std::vector zcorn; - read_vector_from_file(fn, zcorn); - - fn = dir + '/' + "coord.txt"; - ::std::vector coord; - read_vector_from_file(fn, coord); - - fn = dir + '/' + "actnum.txt"; - std::vector actnum; - read_vector_from_file(fn, actnum); - - fn = dir + '/' + "dimens.txt"; - ::std::vector dimens; - read_vector_from_file(fn, dimens); - - struct grdecl grdecl; - grdecl.zcorn = &zcorn[0]; - grdecl.coord = &coord[0]; - grdecl.actnum = &actnum[0]; - - grdecl.dims[0] = dimens[0]; - grdecl.dims[1] = dimens[1]; - grdecl.dims[2] = dimens[2]; - - struct UnstructuredGrid *g = create_grid_cornerpoint(&grdecl, 0.0); - - double vol = 0.0; - for (int c = 0; c < g->number_of_cells; c++) { - vol += g->cell_volumes[c]; - } - std::cout << "Sum volumes = " << vol << '\n'; - - for (int c = 0, i = 0; c < g->number_of_cells; c++) { - for (; i < g->cell_facepos[c + 1]; i++) { - std::cout << "(c,i) = (" << c << "," << g->cell_facetag[i] << ")\n"; - } - } - - return g; -} - -int main() -{ - struct UnstructuredGrid *g; - - g = read_grid(std::string("example")); - - destroy_grid(g); - - return 0; -}