From 171e4a51349354f8b0e03bdb6769bdb82c584c0c Mon Sep 17 00:00:00 2001 From: Wenzheng Chen Date: Thu, 27 Sep 2018 01:51:54 -0400 Subject: [PATCH 1/2] new --- src/fd_grad.cpp | 51 +++++++-- src/fd_interpolate.cpp | 102 +++++++++++++++--- src/fd_partial_derivative.cpp | 46 +++++++-- src/poisson_surface_reconstruction.cpp | 138 ++++++++++++++++--------- 4 files changed, 253 insertions(+), 84 deletions(-) diff --git a/src/fd_grad.cpp b/src/fd_grad.cpp index 9206181..5a24c35 100644 --- a/src/fd_grad.cpp +++ b/src/fd_grad.cpp @@ -1,13 +1,44 @@ #include "fd_grad.h" -void fd_grad( - const int nx, - const int ny, - const int nz, - const double h, - Eigen::SparseMatrix & G) -{ - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// +typedef Eigen::Triplet T; + +void intervalue2(std::vector &coef, int row, int aidx, int bidx, int dir) { + T tmp; + tmp[0] = row; + tmp[1] = bidx * 3 + dir; + tmp[2] = 1; + coef.push_back(tmp); + + tmp[0] = row; + tmp[1] = aidx * 3 + dir; + tmp[2] = -1; + coef.push_back(tmp); +} + +extern int gidx(double i, double j, double k, int nx, int ny); + +void fd_grad(const int nx, const int ny, const int nz, const double h, + Eigen::SparseMatrix & G) { + //////////////////////////////////////////////////////////////////////////// + // Add your code here + int nnum = nx * ny * nz; + G(nnum * 3, nnum * 3); + + // first order + std::vector coef; + for (int i = 0; i < nx - 1; i++) + for (int j = 0; j < ny - 1; j++) + for (int k = 0; k < nz - 1; k++) { + int st = gidx(i, j, k, nx, ny); + int edx = gidx(i + 1, j, k, nx, ny); + int edy = gidx(i, j + 1, k, nx, ny); + int edz = gidx(i, j, k + 1, nx, ny); + + int row = st * 3; + intervalue2(coef, row, st, edx, 0); + intervalue2(coef, row, st+1, edy, 1); + intervalue2(coef, row, st+2, edz, 2); + } + G.setFromTriplets(coef.begin(), coef.end()); + //////////////////////////////////////////////////////////////////////////// } diff --git a/src/fd_interpolate.cpp b/src/fd_interpolate.cpp index 2a32675..7ecb950 100644 --- a/src/fd_interpolate.cpp +++ b/src/fd_interpolate.cpp @@ -1,15 +1,93 @@ #include "fd_interpolate.h" -void fd_interpolate( - const int nx, - const int ny, - const int nz, - const double h, - const Eigen::RowVector3d & corner, - const Eigen::MatrixXd & P, - Eigen::SparseMatrix & W) -{ - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// +typedef Eigen::Triplet T; + +void intervalue(std::vector &coef, int aidx, int bidx, int value) { + T tmp; + for (int i = 0; i < 3; i++) { + tmp[0] = aidx * 3 + i; + tmp[1] = bidx * 3 + i; + tmp[2] = value; + coef.push_back(tmp); + } +} + +int gidx(double i, double j, double k, int nx, int ny) { + return i + nx * (j + k * ny); +} + +void fd_interpolate(const int nx, const int ny, const int nz, const double h, + const Eigen::RowVector3d & corner, const Eigen::MatrixXd & P, + Eigen::SparseMatrix & W) { + //////////////////////////////////////////////////////////////////////////// + // Add your code here + + int pnum = P.rows(); + int nnum = nx * ny * nz; + + Eigen::MatrixXd P2(pnum * 3, 1); + for (int i = 0; i < pnum; i++) { + for (int j = 0; j < 3; j++) { + P2(i * 3 + j, 0) = P(i, j); + } + } + + // size of w should be pnum * 3 x nnum * 3 + + std::vector coef; + for (int p = 0; p < pnum; p++) { + double px = P(p, 0) - corner[0]; + double py = P(p, 1) - corner[1]; + double pz = P(p, 2) - corner[2]; + + double xind = px / h; + double yind = py / h; + double zind = pz / h; + + // 8 equations + double xfloor = std::floor(px); + double yfloor = std::floor(py); + double zfloor = std::floor(pz); + double xceil = xfloor + 1; + double yceil = yfloor + 1; + double zceil = zfloor + 1; + + // left, left, left + double idx = gidx(xfloor, yfloor, zfloor, nx, ny); + double value = (xceil - px) * (yceil - py) * (zceil - pz); + intervalue(coef, p, idx, value); + + idx = gidx(xfloor, yceil, zfloor, nx, ny); + value = (xceil - px) * (py - yfloor) * (zceil - pz); + intervalue(coef, p, idx, value); + + idx = gidx(xceil, yceil, zfloor, nx, ny); + value = (px - xfloor) * (py - yfloor) * (zceil - pz); + intervalue(coef, p, idx, value); + + idx = gidx(xceil, yfloor, zfloor, nx, ny); + value = (px - xfloor) * (yceil - py) * (zceil - pz); + intervalue(coef, p, idx, value); + + /////////////////////////////////////////////////////// + idx = gidx(xfloor, yfloor, zceil, nx, ny); + value = (xceil - px) * (yceil - px) * (pz - zfloor); + intervalue(coef, p, idx, value); + + idx = gidx(xfloor, yceil, zceil, nx, ny); + value = (xceil - px) * (px - yfloor) * (pz - zfloor); + intervalue(coef, p, idx, value); + + idx = gidx(xceil, yceil, zceil, nx, ny); + value = (px - xfloor) * (px - yfloor) * (pz - zfloor); + intervalue(coef, p, idx, value); + + idx = gidx(xceil, yfloor, zceil, nx, ny); + value = (px - xfloor) * (yceil - py) * (pz - zfloor); + intervalue(coef, p, idx, value); + } + + W.setFromTriplets(coef.begin(), coef.end()); + + //////////////////////////////////////////////////////////////////////////// } diff --git a/src/fd_partial_derivative.cpp b/src/fd_partial_derivative.cpp index c3c4deb..25faf4c 100644 --- a/src/fd_partial_derivative.cpp +++ b/src/fd_partial_derivative.cpp @@ -1,14 +1,38 @@ #include "fd_partial_derivative.h" -void fd_partial_derivative( - const int nx, - const int ny, - const int nz, - const double h, - const int dir, - Eigen::SparseMatrix & D) -{ - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// +typedef Eigen::Triplet T; +extern int gidx(double i, double j, double k, int nx, int ny); +extern void intervalue2(std::vector &coef, int row, int aidx, int bidx, + int dir); + +void fd_partial_derivative(const int nx, const int ny, const int nz, + const double h, const int dir, Eigen::SparseMatrix & D) { + //////////////////////////////////////////////////////////////////////////// + // Add your code here + int nnum = nx * ny * nz; + D(nnum * 3, nnum * 3); + + // first order + std::vector coef; + for (int i = 0; i < nx - 1; i++) + for (int j = 0; j < ny - 1; j++) + for (int k = 0; k < nz - 1; k++) { + int st = gidx(i, j, k, nx, ny); + int edx = gidx(i + 1, j, k, nx, ny); + int edy = gidx(i, j + 1, k, nx, ny); + int edz = gidx(i, j, k + 1, nx, ny); + + int row = st * 1; + if (dir == 0) { + intervalue2(coef, row, st, edx, 0); + } else { + if (dir == 1) + intervalue2(coef, row, st, edy, 1); + else + intervalue2(coef, row, st, edz, 2); + } + } + + D.setFromTriplets(coef.begin(), coef.end()); + //////////////////////////////////////////////////////////////////////////// } diff --git a/src/poisson_surface_reconstruction.cpp b/src/poisson_surface_reconstruction.cpp index 4fe7c1e..4d61f31 100644 --- a/src/poisson_surface_reconstruction.cpp +++ b/src/poisson_surface_reconstruction.cpp @@ -2,55 +2,91 @@ #include #include -void poisson_surface_reconstruction( - const Eigen::MatrixXd & P, - const Eigen::MatrixXd & N, - Eigen::MatrixXd & V, - Eigen::MatrixXi & F) -{ - //////////////////////////////////////////////////////////////////////////// - // Construct FD grid, CONGRATULATIONS! You get this for free! - //////////////////////////////////////////////////////////////////////////// - // number of input points - const int n = P.rows(); - // Grid dimensions - int nx, ny, nz; - // Maximum extent (side length of bounding box) of points - double max_extent = - (P.colwise().maxCoeff()-P.colwise().minCoeff()).maxCoeff(); - // padding: number of cells beyond bounding box of input points - const double pad = 8; - // choose grid spacing (h) so that shortest side gets 30+2*pad samples - double h = max_extent/double(30+2*pad); - // Place bottom-left-front corner of grid at minimum of points minus padding - Eigen::RowVector3d corner = P.colwise().minCoeff().array()-pad*h; - // Grid dimensions should be at least 3 - nx = std::max((P.col(0).maxCoeff()-P.col(0).minCoeff()+(2.*pad)*h)/h,3.); - ny = std::max((P.col(1).maxCoeff()-P.col(1).minCoeff()+(2.*pad)*h)/h,3.); - nz = std::max((P.col(2).maxCoeff()-P.col(2).minCoeff()+(2.*pad)*h)/h,3.); - // Compute positions of grid nodes - Eigen::MatrixXd x(nx*ny*nz, 3); - for(int i = 0; i < nx; i++) - { - for(int j = 0; j < ny; j++) - { - for(int k = 0; k < nz; k++) - { - // Convert subscript to index - const auto ind = i + nx*(j + k * ny); - x.row(ind) = corner + h*Eigen::RowVector3d(i,j,k); - } - } - } - Eigen::VectorXd g = Eigen::VectorXd::Zero(nx*ny*nz); - - //////////////////////////////////////////////////////////////////////////// - // Add your code here - //////////////////////////////////////////////////////////////////////////// - - //////////////////////////////////////////////////////////////////////////// - // Run black box algorithm to compute mesh from implicit function: this - // function always extracts g=0, so "pre-shift" your g values by -sigma - //////////////////////////////////////////////////////////////////////////// - igl::copyleft::marching_cubes(g, x, nx, ny, nz, V, F); +#include "fd_interpolate.h" +#include "fd_partial_derivative.h" +#include "fd_grad.h" + +#include + +void poisson_surface_reconstruction(const Eigen::MatrixXd & P, + const Eigen::MatrixXd & N, Eigen::MatrixXd & V, Eigen::MatrixXi & F) { + //////////////////////////////////////////////////////////////////////////// + // Construct FD grid, CONGRATULATIONS! You get this for free! + //////////////////////////////////////////////////////////////////////////// + // number of input points + const int n = P.rows(); + // Grid dimensions + int nx, ny, nz; + // Maximum extent (side length of bounding box) of points + double max_extent = + (P.colwise().maxCoeff() - P.colwise().minCoeff()).maxCoeff(); + // padding: number of cells beyond bounding box of input points + const double pad = 8; + // choose grid spacing (h) so that shortest side gets 30+2*pad samples + double h = max_extent / double(30 + 2 * pad); + // Place bottom-left-front corner of grid at minimum of points minus padding + Eigen::RowVector3d corner = P.colwise().minCoeff().array() - pad * h; + // Grid dimensions should be at least 3 + nx = std::max( + (P.col(0).maxCoeff() - P.col(0).minCoeff() + (2. * pad) * h) / h, + 3.); + ny = std::max( + (P.col(1).maxCoeff() - P.col(1).minCoeff() + (2. * pad) * h) / h, + 3.); + nz = std::max( + (P.col(2).maxCoeff() - P.col(2).minCoeff() + (2. * pad) * h) / h, + 3.); + // Compute positions of grid nodes + Eigen::MatrixXd x(nx * ny * nz, 3); + for (int i = 0; i < nx; i++) { + for (int j = 0; j < ny; j++) { + for (int k = 0; k < nz; k++) { + // Convert subscript to index + const auto ind = i + nx * (j + k * ny); + x.row(ind) = corner + h * Eigen::RowVector3d(i, j, k); + } + } + } + Eigen::VectorXd g = Eigen::VectorXd::Zero(nx * ny * nz); + + //////////////////////////////////////////////////////////////////////////// + // Add your code here + //////////////////////////////////////////////////////////////////////////// + // first, we calculate spaese mtx M, which M_p3_xyz3 * X_xyz3_1 = P_p3_1 + + int pnum = n; + int gnum = nx * ny * nz; + Eigen::SparseMatrix W(pnum * 3, gnum * 3); + fd_interpolate(nx, ny, nz, h, corner, P, W); + + // next, we calculate sparse mtx G, which G_(x-1)(y-1)(z-1)_xyz * X_xyz_3 = D_(x-1)(y-1)(z-1)_3 + + Eigen::SparseMatrix G((nx - 1) * (ny - 1) * (nz - 1) * 3, + nx * ny * nz * 3); + fd_grad(nx, ny, nz, h, G); + + // next, we combine G1, G2, G3 with M + Eigen::SparseMatrix A = W * G; + + // last, we solve Ax = B + BiCGSTAB > solver; + solver.compute(A); + if(solver.info()!=Success) { + // decomposition failed + return; + } + g = solver.solve(b); + if(solver.info()!=Success) { + // solving failed + return; + } + + double sigma = (W * g).mean(); + g = g - sigma; + + //////////////////////////////////////////////////////////////////////////// + // Run black box algorithm to compute mesh from implicit function: this + // function always extracts g=0, so "pre-shift" your g values by -sigma + //////////////////////////////////////////////////////////////////////////// + igl::copyleft::marching_cubes(g, x, nx, ny, nz, V, F); } From eba63354ebbe2da539d0d840f18033d6c9b080a9 Mon Sep 17 00:00:00 2001 From: Wenzheng Chen Date: Thu, 27 Sep 2018 15:17:27 -0400 Subject: [PATCH 2/2] hw1 --- src/fd_grad.cpp | 35 +++++++++------------ src/fd_interpolate.cpp | 43 +++++++++++--------------- src/fd_partial_derivative.cpp | 36 ++++++++++----------- src/poisson_surface_reconstruction.cpp | 30 +++++++++--------- 4 files changed, 65 insertions(+), 79 deletions(-) diff --git a/src/fd_grad.cpp b/src/fd_grad.cpp index 5a24c35..be48f30 100644 --- a/src/fd_grad.cpp +++ b/src/fd_grad.cpp @@ -2,43 +2,38 @@ typedef Eigen::Triplet T; -void intervalue2(std::vector &coef, int row, int aidx, int bidx, int dir) { - T tmp; - tmp[0] = row; - tmp[1] = bidx * 3 + dir; - tmp[2] = 1; - coef.push_back(tmp); - - tmp[0] = row; - tmp[1] = aidx * 3 + dir; - tmp[2] = -1; +void intervalue2(std::vector &coef, int row, int aidx, int bidx) { + T tmp(row, aidx, -1); + T tmp2(row, bidx, 1); coef.push_back(tmp); + coef.push_back(tmp2); } -extern int gidx(double i, double j, double k, int nx, int ny); +extern int gidx(double i, double j, double k, int nx, int ny, int nz); void fd_grad(const int nx, const int ny, const int nz, const double h, Eigen::SparseMatrix & G) { //////////////////////////////////////////////////////////////////////////// // Add your code here - int nnum = nx * ny * nz; - G(nnum * 3, nnum * 3); // first order std::vector coef; for (int i = 0; i < nx - 1; i++) for (int j = 0; j < ny - 1; j++) for (int k = 0; k < nz - 1; k++) { - int st = gidx(i, j, k, nx, ny); - int edx = gidx(i + 1, j, k, nx, ny); - int edy = gidx(i, j + 1, k, nx, ny); - int edz = gidx(i, j, k + 1, nx, ny); + + int st = gidx(i, j, k, nx, ny, nz); + int edx = gidx(i + 1, j, k, nx, ny, nz); + int edy = gidx(i, j + 1, k, nx, ny, nz); + int edz = gidx(i, j, k + 1, nx, ny, nz); int row = st * 3; - intervalue2(coef, row, st, edx, 0); - intervalue2(coef, row, st+1, edy, 1); - intervalue2(coef, row, st+2, edz, 2); + intervalue2(coef, row, st, edx); + intervalue2(coef, row + 1, st, edy); + intervalue2(coef, row + 2, st, edz); } + G.setFromTriplets(coef.begin(), coef.end()); //////////////////////////////////////////////////////////////////////////// } + diff --git a/src/fd_interpolate.cpp b/src/fd_interpolate.cpp index 7ecb950..00e782f 100644 --- a/src/fd_interpolate.cpp +++ b/src/fd_interpolate.cpp @@ -1,19 +1,20 @@ #include "fd_interpolate.h" +#include typedef Eigen::Triplet T; -void intervalue(std::vector &coef, int aidx, int bidx, int value) { - T tmp; +void intervalue(std::vector &coef, double aidx, double bidx, double value) { + for (int i = 0; i < 3; i++) { - tmp[0] = aidx * 3 + i; - tmp[1] = bidx * 3 + i; - tmp[2] = value; + T tmp(aidx * 3 + i, bidx * 3 + i, value); coef.push_back(tmp); } } -int gidx(double i, double j, double k, int nx, int ny) { - return i + nx * (j + k * ny); +double gidx(double i, double j, double k, int nx, int ny, int nz) { + double idx = i + nx * (j + k * ny); + assert(idx < (nx - 1) * (ny - 1) * (nz - 1)); + return idx; } void fd_interpolate(const int nx, const int ny, const int nz, const double h, @@ -23,17 +24,8 @@ void fd_interpolate(const int nx, const int ny, const int nz, const double h, // Add your code here int pnum = P.rows(); - int nnum = nx * ny * nz; - - Eigen::MatrixXd P2(pnum * 3, 1); - for (int i = 0; i < pnum; i++) { - for (int j = 0; j < 3; j++) { - P2(i * 3 + j, 0) = P(i, j); - } - } - - // size of w should be pnum * 3 x nnum * 3 + // size of w should be pnum * 3 x (nx-1)(ny-1)(nz-1) * 3 std::vector coef; for (int p = 0; p < pnum; p++) { double px = P(p, 0) - corner[0]; @@ -53,36 +45,36 @@ void fd_interpolate(const int nx, const int ny, const int nz, const double h, double zceil = zfloor + 1; // left, left, left - double idx = gidx(xfloor, yfloor, zfloor, nx, ny); + double idx = gidx(xfloor, yfloor, zfloor, nx, ny, nz); double value = (xceil - px) * (yceil - py) * (zceil - pz); intervalue(coef, p, idx, value); - idx = gidx(xfloor, yceil, zfloor, nx, ny); + idx = gidx(xfloor, yceil, zfloor, nx, ny, nz); value = (xceil - px) * (py - yfloor) * (zceil - pz); intervalue(coef, p, idx, value); - idx = gidx(xceil, yceil, zfloor, nx, ny); + idx = gidx(xceil, yceil, zfloor, nx, ny, nz); value = (px - xfloor) * (py - yfloor) * (zceil - pz); intervalue(coef, p, idx, value); - idx = gidx(xceil, yfloor, zfloor, nx, ny); + idx = gidx(xceil, yfloor, zfloor, nx, ny, nz); value = (px - xfloor) * (yceil - py) * (zceil - pz); intervalue(coef, p, idx, value); /////////////////////////////////////////////////////// - idx = gidx(xfloor, yfloor, zceil, nx, ny); + idx = gidx(xfloor, yfloor, zceil, nx, ny, nz); value = (xceil - px) * (yceil - px) * (pz - zfloor); intervalue(coef, p, idx, value); - idx = gidx(xfloor, yceil, zceil, nx, ny); + idx = gidx(xfloor, yceil, zceil, nx, ny, nz); value = (xceil - px) * (px - yfloor) * (pz - zfloor); intervalue(coef, p, idx, value); - idx = gidx(xceil, yceil, zceil, nx, ny); + idx = gidx(xceil, yceil, zceil, nx, ny, nz); value = (px - xfloor) * (px - yfloor) * (pz - zfloor); intervalue(coef, p, idx, value); - idx = gidx(xceil, yfloor, zceil, nx, ny); + idx = gidx(xceil, yfloor, zceil, nx, ny, nz); value = (px - xfloor) * (yceil - py) * (pz - zfloor); intervalue(coef, p, idx, value); } @@ -91,3 +83,4 @@ void fd_interpolate(const int nx, const int ny, const int nz, const double h, //////////////////////////////////////////////////////////////////////////// } + diff --git a/src/fd_partial_derivative.cpp b/src/fd_partial_derivative.cpp index 25faf4c..7762dab 100644 --- a/src/fd_partial_derivative.cpp +++ b/src/fd_partial_derivative.cpp @@ -1,35 +1,35 @@ #include "fd_partial_derivative.h" typedef Eigen::Triplet T; -extern int gidx(double i, double j, double k, int nx, int ny); -extern void intervalue2(std::vector &coef, int row, int aidx, int bidx, - int dir); +extern int gidx(double i, double j, double k, int nx, int ny, int nz); +extern void intervalue2(std::vector &coef, int row, int aidx, int bidx); void fd_partial_derivative(const int nx, const int ny, const int nz, const double h, const int dir, Eigen::SparseMatrix & D) { //////////////////////////////////////////////////////////////////////////// // Add your code here - int nnum = nx * ny * nz; - D(nnum * 3, nnum * 3); - // first order std::vector coef; for (int i = 0; i < nx - 1; i++) for (int j = 0; j < ny - 1; j++) for (int k = 0; k < nz - 1; k++) { - int st = gidx(i, j, k, nx, ny); - int edx = gidx(i + 1, j, k, nx, ny); - int edy = gidx(i, j + 1, k, nx, ny); - int edz = gidx(i, j, k + 1, nx, ny); - int row = st * 1; - if (dir == 0) { - intervalue2(coef, row, st, edx, 0); - } else { - if (dir == 1) - intervalue2(coef, row, st, edy, 1); - else - intervalue2(coef, row, st, edz, 2); + int st = gidx(i, j, k, nx, ny, nz); + int edx = gidx(i + 1, j, k, nx, ny, nz); + int edy = gidx(i, j + 1, k, nx, ny, nz); + int edz = gidx(i, j, k + 1, nx, ny, nz); + + int row = st; + switch (dir) { + case 0: + intervalue2(coef, row, st, edx); + break; + case 1: + intervalue2(coef, row + 1, st, edy); + break; + case 2: + intervalue2(coef, row + 2, st, edz); + break; } } diff --git a/src/poisson_surface_reconstruction.cpp b/src/poisson_surface_reconstruction.cpp index 4d61f31..222d2b3 100644 --- a/src/poisson_surface_reconstruction.cpp +++ b/src/poisson_surface_reconstruction.cpp @@ -52,37 +52,35 @@ void poisson_surface_reconstruction(const Eigen::MatrixXd & P, //////////////////////////////////////////////////////////////////////////// // Add your code here //////////////////////////////////////////////////////////////////////////// - // first, we calculate spaese mtx M, which M_p3_xyz3 * X_xyz3_1 = P_p3_1 - + // first, we calculate spaese mtx M, which M_p3_xyz3 * X_(x-1)(y-1)(z-1)3_1 = P_p3_1 int pnum = n; - int gnum = nx * ny * nz; + int gnum = (nx - 1) * (ny - 1) * (nz - 1); Eigen::SparseMatrix W(pnum * 3, gnum * 3); fd_interpolate(nx, ny, nz, h, corner, P, W); - // next, we calculate sparse mtx G, which G_(x-1)(y-1)(z-1)_xyz * X_xyz_3 = D_(x-1)(y-1)(z-1)_3 - + // next, we calculate sparse mtx G, which G_(x-1)(y-1)(z-1)3_xyz * X_xyz_1 = D_(x-1)(y-1)(z-1)_1 Eigen::SparseMatrix G((nx - 1) * (ny - 1) * (nz - 1) * 3, - nx * ny * nz * 3); + nx * ny * nz * 1); fd_grad(nx, ny, nz, h, G); - // next, we combine G1, G2, G3 with M + // next, we combine G with M Eigen::SparseMatrix A = W * G; // last, we solve Ax = B - BiCGSTAB > solver; + Eigen::BiCGSTAB > solver; solver.compute(A); - if(solver.info()!=Success) { - // decomposition failed - return; + + Eigen::VectorXd b(n * 3); + for (int i = 0; i < pnum; i++) { + for (int j = 0; j < 3; j++) { + b(i * 3 + j, 0) = N(i, j); + } } g = solver.solve(b); - if(solver.info()!=Success) { - // solving failed - return; - } double sigma = (W * g).mean(); - g = g - sigma; + for (int i = 0; i < g.cols(); i++) + g(i) = g(i) - sigma; //////////////////////////////////////////////////////////////////////////// // Run black box algorithm to compute mesh from implicit function: this