DGtal 1.3.0
Loading...
Searching...
No Matches
VectorsInHeat.h
1
17#pragma once
18
31#if defined(VectorsInHeat_RECURSES)
32#error Recursive header files inclusion detected in VectorsInHeat.h
33#else // defined(VectorsInHeat_RECURSES)
35#define VectorsInHeat_RECURSES
36
37#if !defined VectorsInHeat_h
39#define VectorsInHeat_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 VectorInHeat
61template <typename TPolygonalCalculus>
63{
64 // ----------------------- Standard services ------------------------------
65public:
66
67 typedef TPolygonalCalculus PolygonalCalculus;
76
80 VectorsInHeat() = delete;
81
85 {
86 myIsInit=false;
87 }
88
92 ~VectorsInHeat() = default;
93
98 VectorsInHeat ( const VectorsInHeat & other ) = delete;
99
104 VectorsInHeat ( VectorsInHeat && other ) = delete;
105
111 VectorsInHeat & operator= ( const VectorsInHeat & other ) = delete;
112
118 VectorsInHeat & operator= ( VectorsInHeat && other ) = delete;
119
120
121
122 // ----------------------- Interface --------------------------------------
123
134 void init( double dt, double lambda = 1.0,
135 bool boundary_with_mixed_solution = false )
136 {
137 myIsInit=true;
138
139 SparseMatrix laplacian = myCalculus->globalLaplaceBeltrami( lambda );
140
141 SparseMatrix connectionLaplacian = myCalculus->globalConnectionLaplace( lambda );
142
143 SparseMatrix mass = myCalculus->globalLumpedMassMatrix();
144 SparseMatrix mass2= myCalculus->doubledGlobalLumpedMassMatrix();
145 myScalarHeatOpe = mass - dt*laplacian;
146 myVectorHeatOpe = mass2 - dt*connectionLaplacian;
147
148 //Prefactorizing
151
152 //empty sources
153 myVectorSource = Vector::Zero(2*myCalculus->nbVertices());
154 myScalarSource = Vector::Zero(myCalculus->nbVertices());
155 myDiracSource = Vector::Zero(myCalculus->nbVertices());
156
157 // Manage boundaries
158 myManageBoundary = false;
159 if ( ! boundary_with_mixed_solution ) return;
160 myBoundary = IntegerVector::Zero(myCalculus->nbVertices());
161 const auto surfmesh = myCalculus->getSurfaceMeshPtr();
162 const auto edges = surfmesh->computeManifoldBoundaryEdges();
163 for ( auto e : edges )
164 {
165 const auto vtcs = surfmesh->edgeVertices( e );
166 myBoundary[ vtcs.first ] = 1;
167 myBoundary[ vtcs.second ] = 1;
168 }
169 myManageBoundary = ! edges.empty();
170 if ( ! myManageBoundary ) return;
171 // Prepare solver for a problem with Dirichlet conditions.
173 // Prefactoring
174 myHeatDirichletSolver.compute( heatOpe_d );
175 }
176
182 void addSource(const Vertex aV,const Vector& ev)
183 {
184 ASSERT_MSG(aV < myCalculus->nbVertices(), "Vertex is not in the surface mesh vertex range");
185 Vector v = myCalculus->Tv(aV).transpose()*ev;
186 v = v.normalized()*ev.norm();
187 myVectorSource( 2*aV ) = v(0);
188 myVectorSource( 2*aV+1 ) = v(1);
189 myScalarSource( aV ) = v.norm();
190 myDiracSource( aV ) = 1;
191 }
192
196 {
197 myVectorSource.clear();
198 myScalarSource.clear();
199 myDiracSource.clear();
200 }
201
206 {
207 FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
208 return myVectorSource;
209 }
210
217 FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
218 return myCalculus->toExtrinsicVector(aV,intrinsicVectorSourceAtVertex(aV));
219 }
220
227 FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
228 Vector s(2);
229 s(0) = myVectorSource(2*aV);
230 s(1) = myVectorSource(2*aV+1);
231 return s;
232 }
233
234
237 std::vector<Vector> compute() const
238 {
239 FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
240 //Heat diffusion
241 Vector vectorHeatDiffusion = myVectorHeatSolver.solve(myVectorSource);
242 Vector scalarHeatDiffusion = myScalarHeatSolver.solve(myScalarSource);
243 Vector diracHeatDiffusion = myScalarHeatSolver.solve(myDiracSource);
244 auto surfmesh = myCalculus->getSurfaceMeshPtr();
245
246
247 // Take care of boundaries
248 if ( myManageBoundary )
249 {
250 Vector bValues = Vector::Zero( myCalculus->nbVertices() );
252 myBoundary, bValues );
253 Vector bSol = myHeatDirichletSolver.solve( bNormSources );
254 Vector heatDiffusionDirichlet
255 = Conditions::dirichletSolution( bSol, myBoundary, bValues );
256 scalarHeatDiffusion = 0.5 * ( scalarHeatDiffusion + heatDiffusionDirichlet );
257 }
258
259 std::vector<Vector> result(surfmesh->nbVertices());
260
261 for (auto v = 0;v<surfmesh->nbVertices();v++){
262 Vector Y(2);
263 Y(0) = vectorHeatDiffusion(2*v);
264 Y(1) = vectorHeatDiffusion(2*v+1);
265 Y = Y.normalized()*(scalarHeatDiffusion(v)/diracHeatDiffusion(v));
266 result[v] = myCalculus->toExtrinsicVector(v,Y);
267 }
268
269 return result;
270 }
271
272
274 bool isValid() const
275 {
276 return myIsInit && myCalculus->isValid();
277 }
278
279 // ----------------------- Private --------------------------------------
280
281private:
282
285
289
293
298
302
305
308
311
312}; // end of class VectorsInHeat
313} // namespace DGtal
314
315// //
317
318#endif // !defined VectorsInHeat_h
319
320#undef VectorsInHeat_RECURSES
321#endif // else defined(VectorsInHeat_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)
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.
This class implements on polygonal surfaces (using Discrete differential calculus on polygonal surfa...
Definition: VectorsInHeat.h:63
Vector extrinsicVectorSourceAtVertex(const Vertex aV)
extrinsicVectorSourceAtVertex get extrinsic source at vertex
Vector intrinsicVectorSourceAtVertex(const Vertex aV)
intrinsicVectorSourceAtVertex get intrinsic source at vertex
PolygonalCalculus::Vector Vector
Definition: VectorsInHeat.h:71
DirichletConditions< LinAlgBackend > Conditions
Definition: VectorsInHeat.h:74
void addSource(const Vertex aV, const Vector &ev)
PolygonalCalculus::LinAlg LinAlgBackend
Definition: VectorsInHeat.h:73
VectorsInHeat & operator=(const VectorsInHeat &other)=delete
Vector myScalarSource
Source vectors.
PolygonalCalculus::Vertex Vertex
Definition: VectorsInHeat.h:72
PolygonalCalculus::Solver Solver
Definition: VectorsInHeat.h:70
PolygonalCalculus::DenseMatrix DenseMatrix
Definition: VectorsInHeat.h:69
void init(double dt, double lambda=1.0, bool boundary_with_mixed_solution=false)
SparseMatrix myScalarHeatOpe
The operators for heat diffusion.
std::vector< Vector > compute() const
PolygonalCalculus::SparseMatrix SparseMatrix
Definition: VectorsInHeat.h:68
Vector vectorSource() const
VectorsInHeat(VectorsInHeat &&other)=delete
VectorsInHeat(const VectorsInHeat &other)=delete
bool isValid() const
IntegerVector myBoundary
The boundary characteristic vector.
bool myIsInit
Validitate flag.
Solver myScalarHeatSolver
Heat solvers.
TPolygonalCalculus PolygonalCalculus
Definition: VectorsInHeat.h:67
Conditions::IntegerVector IntegerVector
Definition: VectorsInHeat.h:75
VectorsInHeat(ConstAlias< PolygonalCalculus > calculus)
Definition: VectorsInHeat.h:84
SparseMatrix myVectorHeatOpe
const PolygonalCalculus * myCalculus
The underlying PolygonalCalculus instance.
Solver myHeatDirichletSolver
Heat solver with Dirichlet boundary conditions.
~VectorsInHeat()=default
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