DGtal  1.4.beta
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 
43 using namespace DGtal;
44 using namespace Z3i;
45 
46 // Using standard 3D digital space.
49 // The following typedefs are useful
51 typedef std::size_t Index;
52 //Polyscope global
53 polyscope::SurfaceMesh *psMesh;
54 SurfMesh surfmesh;
55 float scale = 0.1;
56 
57 void computeLaplace()
58 {
62  PolyDEC calculus(surfmesh);
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 
96 void myCallback()
97 {
98  ImGui::SliderFloat("Phi scale", &scale, 0., 10.);
99  if(ImGui::Button("Compute Laplace problem"))
100  computeLaplace();
101 }
102 
103 int main()
104 {
105  auto params = SH3::defaultParameters() | SHG3::defaultParameters() | SHG3::parametersGeometryEstimation();
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 }
DGtal::SurfaceMesh::Vertex
Index Vertex
Definition: SurfaceMesh.h:108
DGtal::SurfaceMesh::position
RealPoint & position(Vertex v)
Definition: SurfaceMesh.h:637
SH3
Shortcuts< KSpace > SH3
Definition: testArithmeticalDSSComputerOnSurfels.cpp:49
DGtal::PolygonalCalculus::globalLaplaceBeltrami
SparseMatrix globalLaplaceBeltrami(const double lambda=1.0) const
Definition: PolygonalCalculus.h:798
Index
SMesh::Index Index
Definition: fullConvexitySphereGeodesics.cpp:117
K
KSpace K
Definition: testCubicalComplex.cpp:62
DGtal::SurfaceMesh
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:91
DGtal::SurfaceMesh::computeManifoldBoundaryEdges
Edges computeManifoldBoundaryEdges() const
DGtal::Shortcuts::makeDigitalSurface
static CountedPtr< DigitalSurface > makeDigitalSurface(CountedPtr< TPointPredicate > bimage, const KSpace &K, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1209
DGtal::Shortcuts::makePrimalSurfaceMesh
static CountedPtr< SurfaceMesh > makePrimalSurfaceMesh(Cell2Index &c2i, CountedPtr< ::DGtal::DigitalSurface< TContainer > > aSurface)
Definition: Shortcuts.h:2372
DGtal::PolygonalCalculus
Implements differential operators on polygonal surfaces from .
Definition: PolygonalCalculus.h:70
DGtal::Shortcuts
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:104
DGtal::PolygonalCalculus::form0
Form form0() const
Definition: PolygonalCalculus.h:771
SparseMatrix
EigenLinearAlgebraBackend::SparseMatrix SparseMatrix
Definition: testHeatLaplace.cpp:50
SHG3
ShortcutsGeometry< Z3i::KSpace > SHG3
Definition: testIntegralInvariantShortcuts.cpp:47
DGtal::Shortcuts::defaultParameters
static Parameters defaultParameters()
Definition: Shortcuts.h:203
DGtal::Shortcuts::makeImplicitShape3D
static CountedPtr< ImplicitShape3D > makeImplicitShape3D(const Parameters &params=parametersImplicitShape3D())
Definition: Shortcuts.h:282
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::SurfaceMesh::edgeVertices
const VertexPair & edgeVertices(Edge e) const
Definition: SurfaceMesh.h:329
DGtal::ShortcutsGeometry
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: ShortcutsGeometry.h:76
DGtal::Shortcuts::makeDigitizedImplicitShape3D
static CountedPtr< DigitizedImplicitShape3D > makeDigitizedImplicitShape3D(CountedPtr< ImplicitShape3D > shape, Parameters params=parametersDigitizedImplicitShape3D())
Definition: Shortcuts.h:523
main
int main(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147
DGtal::Shortcuts::getKSpace
static KSpace getKSpace(const Point &low, const Point &up, Parameters params=parametersKSpace())
Definition: Shortcuts.h:332
DGtal::DirichletConditions
Aim: A helper class to solve a system with Dirichlet boundary conditions.
Definition: DirichletConditions.h:97
DGtal::Shortcuts::makeBinaryImage
static CountedPtr< BinaryImage > makeBinaryImage(Domain shapeDomain)
Definition: Shortcuts.h:561
init
K init(Point(0, 0, 0), Point(512, 512, 512), true)