DGtal 1.3.0
Searching...
No Matches
GeodesicsInHeat.h
1
17#pragma once
18
31#if defined(GeodesicsInHeat_RECURSES)
32#error Recursive header files inclusion detected in GeodesicsInHeat.h
33#else // defined(GeodesicsInHeat_RECURSES)
35#define GeodesicsInHeat_RECURSES
36
37#if !defined GeodesicsInHeat_h
39#define GeodesicsInHeat_h
40
42// Inclusions
43#include <iostream>
44#include "DGtal/base/Common.h"
45#include "DGtal/base/ConstAlias.h"
46#include "DGtal/math/linalg/DirichletConditions.h"
48
49namespace DGtal
50{
52 // template class GeodesicsInHeat
61 template <typename TPolygonalCalculus>
63 {
64 // ----------------------- Standard services ------------------------------
65 public:
66
67 typedef TPolygonalCalculus PolygonalCalculus;
76
80 GeodesicsInHeat() = delete;
81
85 {
86 myIsInit=false;
87 }
88
92 ~GeodesicsInHeat() = default;
93
98 GeodesicsInHeat ( const GeodesicsInHeat & other ) = delete;
99
104 GeodesicsInHeat ( GeodesicsInHeat && other ) = delete;
105
111 GeodesicsInHeat & operator= ( const GeodesicsInHeat & other ) = delete;
112
119
120
121
122 // ----------------------- Interface --------------------------------------
123
135 void init( double dt, double lambda = 1.0,
136 bool boundary_with_mixed_solution = false )
137 {
138 myIsInit = true;
139 myLambda = lambda;
140
141 SparseMatrix laplacian = myCalculus->globalLaplaceBeltrami( lambda );
142 SparseMatrix mass = myCalculus->globalLumpedMassMatrix();
143 myHeatOpe = mass - dt*laplacian;
144
145 //Prefactorizing
146 myPoissonSolver.compute( laplacian );
147 myHeatSolver.compute ( myHeatOpe );
148
149 //empty source
150 mySource = Vector::Zero(myCalculus->nbVertices());
151
152 // Manage boundaries
153 myManageBoundary = false;
154 if ( ! boundary_with_mixed_solution ) return;
155 myBoundary = IntegerVector::Zero(myCalculus->nbVertices());
156 const auto surfmesh = myCalculus->getSurfaceMeshPtr();
157 const auto edges = surfmesh->computeManifoldBoundaryEdges();
158 for ( auto e : edges )
159 {
160 const auto vtcs = surfmesh->edgeVertices( e );
161 myBoundary[ vtcs.first ] = 1;
162 myBoundary[ vtcs.second ] = 1;
163 }
164 myManageBoundary = ! edges.empty();
165 if ( ! myManageBoundary ) return;
166 // Prepare solver for a problem with Dirichlet conditions.
168 // Prefactoring
169 myHeatDirichletSolver.compute( heatOpe_d );
170 }
171
175 void addSource(const Vertex aV)
176 {
177 ASSERT_MSG(aV < myCalculus->nbVertices(), "Vertex is not in the surface mesh vertex range");
179 mySource( aV ) = 1.0;
180 }
181
185 {
186 mySource = Vector::Zero(myCalculus->nbVertices());
187 }
188
193 {
194 FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
195 return mySource;
196 }
197
198
202 {
203 FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
204 //Heat diffusion
205 Vector heatDiffusion = myHeatSolver.solve(mySource);
206 // Take care of boundaries
207 if ( myManageBoundary )
208 {
209 Vector bValues = Vector::Zero( myCalculus->nbVertices() );
211 myBoundary, bValues );
212 Vector bSol = myHeatDirichletSolver.solve( bSources );
213 Vector heatDiffusionDirichlet
214 = Conditions::dirichletSolution( bSol, myBoundary, bValues );
215 heatDiffusion = 0.5 * ( heatDiffusion + heatDiffusionDirichlet );
216 }
217
218 Vector divergence = Vector::Zero(myCalculus->nbVertices());
219 auto cpt=0;
220 auto surfmesh = myCalculus->getSurfaceMeshPtr();
221
222 // Heat, normalization and divergence per face
223 for(auto f=0; f< myCalculus->nbFaces(); ++f)
224 {
225 Vector faceHeat( myCalculus->degree(f));
226 cpt=0;
227 auto vertices = surfmesh->incidentVertices(f);
228 for(auto v: vertices)
229 {
230 faceHeat(cpt) = heatDiffusion( v );
231 ++cpt;
232 }
233 // ∇heat / ∣∣∇heat∣∣
234 Vector grad = -myCalculus->gradient(f) * faceHeat;
236
237 // div
238 DenseMatrix oneForm = myCalculus->flat(f)*grad;
239 Vector divergenceFace = myCalculus->divergence( f, myLambda ) * oneForm;
240 cpt=0;
241 for(auto v: vertices)
242 {
243 divergence(v) += divergenceFace(cpt);
244 ++cpt;
245 }
246 }
247
248 // Last Poisson solve
249 Vector distVec = myPoissonSolver.solve(divergence);
250
251 //Source val
252 auto sourceval = distVec(myLastSourceIndex);
253
254 //shifting the distances to get 0 at sources
255 return distVec - sourceval*Vector::Ones(myCalculus->nbVertices());
256 }
257
258
260 bool isValid() const
261 {
262 return myIsInit && myCalculus->isValid();
263 }
264
265 // ----------------------- Private --------------------------------------
266
267 private:
268
271
274
277
280
283
286
289
291 double myLambda;
292
296
299
302
303
304 }; // end of class GeodesicsInHeat
305} // namespace DGtal
306
307// //
309
310#endif // !defined GeodesicsInHeat_h
311
312#undef GeodesicsInHeat_RECURSES
313#endif // else defined(GeodesicsInHeat_RECURSES)
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Definition: ConstAlias.h:187
Aim: A helper class to solve a system with Dirichlet boundary conditions.
LinearAlgebraBackend::IntegerVector IntegerVector
static DenseVector dirichletVector(const SparseMatrix &A, const DenseVector &b, const IntegerVector &p, const DenseVector &u)
static SparseMatrix dirichletOperator(const SparseMatrix &A, const IntegerVector &p)
static DenseVector dirichletSolution(const DenseVector &xd, const IntegerVector &p, const DenseVector &u)
This class implements on polygonal surfaces (using Discrete differential calculus on polygonal surfa...
GeodesicsInHeat(GeodesicsInHeat &&other)=delete
DirichletConditions< LinAlgBackend > Conditions
PolygonalCalculus::LinAlg LinAlgBackend
double myLambda
Lambda parameter.
PolygonalCalculus::SparseMatrix SparseMatrix
Vertex myLastSourceIndex
Vertex index to the last source point (to shift the distances)
const PolygonalCalculus * myCalculus
The underlying PolygonalCalculus instance.
GeodesicsInHeat(const GeodesicsInHeat &other)=delete
void init(double dt, double lambda=1.0, bool boundary_with_mixed_solution=false)
Solver myHeatSolver
Heat solver.
bool myIsInit
Validitate flag.
SparseMatrix myHeatOpe
The operator for heat diffusion.
PolygonalCalculus::Vertex Vertex
void addSource(const Vertex aV)
IntegerVector myBoundary
The boundary characteristic vector.
Solver myPoissonSolver
Poisson solver.
Solver myHeatDirichletSolver
Heat solver with Dirichlet boundary conditions.
GeodesicsInHeat & operator=(const GeodesicsInHeat &other)=delete
PolygonalCalculus::Vector Vector
TPolygonalCalculus PolygonalCalculus
PolygonalCalculus::Solver Solver
GeodesicsInHeat(ConstAlias< PolygonalCalculus > calculus)
PolygonalCalculus::DenseMatrix DenseMatrix
Vector mySource
Source vector.
Conditions::IntegerVector IntegerVector
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
MySurfaceMesh::Vertex Vertex
Vertex type.
LinAlg::SolverSimplicialLDLT Solver
Type of a sparse matrix solver.
LinAlg::DenseMatrix DenseMatrix
Type of dense matrix.
LinAlg::DenseVector Vector
Type of Vector.
DGtal is the top-level namespace which contains all DGtal functions and types.
void laplacian(Shape &shape, const Options &options, std::function< double(const RealPoint3D &)> input_function, std::function< double(const RealPoint3D &)> target_function, int argc, char **argv)
Aim: Provide linear algebra backend using Eigen dense and sparse matrix as well as dense vector....
Definition: EigenSupport.h:97