diff --git a/GridCplusplus.hpp b/GridCplusplus.hpp index 04d2a859..601efaf3 100644 --- a/GridCplusplus.hpp +++ b/GridCplusplus.hpp @@ -137,7 +137,7 @@ private: int num_cells = grid.numCells(); int num_nodes = grid.numVertices(); int num_faces = grid.numFaces(); - int dim = Grid::dimensions; + int dim = Grid::dimension; node_coordinates_.resize(dim*num_nodes); for (int n = 0; n < num_nodes; ++n) { for (int dd = 0; dd < dim; ++dd) { @@ -153,7 +153,7 @@ private: face_areas_[f] = grid.faceArea(f); for (int dd = 0; dd < dim; ++dd) { face_centroids_[dim*f + dd] = grid.faceCentroid(f)[dd]; - face_normals_[dim*f + dd] = grid.faceNormals(f)[dd]; + face_normals_[dim*f + dd] = grid.faceNormal(f)[dd]; } } diff --git a/LibMimetic.hpp b/LibMimetic.hpp index 482d2e07..64296a3d 100644 --- a/LibMimetic.hpp +++ b/LibMimetic.hpp @@ -34,6 +34,7 @@ #define SINTEF_LIBMIMETIC_HEADER #include "ifsh.h" +#include "mimetic.h" #include "GridCplusplus.hpp" #include @@ -42,6 +43,7 @@ /// Encapsulates the ifsh (= incompressible flow solver hybrid) solver modules. class Ifsh { +public: /// @brief /// Default constructor, does nothing. Ifsh() @@ -61,7 +63,7 @@ class Ifsh /// @tparam Grid This must conform to the SimpleGrid concept. /// @param grid The grid object. template - void init(const Grid& grid) + void init(const Grid& grid, const double* perm) { // Build C grid structure. grid_.init(grid); @@ -74,25 +76,59 @@ class Ifsh if (!data_) { throw std::runtime_error("Failed to initialize ifsh solver."); } + + // Compute inner products. + int num_cells = grid.numCells(); + ncf_.resize(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + ncf_[cell] = grid.numCellFaces(cell); + } + // Zero gravity (for now) + int ngconn = grid_.c_grid()->cell_facepos[num_cells]; + int ngconn2 = data_->sum_ngconn2; + Binv_.resize(ngconn2); + gpress_.clear(); + gpress_.resize(ngconn, 0.0); // Zero gravity for now! + + grid_t* g = grid_.c_grid(); + mim_ip_simple_all(g->number_of_cells, g->dimensions, + data_->max_ngconn, + &ncf_[0], g->cell_facepos, g->cell_faces, + g->face_cells, g->face_centroids, + g->face_normals, g->face_areas, + g->cell_centroids, g->cell_volumes, + const_cast(perm), &Binv_[0]); } /// @brief /// Assemble the sparse system. - void assemble(const std::vector& sources) + void assemble(const std::vector& sources, + const std::vector& total_mobilities, + const std::vector& omegas) { - flowbc_t* bc = 0; // TODO - double* src = 0; // &sources[0]; // TODO - double* Binv = 0; // TODO - double* gpress = 0; // TODO + // Noflow conditions for now. + int num_faces = grid_.c_grid()->number_of_faces; + std::vector bc_types(num_faces, FLUX); + std::vector bc_vals(num_faces, 0); + flowbc_t bc = { &bc_types[0], &bc_vals[0] }; + + // Source terms from user. + double* src = const_cast(&sources[0]); // Ugly? Yes. Safe? I think so. + + // Inner products are precomputed. + double* Binv = &Binv_[0]; + + // Gravity contribs are precomputed. + double* gpress = &gpress_[0]; // All well related things are zero. well_control_t* wctrl = 0; double* WI = 0; double* wdp = 0; - double* totmob = 0; // TODO - double* omega = 0; // TODO - ifsh_assemble(bc, src, Binv, gpress, wctrl, WI, wdp, totmob, omega, data_); + double* totmob = const_cast(&total_mobilities[0]); + double* omega = const_cast(&omega[0]); + ifsh_assemble(&bc, src, Binv, gpress, wctrl, WI, wdp, totmob, omega, data_); } private: @@ -104,21 +140,19 @@ private: ifsh_data* data_; // Grid. GridCplusplus grid_; + // Number of faces per cell. + std::vector ncf_; + // B^{-1} storage. + std::vector Binv_; + // Face pressures ??? + std::vector gpress_; + // Total mobilities. + std::vector totmob_; + // Gravity coefficients (\omega = sum_{i = 1}^{num phases}f_i \rho_i[TODO: check this]). + std::vector omega_; }; #if 0 -void -ifsh_assemble(flowbc_t *bc, - double *src, - double *Binv, - double *gpress, - well_control_t *wctrl, - double *WI, - double *wdp, - double *totmob, /* \sum_i \lambda_i */ - double *omega, /* \sum_i \rho_i f_i */ - struct ifsh_data *h); - void ifsh_press_flux(grid_t *G, struct ifsh_data *h, double *src, double *cpress, double *fflux);