You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
51 lines
1.6 KiB
51 lines
1.6 KiB
#include <medusa/Medusa.hpp>
|
|
#include <medusa/bits/domains/PolyhedronShape.hpp> // This must be included separately.
|
|
#include <Eigen/Sparse>
|
|
|
|
// Example using PolyhedronShape with realistic 3D objects.
|
|
// http://e6.ijs.si/medusa/wiki/index.php/Realistic_3D_models
|
|
// Needs CGAL (https://www.cgal.org/).
|
|
|
|
using namespace mm; // NOLINT
|
|
|
|
int main() {
|
|
std::string filename = "../../test/testdata/bunny.off";
|
|
PolytopeShape<Vec3d> bunny = PolytopeShape<Vec3d>::fromOFF(filename);
|
|
|
|
double dx = 2.5;
|
|
DomainDiscretization<Vec3d> domain = bunny.discretizeBoundaryWithStep(dx);
|
|
GeneralFill<Vec3d> fill;
|
|
fill.seed(0);
|
|
fill(domain, dx);
|
|
|
|
int N = domain.size();
|
|
prn(N);
|
|
domain.findSupport(FindClosest(25));
|
|
RBFFD<Polyharmonic<double, 3>, Vec3d, ScaleToClosest> approx({}, 2);
|
|
auto storage = domain.template computeShapes<sh::lap|sh::d1>(approx);
|
|
|
|
Eigen::SparseMatrix<double, Eigen::RowMajor> M(N, N);
|
|
Eigen::VectorXd rhs(N); rhs.setZero();
|
|
|
|
auto op = storage.implicitOperators(M, rhs);
|
|
M.reserve(storage.supportSizes());
|
|
|
|
for (int i : domain.interior()) {
|
|
0.1*op.grad(i, -1) + 2.0*op.lap(i) = -1.0;
|
|
}
|
|
for (int i : domain.boundary()) {
|
|
op.value(i) = 0.0;
|
|
}
|
|
|
|
Eigen::BiCGSTAB<decltype(M), Eigen::IncompleteLUT<double>> solver;
|
|
solver.preconditioner().setDroptol(1e-4);
|
|
solver.preconditioner().setFillfactor(20);
|
|
solver.compute(M);
|
|
ScalarFieldd u = solver.solve(rhs);
|
|
|
|
HDF hdf("bunny_poisson.h5", HDF::DESTROY);
|
|
hdf.writeDouble2DArray("pos", domain.positions());
|
|
hdf.writeEigen("sol", u);
|
|
|
|
return 0;
|
|
}
|
|
|