DGtal  1.4.beta
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)
34 
35 #define VectorsInHeat_RECURSES
36 
37 #if !defined VectorsInHeat_h
38 
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 
49 namespace DGtal
50 {
52 // template class VectorInHeat
61 template <typename TPolygonalCalculus>
63 {
64  // ----------------------- Standard services ------------------------------
65 public:
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 
195  void clearSource()
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 
281 private:
282 
285 
289 
293 
298 
302 
305 
307  bool myIsInit;
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)
DGtal::VectorsInHeat::IntegerVector
Conditions::IntegerVector IntegerVector
Definition: VectorsInHeat.h:75
DGtal::VectorsInHeat::Conditions
DirichletConditions< LinAlgBackend > Conditions
Definition: VectorsInHeat.h:74
DGtal::VectorsInHeat::SparseMatrix
PolygonalCalculus::SparseMatrix SparseMatrix
Definition: VectorsInHeat.h:68
DGtal::VectorsInHeat::compute
std::vector< Vector > compute() const
Definition: VectorsInHeat.h:237
DGtal::VectorsInHeat::myManageBoundary
bool myManageBoundary
Definition: VectorsInHeat.h:301
DGtal::ConstAlias
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Definition: ConstAlias.h:186
DGtal::VectorsInHeat::operator=
VectorsInHeat & operator=(const VectorsInHeat &other)=delete
DGtal::VectorsInHeat::myVectorHeatOpe
SparseMatrix myVectorHeatOpe
Definition: VectorsInHeat.h:288
DGtal::DirichletConditions::dirichletVector
static DenseVector dirichletVector(const SparseMatrix &A, const DenseVector &b, const IntegerVector &p, const DenseVector &u)
Definition: DirichletConditions.h:170
DGtal::PolygonalCalculus::SparseMatrix
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
Definition: PolygonalCalculus.h:102
DGtal::VectorsInHeat::myHeatDirichletSolver
Solver myHeatDirichletSolver
Heat solver with Dirichlet boundary conditions.
Definition: VectorsInHeat.h:310
DGtal::VectorsInHeat::myDiracSource
Vector myDiracSource
Definition: VectorsInHeat.h:296
DGtal::VectorsInHeat::PolygonalCalculus
TPolygonalCalculus PolygonalCalculus
Definition: VectorsInHeat.h:67
DGtal::VectorsInHeat::myScalarHeatSolver
Solver myScalarHeatSolver
Heat solvers.
Definition: VectorsInHeat.h:291
DGtal::SurfaceMesh::computeManifoldBoundaryEdges
Edges computeManifoldBoundaryEdges() const
DGtal::VectorsInHeat::addSource
void addSource(const Vertex aV, const Vector &ev)
Definition: VectorsInHeat.h:182
DGtal::DirichletConditions::dirichletSolution
static DenseVector dirichletSolution(const DenseVector &xd, const IntegerVector &p, const DenseVector &u)
Definition: DirichletConditions.h:203
DGtal::DirichletConditions::dirichletOperator
static SparseMatrix dirichletOperator(const SparseMatrix &A, const IntegerVector &p)
Definition: DirichletConditions.h:130
DGtal::VectorsInHeat::init
void init(double dt, double lambda=1.0, bool boundary_with_mixed_solution=false)
Definition: VectorsInHeat.h:134
DGtal::VectorsInHeat::intrinsicVectorSourceAtVertex
Vector intrinsicVectorSourceAtVertex(const Vertex aV)
intrinsicVectorSourceAtVertex get intrinsic source at vertex
Definition: VectorsInHeat.h:226
DGtal::SurfaceMesh::nbVertices
Size nbVertices() const
Definition: SurfaceMesh.h:280
DGtal::PolygonalCalculus::DenseMatrix
LinAlg::DenseMatrix DenseMatrix
Type of dense matrix.
Definition: PolygonalCalculus.h:100
laplacian
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)
Definition: sphereCotangentLaplaceOperator.cpp:87
DGtal::VectorsInHeat
This class implements on polygonal surfaces (using Discrete differential calculus on polygonal surfa...
Definition: VectorsInHeat.h:62
DGtal::VectorsInHeat::myScalarHeatOpe
SparseMatrix myScalarHeatOpe
The operators for heat diffusion.
Definition: VectorsInHeat.h:287
DGtal::VectorsInHeat::~VectorsInHeat
~VectorsInHeat()=default
DGtal::VectorsInHeat::VectorsInHeat
VectorsInHeat(ConstAlias< PolygonalCalculus > calculus)
Definition: VectorsInHeat.h:84
DGtal::VectorsInHeat::myBoundary
IntegerVector myBoundary
The boundary characteristic vector.
Definition: VectorsInHeat.h:304
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::EigenLinearAlgebraBackend
Aim: Provide linear algebra backend using Eigen dense and sparse matrix as well as dense vector....
Definition: EigenSupport.h:96
DGtal::SurfaceMesh::edgeVertices
const VertexPair & edgeVertices(Edge e) const
Definition: SurfaceMesh.h:329
DGtal::VectorsInHeat::vectorSource
Vector vectorSource() const
Definition: VectorsInHeat.h:205
DGtal::PolygonalCalculus::Solver
LinAlg::SolverSimplicialLDLT Solver
Type of a sparse matrix solver.
Definition: PolygonalCalculus.h:107
DGtal::VectorsInHeat::LinAlgBackend
PolygonalCalculus::LinAlg LinAlgBackend
Definition: VectorsInHeat.h:73
DGtal::VectorsInHeat::Vector
PolygonalCalculus::Vector Vector
Definition: VectorsInHeat.h:71
DGtal::VectorsInHeat::myIsInit
bool myIsInit
Validitate flag.
Definition: VectorsInHeat.h:307
DGtal::DirichletConditions::IntegerVector
LinearAlgebraBackend::IntegerVector IntegerVector
Definition: DirichletConditions.h:105
DGtal::VectorsInHeat::myCalculus
const PolygonalCalculus * myCalculus
The underlying PolygonalCalculus instance.
Definition: VectorsInHeat.h:284
DGtal::VectorsInHeat::myVectorSource
Vector myVectorSource
Definition: VectorsInHeat.h:297
DGtal::VectorsInHeat::clearSource
void clearSource()
Definition: VectorsInHeat.h:195
DGtal::VectorsInHeat::Solver
PolygonalCalculus::Solver Solver
Definition: VectorsInHeat.h:70
DGtal::PolygonalCalculus::Vertex
MySurfaceMesh::Vertex Vertex
Vertex type.
Definition: PolygonalCalculus.h:85
DGtal::DirichletConditions
Aim: A helper class to solve a system with Dirichlet boundary conditions.
Definition: DirichletConditions.h:97
DGtal::PolygonalCalculus::Vector
LinAlg::DenseVector Vector
Type of Vector.
Definition: PolygonalCalculus.h:96
DGtal::VectorsInHeat::Vertex
PolygonalCalculus::Vertex Vertex
Definition: VectorsInHeat.h:72
DGtal::VectorsInHeat::VectorsInHeat
VectorsInHeat()=delete
DGtal::VectorsInHeat::myVectorHeatSolver
Solver myVectorHeatSolver
Definition: VectorsInHeat.h:292
DGtal::VectorsInHeat::extrinsicVectorSourceAtVertex
Vector extrinsicVectorSourceAtVertex(const Vertex aV)
extrinsicVectorSourceAtVertex get extrinsic source at vertex
Definition: VectorsInHeat.h:216
DGtal::VectorsInHeat::myScalarSource
Vector myScalarSource
Source vectors.
Definition: VectorsInHeat.h:295
DGtal::VectorsInHeat::isValid
bool isValid() const
Definition: VectorsInHeat.h:274
DGtal::VectorsInHeat::DenseMatrix
PolygonalCalculus::DenseMatrix DenseMatrix
Definition: VectorsInHeat.h:69