DGtal 1.3.0
Loading...
Searching...
No Matches
dgtalCalculus-poisson.cpp
1
25#include <iostream>
26#include <DGtal/base/Common.h>
27#include <DGtal/helpers/StdDefs.h>
28#include <DGtal/helpers/Shortcuts.h>
29#include <DGtal/helpers/ShortcutsGeometry.h>
30#include <DGtal/shapes/SurfaceMesh.h>
31#include <DGtal/geometry/surfaces/DigitalSurfaceRegularization.h>
32#include <DGtal/dec/PolygonalCalculus.h>
33#include <DGtal/math/linalg/DirichletConditions.h>
34
35#include <polyscope/polyscope.h>
36#include <polyscope/surface_mesh.h>
37#include <polyscope/point_cloud.h>
38
39
40#include <Eigen/Dense>
41#include <Eigen/Sparse>
42
43using namespace DGtal;
44using namespace Z3i;
45
46// Using standard 3D digital space.
49// The following typedefs are useful
51typedef std::size_t Index;
52//Polyscope global
53polyscope::SurfaceMesh *psMesh;
54SurfMesh surfmesh;
55float scale = 0.1;
56
57void computeLaplace()
58{
62 PolyDEC calculus(surfmesh);
63 PolyDEC::SparseMatrix L = calculus.globalLaplaceBeltrami();
64 PolyDEC::Form g = calculus.form0();
65 DC::IntegerVector b = DC::IntegerVector::Zero( g.rows() );
66
67 //We set values on the boundary
68 auto boundaryEdges = surfmesh.computeManifoldBoundaryEdges();
69 std::cout<< "Number of boundary edges= "<<boundaryEdges.size()<<std::endl;
70
71 auto pihVertex=[&](const SurfMesh::Vertex &v){return cos(scale*(surfmesh.position(v)[0]))*(scale*surfmesh.position(v)[1]);};
72
73 for(auto &e: boundaryEdges)
74 {
75 auto adjVertices = surfmesh.edgeVertices(e);
76 g(adjVertices.first) = pihVertex(adjVertices.first);
77 g(adjVertices.second) = pihVertex(adjVertices.second);
78 b(adjVertices.first) = 1;
79 b(adjVertices.second) = 1;
80 }
81
82 // Solve Δu=0 with g as boundary conditions
83 PolyDEC::Solver solver;
84 PolyDEC::SparseMatrix L_dirichlet = DC::dirichletOperator( L, b );
85 solver.compute( L_dirichlet );
86 ASSERT(solver.info()==Eigen::Success);
87 PolyDEC::Form g_dirichlet = DC::dirichletVector( L, g, b, g );
88 PolyDEC::Form x_dirichlet = solver.solve( g_dirichlet );
89 PolyDEC::Form u = DC::dirichletSolution( x_dirichlet, b, g );
91
92 psMesh->addVertexScalarQuantity("g", g);
93 psMesh->addVertexScalarQuantity("u", u);
94}
95
96void myCallback()
97{
98 ImGui::SliderFloat("Phi scale", &scale, 0., 10.);
99 if(ImGui::Button("Compute Laplace problem"))
100 computeLaplace();
101}
102
103int main()
104{
106
107 auto h=.7 ; //gridstep
108 params( "polynomial", "0.1*y*y -0.1*x*x - 2.0*z" )( "gridstep", h );
109 auto implicit_shape = SH3::makeImplicitShape3D ( params );
110 auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
111 auto K = SH3::getKSpace( params );
112 auto binary_image = SH3::makeBinaryImage( digitized_shape, params );
113 auto surface = SH3::makeDigitalSurface( binary_image, K, params );
114 auto primalSurface = SH3::makePrimalSurfaceMesh(surface);
115
116 //Need to convert the faces
117 std::vector<std::vector<SH3::SurfaceMesh::Vertex>> faces;
118 std::vector<RealPoint> positions = primalSurface->positions();
119 for(auto face= 0 ; face < primalSurface->nbFaces(); ++face)
120 faces.push_back(primalSurface->incidentVertices( face ));
121
122 surfmesh = SurfMesh(positions.begin(),
123 positions.end(),
124 faces.begin(),
125 faces.end());
126 psMesh = polyscope::registerSurfaceMesh("digital surface", positions, faces);
127
128 // Initialize polyscope
129 polyscope::init();
130
131 // Set the callback function
132 polyscope::state::userCallback = myCallback;
133 polyscope::show();
134 return EXIT_SUCCESS;
135}
Aim: A helper class to solve a system with Dirichlet boundary conditions.
Implements differential operators on polygonal surfaces from .
SparseMatrix globalLaplaceBeltrami(const double lambda=1.0) const
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
static Parameters parametersGeometryEstimation()
static Parameters defaultParameters()
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:105
static KSpace getKSpace(const Point &low, const Point &up, Parameters params=parametersKSpace())
Definition: Shortcuts.h:332
static CountedPtr< DigitizedImplicitShape3D > makeDigitizedImplicitShape3D(CountedPtr< ImplicitShape3D > shape, Parameters params=parametersDigitizedImplicitShape3D())
Definition: Shortcuts.h:523
static CountedPtr< DigitalSurface > makeDigitalSurface(CountedPtr< TPointPredicate > bimage, const KSpace &K, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1209
static Parameters defaultParameters()
Definition: Shortcuts.h:203
static CountedPtr< SurfaceMesh > makePrimalSurfaceMesh(Cell2Index &c2i, CountedPtr< ::DGtal::DigitalSurface< TContainer > > aSurface)
Definition: Shortcuts.h:2372
static CountedPtr< BinaryImage > makeBinaryImage(Domain shapeDomain)
Definition: Shortcuts.h:561
static CountedPtr< ImplicitShape3D > makeImplicitShape3D(const Parameters &params=parametersImplicitShape3D())
Definition: Shortcuts.h:282
SMesh::Index Index
DGtal is the top-level namespace which contains all DGtal functions and types.
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:92
const VertexPair & edgeVertices(Edge e) const
Definition: SurfaceMesh.h:329
Edges computeManifoldBoundaryEdges() const
RealPoint & position(Vertex v)
Definition: SurfaceMesh.h:637
int main(int argc, char **argv)
KSpace K