TTK
Loading...
Searching...
No Matches
HarmonicField.cpp
Go to the documentation of this file.
1#include <HarmonicField.h>
2#include <set>
3
4#ifdef TTK_ENABLE_EIGEN
5#include <Eigen/Sparse>
6#endif // TTK_ENABLE_EIGEN
7
9 ttk::HarmonicField::findBestSolver(const SimplexId vertexNumber,
10 const SimplexId edgeNumber) const {
11
12 // for switching between Cholesky factorization and Iterate
13 // (conjugate gradients) method
14 const SimplexId threshold = 500000;
15
16 // compare threshold to number of non-zero values in laplacian matrix
17 if(2 * edgeNumber + vertexNumber > threshold) {
19 }
21}
22
23template <typename SparseMatrixType,
24 typename SparseVectorType,
25 typename SolverType>
26int ttk::HarmonicField::solve(SparseMatrixType const &lap,
27 SparseMatrixType const &penalty,
28 SparseVectorType const &constraints,
29 SparseMatrixType &sol) const {
30 SolverType solver(lap - penalty);
31 sol = solver.solve(penalty * constraints);
32 return solver.info();
33}
34
35// main routine
36template <class T, class TriangulationType>
37int ttk::HarmonicField::execute(const TriangulationType &triangulation,
38 const SimplexId constraintNumber,
39 const SimplexId *const sources,
40 const T *const constraints,
41 T *const outputScalarField,
42 const bool useCotanWeights,
43 const SolvingMethodUserType solvingMethod,
44 const double logAlpha) const {
45
46#ifdef TTK_ENABLE_EIGEN
47
48#ifdef TTK_ENABLE_OPENMP
49 Eigen::setNbThreads(threadNumber_);
50#endif // TTK_ENABLE_OPENMP
51
52 if(constraintNumber < 1) {
53 this->printErr("Cannot solve Laplace problem with no boundary constraints");
54 return 1;
55 }
56
57 using SpMat = Eigen::SparseMatrix<T>;
58 using SpVec = Eigen::SparseVector<T>;
59 using TripletType = Eigen::Triplet<T>;
60
61 Timer tm;
62
63 const auto vertexNumber = triangulation.getNumberOfVertices();
64 const auto edgeNumber = triangulation.getNumberOfEdges();
65
66 // find the right solving method
67 auto findSolvingMethod = [&]() -> SolvingMethodType {
69 switch(solvingMethod) {
70 case SolvingMethodUserType::AUTO:
71 res = findBestSolver(vertexNumber, edgeNumber);
72 break;
73 case SolvingMethodUserType::CHOLESKY:
74 res = SolvingMethodType::CHOLESKY;
75 break;
76 case SolvingMethodUserType::ITERATIVE:
77 res = SolvingMethodType::ITERATIVE;
78 break;
79 }
80 return res;
81 };
82
83 auto sm = findSolvingMethod();
84 std::string begMsg{"Beginning computation... ("};
85 if(useCotanWeights) {
86 begMsg.append("cotan weights, ");
87 } else {
88 begMsg.append("discrete laplacian, ");
89 }
90 if(sm == SolvingMethodType::ITERATIVE) {
91 begMsg.append("iterative method)");
92 } else {
93 begMsg.append("Cholesky method)");
94 }
95
96 this->printMsg(begMsg);
97
98 // filter unique constraint identifiers
99 std::set<SimplexId> uniqueIdentifiersSet;
100 for(SimplexId i = 0; i < constraintNumber; ++i) {
101 uniqueIdentifiersSet.insert(sources[i]);
102 }
103
104 // (constraint id, scalar field value)
105 std::vector<std::pair<SimplexId, T>> idValues{};
106
107 for(const auto id : uniqueIdentifiersSet) {
108 for(SimplexId j = 0; j < constraintNumber; ++j) {
109 if(id == sources[j]) {
110 idValues.emplace_back(id, constraints[j]);
111 break;
112 }
113 }
114 }
115
116 // unique constraint number
117 size_t uniqueConstraintNumber = idValues.size();
118
119 // graph laplacian of current mesh
120 SpMat lap;
121 if(useCotanWeights) {
122 Laplacian::cotanWeights<T>(lap, *this, triangulation);
123 } else {
124 Laplacian::discreteLaplacian<T>(lap, *this, triangulation);
125 }
126
127 // constraints vector
128 SpVec constraintsMat(vertexNumber);
129 for(const auto &pair : idValues) {
130 // put constraint at identifier index
131 constraintsMat.coeffRef(pair.first) = pair.second;
132 }
133
134 // penalty matrix
135 SpMat penalty(vertexNumber, vertexNumber);
136 // penalty value
137 const T alpha = Geometry::powIntTen(logAlpha);
138
139 std::vector<TripletType> triplets;
140 triplets.reserve(uniqueConstraintNumber);
141 for(const auto &pair : idValues) {
142 triplets.emplace_back(TripletType(pair.first, pair.first, alpha));
143 }
144#ifndef __clang_analyzer__
145 penalty.setFromTriplets(triplets.begin(), triplets.end());
146#endif // __clang_analyzer__
147
148 int res = 0;
149 SpMat sol;
150
151 switch(sm) {
152 case SolvingMethodType::CHOLESKY:
153 res = solve<SpMat, SpVec, Eigen::SimplicialCholesky<SpMat>>(
154 lap, penalty, constraintsMat, sol);
155 break;
156 case SolvingMethodType::ITERATIVE:
157 res = solve<SpMat, SpVec,
158 Eigen::ConjugateGradient<SpMat, Eigen::Upper | Eigen::Lower>>(
159 lap, penalty, constraintsMat, sol);
160 break;
161 }
162
163 auto info = static_cast<Eigen::ComputationInfo>(res);
164 switch(info) {
165 case Eigen::ComputationInfo::NumericalIssue:
166 this->printMsg("Numerical Issue!", ttk::debug::Priority::ERROR);
167 break;
168 case Eigen::ComputationInfo::NoConvergence:
169 this->printMsg("No Convergence!", ttk::debug::Priority::ERROR);
170 break;
171 case Eigen::ComputationInfo::InvalidInput:
172 this->printMsg("Invalid Input!", ttk::debug::Priority::ERROR);
173 break;
174 default:
175 break;
176 }
177
178 // convert to dense Eigen matrix
179 Eigen::Matrix<T, Eigen::Dynamic, 1> solDense(sol);
180
181 // copy solver solution into output array
182#ifdef TTK_ENABLE_OPENMP
183#pragma omp parallel for num_threads(threadNumber_)
184#endif // TTK_ENABLE_OPENMP
185 for(SimplexId i = 0; i < vertexNumber; ++i) {
186 // cannot avoid copy here...
187 outputScalarField[i] = -solDense(i, 0);
188 }
189
190 this->printMsg("Complete", 1.0, tm.getElapsedTime(), this->threadNumber_);
191
192#else
193 TTK_FORCE_USE(triangulation);
194 TTK_FORCE_USE(constraintNumber);
195 TTK_FORCE_USE(sources);
196 TTK_FORCE_USE(constraints);
197 TTK_FORCE_USE(outputScalarField);
198 TTK_FORCE_USE(useCotanWeights);
199 TTK_FORCE_USE(solvingMethod);
200 TTK_FORCE_USE(logAlpha);
201
202 this->printMsg(
203 std::vector<std::string>{
204 "Eigen support disabled, computation skipped!",
205 "Please re-compile TTK with Eigen support to enable this feature."},
207#endif // TTK_ENABLE_EIGEN
208
209 this->printMsg(ttk::debug::Separator::L1); // horizontal '=' separator
210 return 0;
211}
212
213// explicit template specializations for double and float types
214#define HARMONICFIELD_SPECIALIZE(TYPE) \
215 template int ttk::HarmonicField::execute<TYPE>( \
216 const Triangulation &, const SimplexId, const SimplexId *const, \
217 const TYPE *const, TYPE *const, const bool, const SolvingMethodUserType, \
218 const double) const
219
#define TTK_FORCE_USE(x)
Force the compiler to use the function/method parameter.
Definition: BaseClass.h:57
#define HARMONICFIELD_SPECIALIZE(TYPE)
int execute(const TriangulationType &triangulation, const SimplexId constraintNumber, const SimplexId *const sources, const T *const constraints, T *const outputScalarField, const bool useCotanWeights=true, const SolvingMethodUserType solvingMethod=SolvingMethodUserType::AUTO, const double logAlpha=5.0) const
double getElapsedTime()
Definition: Timer.h:15
T powIntTen(const int n)
Compute the nth power of ten.
Definition: Geometry.h:360
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
printMsg(debug::output::GREEN+"                           "+debug::output::ENDCOLOR+debug::output::GREEN+"▒"+debug::output::ENDCOLOR+debug::output::GREEN+"▒▒▒▒▒▒▒▒▒▒▒▒▒░"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)