diff --git a/opm/core/fluid/BlackoilPropertiesBasic.hpp b/opm/core/fluid/BlackoilPropertiesBasic.hpp index 44a76013b..cc34a3cfe 100644 --- a/opm/core/fluid/BlackoilPropertiesBasic.hpp +++ b/opm/core/fluid/BlackoilPropertiesBasic.hpp @@ -25,6 +25,7 @@ #include #include #include +#include namespace Opm { diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index 981ccd9ba..f0de92b07 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -42,12 +42,15 @@ namespace Opm // to pressure in first perforation cell. for (int w = 0; w < nw; ++w) { const WellControls* ctrl = wells->ctrls[w]; - if (ctrl->type[ctrl->current] == BHP) { - bhp_[w] = ctrl->target[ctrl->current]; - } else { + + if ((ctrl->current < 0) || // SHUT + (ctrl->type[ctrl->current] != BHP)) { const int cell = wells->well_cells[wells->well_connpos[w]]; bhp_[w] = state.pressure()[cell]; } + else { + bhp_[w] = ctrl->target[ctrl->current]; + } } perfrates_.resize(wells->well_connpos[nw], 0.0); perfpress_.resize(wells->well_connpos[nw], -1e100); diff --git a/opm/core/transport/reorder/DGBasis.cpp b/opm/core/transport/reorder/DGBasis.cpp index 8527a1caa..8ffb95ec4 100644 --- a/opm/core/transport/reorder/DGBasis.cpp +++ b/opm/core/transport/reorder/DGBasis.cpp @@ -25,12 +25,30 @@ namespace Opm { + // ---------------- Methods for class DGBasisInterface ---------------- + /// Virtual destructor. DGBasisInterface::~DGBasisInterface() { } + /// Evaluate function f = sum_i c_i b_i at the point x. + /// Note that this function is not virtual, but implemented in + /// terms of the virtual functions of the class. + /// \param[in] cell Cell index + /// \param[in] coefficients Coefficients {c_i} for a single cell. + /// \param[in] x Point at which to compute f(x). + double DGBasisInterface::evalFunc(const int cell, + const double* coefficients, + const double* x) const + { + bvals_.resize(numBasisFunc()); + eval(cell, x, &bvals_[0]); + return std::inner_product(bvals_.begin(), bvals_.end(), coefficients, 0.0); + } + + // ---------------- Methods for class DGBasisBoundedTotalDegree ---------------- @@ -154,6 +172,12 @@ namespace Opm } } + /// Compute the average of the function f = sum_i c_i b_i. + /// \param[in] coefficients Coefficients {c_i} for a single cell. + double DGBasisBoundedTotalDegree::functionAverage(const double* coefficients) const + { + return coefficients[0]; + } @@ -300,11 +324,18 @@ namespace Opm double* coefficients) const { const int nb = numBasisFunc(); - const double average = std::accumulate(coefficients, coefficients + nb, 0.0)/double(nb); + const double aver = functionAverage(coefficients); for (int ix = 0; ix < nb; ++ix) { - coefficients[ix] = factor*(coefficients[ix] - average) + average; + coefficients[ix] = factor*(coefficients[ix] - aver) + aver; } } + /// Compute the average of the function f = sum_i c_i b_i. + /// \param[in] coefficients Coefficients {c_i} for a single cell. + double DGBasisMultilin::functionAverage(const double* coefficients) const + { + const int nb = numBasisFunc(); + return std::accumulate(coefficients, coefficients + nb, 0.0)/double(nb); + } } // namespace Opm diff --git a/opm/core/transport/reorder/DGBasis.hpp b/opm/core/transport/reorder/DGBasis.hpp index a74108148..f1fc72341 100644 --- a/opm/core/transport/reorder/DGBasis.hpp +++ b/opm/core/transport/reorder/DGBasis.hpp @@ -20,6 +20,8 @@ #ifndef OPM_DGBASIS_HEADER_INCLUDED #define OPM_DGBASIS_HEADER_INCLUDED +#include + struct UnstructuredGrid; namespace Opm @@ -75,6 +77,23 @@ namespace Opm /// \param[out] coefficients Coefficients {c_i} for a single cell. virtual void multiplyGradient(const double factor, double* coefficients) const = 0; + + /// Evaluate function f = sum_i c_i b_i at the point x. + /// Note that this function is not virtual, but implemented in + /// terms of the virtual functions of the class. + /// \param[in] cell Cell index + /// \param[in] coefficients Coefficients {c_i} for a single cell. + /// \param[in] x Point at which to compute f(x). + double evalFunc(const int cell, + const double* coefficients, + const double* x) const; + + /// Compute the average of the function f = sum_i c_i b_i. + /// \param[in] coefficients Coefficients {c_i} for a single cell. + virtual double functionAverage(const double* coefficients) const = 0; + + private: + mutable std::vector bvals_; // For evalFunc(). }; @@ -144,6 +163,10 @@ namespace Opm virtual void multiplyGradient(const double factor, double* coefficients) const; + /// Compute the average of the function f = sum_i c_i b_i. + /// \param[in] coefficients Coefficients {c_i} for a single cell. + virtual double functionAverage(const double* coefficients) const; + private: const UnstructuredGrid& grid_; const int degree_; @@ -217,6 +240,10 @@ namespace Opm virtual void multiplyGradient(const double factor, double* coefficients) const; + /// Compute the average of the function f = sum_i c_i b_i. + /// \param[in] coefficients Coefficients {c_i} for a single cell. + virtual double functionAverage(const double* coefficients) const; + private: const UnstructuredGrid& grid_; const int degree_; diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 132041a4f..54a6e5858 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -411,10 +411,10 @@ namespace Opm { switch (limiter_method_) { case MinUpwindFace: - applyMinUpwindFaceLimiter(cell, tof); + applyMinUpwindLimiter(cell, true, tof); break; case MinUpwindAverage: - applyMinUpwindAverageLimiter(cell, tof); + applyMinUpwindLimiter(cell, false, tof); break; default: THROW("Limiter type not implemented: " << limiter_method_); @@ -424,48 +424,32 @@ namespace Opm - void TransportModelTracerTofDiscGal::applyMinUpwindFaceLimiter(const int cell, double* tof) + void TransportModelTracerTofDiscGal::applyMinUpwindLimiter(const int cell, const bool face_min, double* tof) { if (basis_func_->degree() != 1) { THROW("This limiter only makes sense for our DG1 implementation."); } // Limiter principles: - // 1. Let M be the minimum TOF value on the upstream faces, - // evaluated in the upstream cells. Then the value at all - // points in this cell shall be at least M. - // Upstream faces whose flux does not exceed the relative - // flux threshold are not considered for this minimum. + // 1. Let M be either: + // - the minimum TOF value of all upstream faces, + // evaluated in the upstream cells + // (chosen if face_min is true). + // or: + // - the minimum average TOF value of all upstream cells + // (chosen if face_min is false). + // Then the value at all points in this cell shall be at + // least M. Upstream faces whose flux does not exceed the + // relative flux threshold are not considered for this + // minimum. // 2. The TOF shall not be below zero in any point. - // Find total upstream/downstream fluxes. - double upstream_flux = 0.0; - double downstream_flux = 0.0; - for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { - const int face = grid_.cell_faces[hface]; - double flux = 0.0; - if (cell == grid_.face_cells[2*face]) { - flux = darcyflux_[face]; - } else { - flux = -darcyflux_[face]; - } - if (flux < 0.0) { - upstream_flux += flux; - } else { - downstream_flux += flux; - } - } - // In the presence of sources, significant fluxes may be missing from the computed fluxes, - // setting the total flux to the (positive) maximum avoids this: since source is either - // inflow or outflow, not both, either upstream_flux or downstream_flux must be correct. - const double total_flux = std::max(-upstream_flux, downstream_flux); - - // Find minimum tof on upstream faces and for this cell. - const int dim = grid_.dimensions; + // Find minimum tof on upstream faces/cells and for this cell. const int num_basis = basis_func_->numBasisFunc(); double min_upstream_tof = 1e100; double min_here_tof = 1e100; int num_upstream_faces = 0; + const double total_flux = totalFlux(cell); for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { const int face = grid_.cell_faces[hface]; double flux = 0.0; @@ -478,30 +462,22 @@ namespace Opm upstream_cell = grid_.face_cells[2*face]; } const bool upstream = (flux < -total_flux*limiter_relative_flux_threshold_); + const bool interior = (upstream_cell >= 0); + + // Find minimum tof in this cell and upstream. + // The meaning of minimum upstream tof depends on method. + min_here_tof = std::min(min_here_tof, minCornerVal(cell, face)); if (upstream) { ++num_upstream_faces; - } - bool interior = (upstream_cell >= 0); - - // Evaluate the solution in all corners. - for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) { - const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode]; - basis_func_->eval(cell, nc, &basis_[0]); - const double tof_here = std::inner_product(basis_.begin(), basis_.end(), - tof_coeff_ + num_basis*cell, 0.0); - min_here_tof = std::min(min_here_tof, tof_here); - if (upstream) { - if (interior) { - basis_func_->eval(upstream_cell, nc, &basis_nb_[0]); - const double tof_upstream - = std::inner_product(basis_nb_.begin(), basis_nb_.end(), - tof_coeff_ + num_basis*upstream_cell, 0.0); - min_upstream_tof = std::min(min_upstream_tof, tof_upstream); + double upstream_tof = 0.0; + if (interior) { + if (face_min) { + upstream_tof = minCornerVal(upstream_cell, face); } else { - // Allow tof down to 0 on inflow boundaries. - min_upstream_tof = std::min(min_upstream_tof, 0.0); + upstream_tof = basis_func_->functionAverage(tof_coeff_ + num_basis*upstream_cell); } } + min_upstream_tof = std::min(min_upstream_tof, upstream_tof); } } @@ -513,7 +489,7 @@ namespace Opm if (min_upstream_tof < 0.0) { min_upstream_tof = 0.0; } - const double tof_c = tof_coeff_[num_basis*cell]; + const double tof_c = basis_func_->functionAverage(tof_coeff_ + num_basis*cell); double limiter = (tof_c - min_upstream_tof)/(tof_c - min_here_tof); if (tof_c < min_upstream_tof) { // Handle by setting a flat solution. @@ -535,110 +511,6 @@ namespace Opm - void TransportModelTracerTofDiscGal::applyMinUpwindAverageLimiter(const int cell, double* tof) - { - if (basis_func_->degree() != 1) { - THROW("This limiter only makes sense for our DG1 implementation."); - } - - // Limiter principles: - // 1. Let M be the average TOF value of the upstream cells. - /// Then the value at all points in this cell shall be at least M. - // Upstream faces whose flux does not exceed the relative - // flux threshold are not considered for this minimum. - // 2. The TOF shall not be below zero in any point. - - // Find total upstream/downstream fluxes. - double upstream_flux = 0.0; - double downstream_flux = 0.0; - for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { - const int face = grid_.cell_faces[hface]; - double flux = 0.0; - if (cell == grid_.face_cells[2*face]) { - flux = darcyflux_[face]; - } else { - flux = -darcyflux_[face]; - } - if (flux < 0.0) { - upstream_flux += flux; - } else { - downstream_flux += flux; - } - } - // In the presence of sources, significant fluxes may be missing from the computed fluxes, - // setting the total flux to the (positive) maximum avoids this: since source is either - // inflow or outflow, not both, either upstream_flux or downstream_flux must be correct. - const double total_flux = std::max(-upstream_flux, downstream_flux); - - // Find minimum tof on upstream faces and for this cell. - const int dim = grid_.dimensions; - const int num_basis = basis_func_->numBasisFunc(); - double min_upstream_tof = 1e100; - double min_here_tof = 1e100; - int num_upstream_faces = 0; - for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { - const int face = grid_.cell_faces[hface]; - double flux = 0.0; - int upstream_cell = -1; - if (cell == grid_.face_cells[2*face]) { - flux = darcyflux_[face]; - upstream_cell = grid_.face_cells[2*face+1]; - } else { - flux = -darcyflux_[face]; - upstream_cell = grid_.face_cells[2*face]; - } - const bool upstream = (flux < -total_flux*limiter_relative_flux_threshold_); - if (upstream) { - ++num_upstream_faces; - } - bool interior = (upstream_cell >= 0); - - // Evaluate the solution in all corners. - for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) { - const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode]; - basis_func_->eval(cell, nc, &basis_[0]); - const double tof_here = std::inner_product(basis_.begin(), basis_.end(), - tof_coeff_ + num_basis*cell, 0.0); - min_here_tof = std::min(min_here_tof, tof_here); - if (upstream) { - if (interior) { - min_upstream_tof = std::min(min_upstream_tof, tof_coeff_[num_basis*upstream_cell]); - } else { - // Allow tof down to 0 on inflow boundaries. - min_upstream_tof = std::min(min_upstream_tof, 0.0); - } - } - } - } - - // Compute slope multiplier (limiter). - if (num_upstream_faces == 0) { - min_upstream_tof = 0.0; - min_here_tof = 0.0; - } - if (min_upstream_tof < 0.0) { - min_upstream_tof = 0.0; - } - const double tof_c = tof_coeff_[num_basis*cell]; - double limiter = (tof_c - min_upstream_tof)/(tof_c - min_here_tof); - if (tof_c < min_upstream_tof) { - // Handle by setting a flat solution. - std::cout << "Trouble in cell " << cell << std::endl; - limiter = 0.0; - basis_func_->addConstant(min_upstream_tof - tof_c, tof + num_basis*cell); - } - ASSERT(limiter >= 0.0); - - // Actually do the limiting (if applicable). - if (limiter < 1.0) { - // std::cout << "Applying limiter in cell " << cell << ", limiter = " << limiter << std::endl; - basis_func_->multiplyGradient(limiter, tof + num_basis*cell); - } else { - // std::cout << "Not applying limiter in cell " << cell << "!" << std::endl; - } - } - - void TransportModelTracerTofDiscGal::applyLimiterAsPostProcess() @@ -675,4 +547,51 @@ namespace Opm + double TransportModelTracerTofDiscGal::totalFlux(const int cell) const + { + // Find total upstream/downstream fluxes. + double upstream_flux = 0.0; + double downstream_flux = 0.0; + for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { + const int face = grid_.cell_faces[hface]; + double flux = 0.0; + if (cell == grid_.face_cells[2*face]) { + flux = darcyflux_[face]; + } else { + flux = -darcyflux_[face]; + } + if (flux < 0.0) { + upstream_flux += flux; + } else { + downstream_flux += flux; + } + } + // In the presence of sources, significant fluxes may be missing from the computed fluxes, + // setting the total flux to the (positive) maximum avoids this: since source is either + // inflow or outflow, not both, either upstream_flux or downstream_flux must be correct. + return std::max(-upstream_flux, downstream_flux); + } + + + + + double TransportModelTracerTofDiscGal::minCornerVal(const int cell, const int face) const + { + // Evaluate the solution in all corners. + const int dim = grid_.dimensions; + const int num_basis = basis_func_->numBasisFunc(); + double min_cornerval = 1e100; + for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) { + const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode]; + basis_func_->eval(cell, nc, &basis_[0]); + const double tof_corner = std::inner_product(basis_.begin(), basis_.end(), + tof_coeff_ + num_basis*cell, 0.0); + min_cornerval = std::min(min_cornerval, tof_corner); + } + return min_cornerval; + } + + + + } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index bcee06dfc..cf55ca8cb 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -119,8 +119,8 @@ namespace Opm std::vector orig_jac_; // single-cell jacobian (copy) // Below: storage for quantities needed by solveSingleCell(). std::vector coord_; - std::vector basis_; - std::vector basis_nb_; + mutable std::vector basis_; + mutable std::vector basis_nb_; std::vector grad_basis_; std::vector velocity_; @@ -130,10 +130,11 @@ namespace Opm // (will read data from tof_coeff_, it is ok to call // with tof_coeff as tof argument. void applyLimiter(const int cell, double* tof); - void applyMinUpwindFaceLimiter(const int cell, double* tof); - void applyMinUpwindAverageLimiter(const int cell, double* tof); + void applyMinUpwindLimiter(const int cell, const bool face_min, double* tof); void applyLimiterAsPostProcess(); void applyLimiterAsSimultaneousPostProcess(); + double totalFlux(const int cell) const; + double minCornerVal(const int cell, const int face) const; }; } // namespace Opm diff --git a/opm/core/transport/reorder/tarjan.c b/opm/core/transport/reorder/tarjan.c index 9150c2092..d3245610b 100644 --- a/opm/core/transport/reorder/tarjan.c +++ b/opm/core/transport/reorder/tarjan.c @@ -20,8 +20,8 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ -#include #include +#include #ifdef MATLAB_MEX_FILE #include "tarjan.h" @@ -30,7 +30,13 @@ SOFTWARE. #endif +static void +clear_vector(size_t n, int *v) +{ + size_t i; + for (i = 0; i < n; i++) { v[i] = 0; } +} static int min(int a, int b){ return a < b? a : b;} @@ -73,19 +79,21 @@ tarjan (int nv, const int *ia, const int *ja, int *vert, int *comp, int *stack = comp + nv; int *bottom = stack; int *cstack = vert + nv-1; + +#if !defined(NDEBUG) int *cbottom = cstack; +#endif + int t = 0; int pos = 0; int *time = work; - int *link = (int *) time + nv; - int *status = (int *) link + nv; /* dual usage... */ + int *link = time + nv; + int *status = link + nv; /* dual usage... */ - (void) cbottom; - - memset(work, 0, 3*nv * sizeof *work); - memset(vert, 0, nv * sizeof *vert ); - memset(comp, 0, (nv+1) * sizeof *comp ); + clear_vector(3 * ((size_t) nv), work); + clear_vector(1 * ((size_t) nv), vert); + clear_vector(1 + ((size_t) nv), comp); /* Init status all vertices */ for (i=0; i. - */ + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ #include diff --git a/opm/core/wells/WellCollection.hpp b/opm/core/wells/WellCollection.hpp index 1957576a6..4bda76d9b 100644 --- a/opm/core/wells/WellCollection.hpp +++ b/opm/core/wells/WellCollection.hpp @@ -1,22 +1,22 @@ /* -Copyright 2011 SINTEF ICT, Applied Mathematics. + Copyright 2011 SINTEF ICT, Applied Mathematics. + This file is part of the Open Porous Media project (OPM). -This file is part of The Open Reservoir Simulator Project (OpenRS). + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. -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. + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. -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 OPM. If not, see . +*/ -You should have received a copy of the GNU General Public License -along with OpenRS. If not, see . - */ #ifndef OPM_WELLCOLLECTION_HPP diff --git a/tests/test_dgbasis.cpp b/tests/test_dgbasis.cpp new file mode 100644 index 000000000..d5eb8c8c7 --- /dev/null +++ b/tests/test_dgbasis.cpp @@ -0,0 +1,160 @@ +/* + Copyright 2013 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif +#define NVERBOSE // to suppress our messages when throwing + +#define BOOST_TEST_MODULE DGBasisTest +#include + +#include +#include +#include +#include + +using namespace Opm; + + +namespace +{ + + + bool aequal(double a, double b) + { + const double eps = 1e-15; + return std::fabs(a - b) < eps; + } + +} // anonymous namespace + + +namespace cart2d +{ + + static void test() + { + // Set up 2d 1-cell cartesian case. + GridManager g(1, 1); + const UnstructuredGrid& grid = *g.c_grid(); + + // Test DGBasisBoundedTotalDegree, degree 0. + { + DGBasisBoundedTotalDegree b(grid, 0); + BOOST_CHECK_EQUAL(b.numBasisFunc(), 1); + std::vector bx(b.numBasisFunc(), 0.0); + b.eval(0, grid.cell_centroids, &bx[0]); + BOOST_CHECK(aequal(bx[0], 1.0)); + double x[2] = { 0.123, 0.456 }; + b.eval(0, x, &bx[0]); + BOOST_CHECK(aequal(bx[0], 1.0)); + std::vector c(b.numBasisFunc(), 0.0); + b.addConstant(0.789, &c[0]); + BOOST_CHECK(aequal(c[0], 0.789)); + b.multiplyGradient(1.234, &c[0]); + BOOST_CHECK(aequal(c[0], 0.789)); + } + // Test DGBasisBoundedTotalDegree, degree 1. + { + DGBasisBoundedTotalDegree b(grid, 1); + BOOST_CHECK_EQUAL(b.numBasisFunc(), 3); + std::vector bx(b.numBasisFunc(), 0.0); + b.eval(0, grid.cell_centroids, &bx[0]); + BOOST_CHECK(aequal(bx[0], 1.0)); + BOOST_CHECK(aequal(bx[1], 0.0)); + BOOST_CHECK(aequal(bx[2], 0.0)); + double x[2] = { 0.123, 0.456 }; + b.eval(0, x, &bx[0]); + BOOST_CHECK(aequal(bx[0], 1.0)); + BOOST_CHECK(aequal(bx[1], 0.123 - 0.5)); + BOOST_CHECK(aequal(bx[2], 0.456 - 0.5)); + std::vector c(b.numBasisFunc(), 0.0); + c[0] = 1.0; c[1] = 2.0; c[2] = 3.0; + b.addConstant(0.789, &c[0]); + BOOST_CHECK(aequal(c[0], 1.789)); + BOOST_CHECK(aequal(c[1], 2.0)); + BOOST_CHECK(aequal(c[2], 3.0)); + const double fx = c[0]*bx[0] + c[1]*bx[1] + c[2]*bx[2]; + b.multiplyGradient(1.234, &c[0]); + BOOST_CHECK(aequal(c[0], 1.789)); + BOOST_CHECK(aequal(c[1], 2.0*1.234)); + BOOST_CHECK(aequal(c[2], 3.0*1.234)); + const double fx2 = c[0]*bx[0] + c[1]*bx[1] + c[2]*bx[2]; + BOOST_CHECK(aequal(fx2 - c[0], 1.234*(fx - c[0]))); + } + // Test DGBasisMultilin, degree 0. + { + DGBasisMultilin b(grid, 0); + BOOST_CHECK_EQUAL(b.numBasisFunc(), 1); + std::vector bx(b.numBasisFunc(), 0.0); + b.eval(0, grid.cell_centroids, &bx[0]); + BOOST_CHECK(aequal(bx[0], 1.0)); + double x[2] = { 0.123, 0.456 }; + b.eval(0, x, &bx[0]); + BOOST_CHECK(aequal(bx[0], 1.0)); + std::vector c(b.numBasisFunc(), 0.0); + b.addConstant(0.789, &c[0]); + BOOST_CHECK(aequal(c[0], 0.789)); + b.multiplyGradient(1.234, &c[0]); + BOOST_CHECK(aequal(c[0], 0.789)); + } + // Test DGBasisMultilin, degree 1. + { + DGBasisMultilin b(grid, 1); + BOOST_CHECK_EQUAL(b.numBasisFunc(), 4); + std::vector bx(b.numBasisFunc(), 0.0); + b.eval(0, grid.cell_centroids, &bx[0]); + BOOST_CHECK(aequal(bx[0], 0.25)); + BOOST_CHECK(aequal(bx[1], 0.25)); + BOOST_CHECK(aequal(bx[2], 0.25)); + BOOST_CHECK(aequal(bx[3], 0.25)); + double x[2] = { 0.123, 0.456 }; + b.eval(0, x, &bx[0]); + const double xm[2] = { 1.0 - x[0], x[0] }; + const double ym[2] = { 1.0 - x[1], x[1] }; + BOOST_CHECK(aequal(bx[0], xm[0]*ym[0])); + BOOST_CHECK(aequal(bx[1], xm[0]*ym[1])); + BOOST_CHECK(aequal(bx[2], xm[1]*ym[0])); + BOOST_CHECK(aequal(bx[3], xm[1]*ym[1])); + std::vector c(b.numBasisFunc(), 0.0); + c[0] = -1.567; c[1] = 1.42; c[2] = 0.59; c[3] = 3.225; + std::vector corig = c; + b.addConstant(0.789, &c[0]); + BOOST_CHECK(aequal(c[0], corig[0] + 0.25*0.789)); + BOOST_CHECK(aequal(c[1], corig[1] + 0.25*0.789)); + BOOST_CHECK(aequal(c[2], corig[2] + 0.25*0.789)); + BOOST_CHECK(aequal(c[3], corig[3] + 0.25*0.789)); + const double fx = c[0]*bx[0] + c[1]*bx[1] + c[2]*bx[2] + c[3]*bx[3]; + const double fc = 0.25*(c[0] + c[1] + c[2] + c[3]); + b.multiplyGradient(1.234, &c[0]); + const double fx2 = c[0]*bx[0] + c[1]*bx[1] + c[2]*bx[2] + c[3]*bx[3]; + BOOST_CHECK(aequal(fx2 - fc, 1.234*(fx - fc))); + } + + } +} // namespace cart2d + + +BOOST_AUTO_TEST_CASE(test_dgbasis) +{ + cart2d::test(); +}