Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
RBFWithQuadratic.cpp
Go to the documentation of this file.
1
3
7
8#include <igl/Timer.h>
9
10#include <Eigen/Dense>
11
12#include <iostream>
13#include <fstream>
14#include <array>
16
17// #define VERBOSE
18
19using namespace polyfem;
20using namespace polyfem::assembler;
21using namespace polyfem::basis;
22using namespace polyfem::quadrature;
23using namespace polyfem::utils;
24
25namespace
26{
27
28 // Harmonic kernel
29 double kernel(const bool is_volume, const double r)
30 {
31 if (r < 1e-8)
32 {
33 return 0;
34 }
35
36 if (is_volume)
37 {
38 return 1 / r;
39 }
40 else
41 {
42 return log(r);
43 }
44 }
45
46 double kernel_prime(const bool is_volume, const double r)
47 {
48 if (r < 1e-8)
49 {
50 return 0;
51 }
52
53 if (is_volume)
54 {
55 return -1 / (r * r);
56 }
57 else
58 {
59 return 1 / r;
60 }
61 }
62
63 // Biharmonic kernel (2d only)
64 // double kernel(const bool is_volume, const double r) {
65 // assert(!is_volume);
66 // if (r < 1e-8) { return 0; }
67
68 // return r * r * (log(r)-1);
69 // }
70
71 // double kernel_prime(const bool is_volume, const double r) {
72 // assert(!is_volume);
73 // if (r < 1e-8) { return 0; }
74
75 // return r * ( 2 * log(r) - 1);
76 // }
77
78} // anonymous namespace
79
81
82// output is std::array<Eigen::MatrixXd, 5> &strong rhs(q(x_i) er)
83void RBFWithQuadratic::setup_monomials_strong_2d(const int dim, const LinearAssembler &assembler, const Eigen::MatrixXd &pts, const QuadratureVector &da, std::array<Eigen::MatrixXd, 5> &strong)
84{
85 // a(u,v) = a(q er, phi_j es) = <rhs(q(x_i) er) , phi_j(x_i) es >
86 // (not a(phi_j es, q er))
87
89 AutodiffHessianPt pt(dim);
90 for (int i = 0; i < 5; ++i)
91 {
92 strong[i].resize(dim * dim, pts.rows());
93 strong[i].setZero();
94 }
95
96 Eigen::MatrixXd tmp;
97
98 for (int i = 0; i < pts.rows(); ++i)
99 {
100 // loop for er
101 for (int d = 0; d < dim; ++d)
102 {
103 pt((d + 1) % dim) = AutodiffScalarHessian(0);
104 // for d = 0 pt(q, 0), for d = 1 pt=(0, q)
105
106 // x
107 pt(d) = AutodiffScalarHessian(0, pts(i, 0)); // pt=(x, 0) or pt=(0, x)
108 tmp = assembler.compute_rhs(pt); // in R^dim
109 for (int d1 = 0; d1 < dim; ++d1)
110 strong[0](d * dim + d1, i) = tmp(d1) * da(i);
111
112 // y
113 pt(d) = AutodiffScalarHessian(1, pts(i, 1));
114 tmp = assembler.compute_rhs(pt);
115 for (int d1 = 0; d1 < dim; ++d1)
116 strong[1](d * dim + d1, i) = tmp(d1) * da(i);
117
118 // xy
119 pt(d) = AutodiffScalarHessian(0, pts(i, 0)) * AutodiffScalarHessian(1, pts(i, 1));
120 tmp = assembler.compute_rhs(pt);
121 for (int d1 = 0; d1 < dim; ++d1)
122 strong[2](d * dim + d1, i) = tmp(d1) * da(i);
123
124 // x^2
125 pt(d) = AutodiffScalarHessian(0, pts(i, 0)) * AutodiffScalarHessian(0, pts(i, 0));
126 tmp = assembler.compute_rhs(pt);
127 for (int d1 = 0; d1 < dim; ++d1)
128 strong[3](d * dim + d1, i) = tmp(d1) * da(i);
129
130 // y^2
131 pt(d) = AutodiffScalarHessian(1, pts(i, 1)) * AutodiffScalarHessian(1, pts(i, 1));
132 tmp = assembler.compute_rhs(pt);
133 for (int d1 = 0; d1 < dim; ++d1)
134 strong[4](d * dim + d1, i) = tmp(d1) * da(i);
135 }
136 }
137}
138
139void RBFWithQuadratic::setup_monomials_vals_2d(const int star_index, const Eigen::MatrixXd &pts, ElementAssemblyValues &vals)
140{
141 assert(star_index + 5 <= vals.basis_values.size());
142 // x
143 vals.basis_values[star_index + 0].val = pts.col(0);
144 vals.basis_values[star_index + 0].grad = Eigen::MatrixXd(pts.rows(), pts.cols());
145 vals.basis_values[star_index + 0].grad.col(0).setOnes();
146 vals.basis_values[star_index + 0].grad.col(1).setZero();
147
148 // y
149 vals.basis_values[star_index + 1].val = pts.col(1);
150 vals.basis_values[star_index + 1].grad = Eigen::MatrixXd(pts.rows(), pts.cols());
151 vals.basis_values[star_index + 1].grad.col(0).setZero();
152 vals.basis_values[star_index + 1].grad.col(1).setOnes();
153
154 // xy
155 vals.basis_values[star_index + 2].val = pts.col(0).array() * pts.col(1).array();
156 vals.basis_values[star_index + 2].grad = Eigen::MatrixXd(pts.rows(), pts.cols());
157 vals.basis_values[star_index + 2].grad.col(0) = pts.col(1);
158 vals.basis_values[star_index + 2].grad.col(1) = pts.col(0);
159
160 // x^2
161 vals.basis_values[star_index + 3].val = pts.col(0).array() * pts.col(0).array();
162 vals.basis_values[star_index + 3].grad = Eigen::MatrixXd(pts.rows(), pts.cols());
163 vals.basis_values[star_index + 3].grad.col(0) = 2 * pts.col(0);
164 vals.basis_values[star_index + 3].grad.col(1).setZero();
165
166 // y^2
167 vals.basis_values[star_index + 4].val = pts.col(1).array() * pts.col(1).array();
168 vals.basis_values[star_index + 4].grad = Eigen::MatrixXd(pts.rows(), pts.cols());
169 vals.basis_values[star_index + 4].grad.col(0).setZero();
170 vals.basis_values[star_index + 4].grad.col(1) = 2 * pts.col(1);
171
172 for (size_t i = star_index; i < star_index + 5; ++i)
173 {
174 vals.basis_values[i].grad_t_m = vals.basis_values[i].grad;
175 }
176
177 // for(size_t i = star_index; i < star_index + 5; ++i)
178 // {
179 // vals.basis_values[i].grad_t_m = Eigen::MatrixXd(pts.rows(), pts.cols());
180 // for(int k = 0; k < vals.jac_it.size(); ++k)
181 // vals.basis_values[i].grad_t_m.row(k) = vals.basis_values[i].grad.row(k) * vals.jac_it[k];
182 // }
183}
184
186 const LinearAssembler &assembler,
187 const Eigen::MatrixXd &centers,
188 const Eigen::MatrixXd &collocation_points,
189 const Eigen::MatrixXd &local_basis_integral,
190 const Quadrature &quadr,
191 Eigen::MatrixXd &rhs,
192 bool with_constraints)
193 : centers_(centers)
194{
195 // centers_.resize(0, centers.cols());
196 compute_weights(assembler, collocation_points, local_basis_integral, quadr, rhs, with_constraints);
197}
198
199// -----------------------------------------------------------------------------
200
201void RBFWithQuadratic::basis(const int local_index, const Eigen::MatrixXd &samples, Eigen::MatrixXd &val) const
202{
203 Eigen::MatrixXd tmp;
204 bases_values(samples, tmp);
205 val = tmp.col(local_index);
206}
207
208// -----------------------------------------------------------------------------
209
210void RBFWithQuadratic::grad(const int local_index, const Eigen::MatrixXd &samples, Eigen::MatrixXd &val) const
211{
212 Eigen::MatrixXd tmp;
213 const int dim = centers_.cols();
214 val.resize(samples.rows(), dim);
215 for (int d = 0; d < dim; ++d)
216 {
217 bases_grads(d, samples, tmp);
218 val.col(d) = tmp.col(local_index);
219 }
220}
221
223
224void RBFWithQuadratic::bases_values(const Eigen::MatrixXd &samples, Eigen::MatrixXd &val) const
225{
226 // Compute A
227 Eigen::MatrixXd A;
228 compute_kernels_matrix(samples, A);
229
230 // Multiply by the weights
231 val = A * weights_;
232}
233
234// -----------------------------------------------------------------------------
235
236void RBFWithQuadratic::bases_grads(const int axis, const Eigen::MatrixXd &samples, Eigen::MatrixXd &val) const
237{
238 const int num_kernels = centers_.rows();
239 const int dim = (is_volume() ? 3 : 2);
240
241 // Compute ∇xA
242 Eigen::MatrixXd A_prime(samples.rows(), num_kernels + 1 + dim + dim * (dim + 1) / 2);
243 A_prime.setZero();
244
245 for (int j = 0; j < num_kernels; ++j)
246 {
247 A_prime.col(j) = (samples.rowwise() - centers_.row(j)).rowwise().norm().unaryExpr([this](double x) { return kernel_prime(is_volume(), x) / x; });
248 A_prime.col(j) = (samples.col(axis).array() - centers_(j, axis)) * A_prime.col(j).array();
249 }
250 // Linear terms
251 A_prime.middleCols(num_kernels + 1 + axis, 1).setOnes();
252 // Mixed terms
253 if (dim == 2)
254 {
255 A_prime.col(num_kernels + 1 + dim) = samples.col(1 - axis);
256 }
257 else
258 {
259 A_prime.col(num_kernels + 1 + dim + axis) = samples.col((axis + 1) % dim);
260 A_prime.col(num_kernels + 1 + dim + (axis + 2) % dim) = samples.col((axis + 2) % dim);
261 }
262 // Quadratic terms
263 A_prime.rightCols(dim).col(axis) = 2.0 * samples.col(axis);
264
265 // Apply weights
266 val = A_prime * weights_;
267}
268
270//
271// For each FEM basis φ that is nonzero on the element E, we want to
272// solve the least square system A w = rhs, where:
273// ┏ ┓
274// ┃ ψ_k(pi) ... 1 xi yi xi*yi xi^2 yi^2 ┃
275// A = ┃ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┃ ∊ ℝ^{#S x (#K+1+dim+dim*(dim+1)/2)}
276// ┃ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┃
277// ┗ ┛
278// ┏ ┓^⊤
279// w = ┃ w_k ... a00 a10 a01 a11 a20 a02 ┃ ∊ ℝ^{#K+1+dim+dim*(dim+1)/2}
280// ┗ ┛
281// - A is the RBF kernels evaluated over the collocation points (#S)
282// - b is the expected value of the basis sampled on the boundary (#S)
283// - w is the weight of the kernels defining the basis
284// - pi = (xi, yi) is the i-th collocation point
285//
286// Moreover, we want to impose a constraint on the weight vector w so that each
287// monomial Q(x,y) = x^α*y^β with α+β <= 2 is in the span of the FEM bases {φ_j}_j.
288//
289// In the case of Laplace's equation, we recall the weak form of the PDE as:
290//
291// Find u such that: ∫_Ω Δu v = - ∫_Ω ∇u·∇v ∀ v
292//
293// For our bases to exactly represent a monomial Q(x,y), it means that its
294// approximation by the finite element bases {φ_j}_j must be equal to Q(x,y).
295// In particular, for any φ_j that is nonzero on the polyhedral element E, we must have:
296//
297// ∫_{𝘅 in Ω} ΔQ(𝘅) φ_j(𝘅) d𝘅 = - ∫_{𝘅 \in Ω} ∇Q(𝘅)·∇φ_j(𝘅) d𝘅 (1)
298//
299// Now, for each of the 5 non-constant monomials (9 in 3D), we need to compute
300// Δ(x^α*y^β). For (α,β) ∊ {(1,0), (0,1), (1,1), (2,0), (0,2)}, this yields
301// the following equalities:
302//
303// Δx = 0 (2a)
304// Δy = 0 (2b)
305// Δxy = 0 (2c)
306// Δx² = 1 (2d)
307// Δy² = 1 (2e)
308//
309// If we plug these back into (1), and split the integral between the polyhedral
310// element E and Ω\E, we obtain the following constraints:
311//
312// ∫_E ∇Q·∇φ_j + ∫_E ΔQ φ_j = - ∫_{Ω\E} ∇Q·∇φ_j - ∫_{Ω\E} ΔQ φ_j (3)
313//
314// Note that the right-hand side of (3) is already known, since no two polyhedral
315// cells are adjacent to each other, and the bases overlapping a polyhedron vanish
316// on the boundary of the domain ∂Ω. This right-hand side is computed in advance
317// and passed to our functions as in argument `local_basis_integral`.
318//
319// The left-hand side of equation (3) reduces to the following (in 2D):
320//
321// ∫_E ∇x(φ_j) = c10 (4a)
322// ∫_E ∇y(φ_j) = c01 (4b)
323// ∫_E (y·∇x(φ_j) + y·∇x(φ_jj)) = c11 (4c)
324// ∫_E 2x·∇x(φ_j) + ∫_E 2 φ_j = c20 (4d)
325// ∫_E 2y·∇y(φ_j) + ∫_E 2 φ_j = c02 (4e)
326//
327// The next step is to express the basis φ_j in terms of the harmonic kernels and
328// quadratic polynomials:
329//
330// φ_j(x,y) = Σ_k w_k ψ_k(x,y) + a00 + a10*x + a01*y + a11*x*y + a20*x² + a02*y²
331//
332// The five equations in (4) become:
333//
334// Σ_j w_k ∫∇x(ψ_k) = ∫ Δ q10 (Σ_j w_k (ψ_k) + a00) + Σ_j w_k ∫∇q10 . ∇(ψ_k + a00)
335// Σ_j w_k ∫∇x(ψ_k) + a10 |E| + a11 ∫y + a20 ∫2x = c10
336// Σ_j w_k ∫∇y(ψ_k) + a01 |E| + a11 ∫x + a02 ∫2y = c01
337// Σ_j w_k (∫y·∇x(ψ_k) + ∫x·∇y(ψ_k)) + a10 ∫y + a01 ∫x + a11 (∫x²+∫y²) + a20 2∫xy + a02 2∫xy = c11
338// Σ_j w_k (2∫x·∇x(ψ_k) + 2ψ_k) + a10 4∫x + a01 2∫y + a11 4∫xy + a20 6∫x² + a02 2∫y² = c20
339// Σ_j w_k (2∫y·∇y(ψ_k) + 2ψ_k) + a10 2∫x + a01 4∫y + a11 4∫xy + a20 2∫x² + a02 6∫y² = c02
340// Σ_j w_k (2∫y·∇y(ψ_k) + 2ψ_k) = ∫ Δ q20 (Σ_j w_k (ψ_k) + a00) + Σ_j w_k ∫∇q20 . ∇(ψ_k + a00) = ∫ -2 (Σ_j w_k (ψ_k) + a00) + Σ_j w_k ∫2x ∇x(ψ_k)
341//
342// This system gives us a relationship between the fives a10, a01, a11, a20, a02
343// and the rest of the w_k + a constant translation term. We can write down the
344// corresponding system:
345//
346// a10 a01 a11 a20 a02
347// ┏ ┓ ┏ ┓
348// ┃ |E| ∫y 2∫x ┃ ┃ w_k ┃
349// ┃ ┃ ┃ ┊ ┃
350// ┃ |E| ∫x 2∫y ┃ ┃ ┊ ┃
351// ┃ ┃ ┃ ┊ ┃
352// M = ┃ ∫y ∫x ∫x²+∫y² 2∫xy 2∫xy ┃ = \tilde{L} ┃ ┊ ┃ + \tilde{t}
353// ┃ ┃ ┃ ┊ ┃
354// ┃ 4∫x 2∫y 4∫xy 6∫x² 2∫y² ┃ ┃ ┊ ┃
355// ┃ ┃ ┃w_#K ┃
356// ┃ 2∫x 4∫y 4∫xy 2∫x² 6∫y² ┃ ┃ a00 ┃
357// ┗ ┛ ┗ ┛
358//
359// Now, if we want to express w as w = Lv + t, and solve our least-square
360// system as before, we need to invert M and compute L and t in terms of
361// \tilde{L} and \tilde{t}
362//
363// ┏ ┓
364// ┃ 1 ┃
365// ┃ 1 ┃
366// ┃ · ┃
367// L = ┃ · ┃ ∊ ℝ^{ (#K+1+dim+dim*(dim+1)/2) x (#K+1}) }
368// ┃ 1 ┃
369// ┃ M^{-1} \tilde{L} ┃
370// ┗ ┛
371// ┏ ┓
372// ┃ 0 ┃
373// ┃ ┊ ┃
374// t = ┃ ┊ ┃ ∊ ℝ^{#K+1+dim+dim*(dim+1)/2}
375// ┃ 0 ┃
376// ┃ M^{-1} \tilde{t} ┃
377// ┗ ┛
378// After solving the new least square system A L v = rhs - A t, we can retrieve
379// w = L v
380//
382
383void RBFWithQuadratic::compute_kernels_matrix(const Eigen::MatrixXd &samples, Eigen::MatrixXd &A) const
384{
385 // Compute A
386 const int num_kernels = centers_.rows();
387 const int dim = (is_volume() ? 3 : 2);
388
389 A.resize(samples.rows(), num_kernels + 1 + dim + dim * (dim + 1) / 2);
390 for (int j = 0; j < num_kernels; ++j)
391 {
392 A.col(j) = (samples.rowwise() - centers_.row(j)).rowwise().norm().unaryExpr([this](double x) { return kernel(is_volume(), x); });
393 }
394 A.col(num_kernels).setOnes(); // constant term
395 A.middleCols(num_kernels + 1, dim) = samples; // linear terms
396 if (dim == 2)
397 {
398 A.middleCols(num_kernels + dim + 1, 1) = samples.rowwise().prod(); // mixed terms
399 }
400 else if (dim == 3)
401 {
402 A.middleCols(num_kernels + dim + 1, 3) = samples;
403 A.middleCols(num_kernels + dim + 1 + 0, 1).array() *= samples.col(1).array();
404 A.middleCols(num_kernels + dim + 1 + 1, 1).array() *= samples.col(2).array();
405 A.middleCols(num_kernels + dim + 1 + 2, 1).array() *= samples.col(0).array();
406 }
407 A.rightCols(dim) = samples.array().square(); // quadratic terms
408}
409
410// -----------------------------------------------------------------------------
411
413 const int num_bases,
414 const Quadrature &quadr,
415 const Eigen::MatrixXd &local_basis_integral,
416 Eigen::MatrixXd &L,
417 Eigen::MatrixXd &t) const
418{
419 const int num_kernels = centers_.rows();
420 const int dim = centers_.cols();
421 assert(dim == 2);
422
423 // K_cst = ∫ψ_k
424 // K_lin = ∫∇x(ψ_k), ∫∇y(ψ_k)
425 // K_mix = ∫y·∇x(ψ_k), ∫x·∇y(ψ_k)
426 // K_sqr = ∫x·∇x(ψ_k), ∫y·∇y(ψ_k)
427 Eigen::VectorXd K_cst = Eigen::VectorXd::Zero(num_kernels);
428 Eigen::MatrixXd K_lin = Eigen::MatrixXd::Zero(num_kernels, dim);
429 Eigen::MatrixXd K_mix = Eigen::MatrixXd::Zero(num_kernels, dim);
430 Eigen::MatrixXd K_sqr = Eigen::MatrixXd::Zero(num_kernels, dim);
431 for (int j = 0; j < num_kernels; ++j)
432 {
433 // ∫∇x(ψ_k)(p) = Σ_q (xq - xk) * 1/r * h'(r) * wq
434 // - xq is the x coordinate of the q-th quadrature point
435 // - wq is the q-th quadrature weight
436 // - r is the distance from pq to the kernel center
437 // - h is the RBF kernel (scalar function)
438 for (int q = 0; q < quadr.points.rows(); ++q)
439 {
440 const RowVectorNd p = quadr.points.row(q) - centers_.row(j);
441 const double r = p.norm();
442 const RowVectorNd gradPhi = p * kernel_prime(is_volume(), r) / r * quadr.weights(q);
443 K_cst(j) += kernel(is_volume(), r) * quadr.weights(q);
444 K_lin.row(j) += gradPhi;
445 K_mix(j, 0) += quadr.points(q, 1) * gradPhi(0);
446 K_mix(j, 1) += quadr.points(q, 0) * gradPhi(1);
447 K_sqr.row(j) += (quadr.points.row(q).array() * gradPhi.array()).matrix();
448 }
449 }
450
451 // I_lin = ∫x, ∫y
452 // I_mix = ∫xy
453 // I_sqr = ∫x², ∫y²
454 Eigen::RowVectorXd I_lin = (quadr.points.array().colwise() * quadr.weights.array()).colwise().sum();
455 Eigen::RowVectorXd I_mix = (quadr.points.rowwise().prod().array() * quadr.weights.array()).colwise().sum();
456 Eigen::RowVectorXd I_sqr = (quadr.points.array().square().colwise() * quadr.weights.array()).colwise().sum();
457 double volume = quadr.weights.sum();
458
459 // Compute M
460 Eigen::Matrix<double, 5, 5> M;
461 M << volume, 0, I_lin(1), 2 * I_lin(0), 0,
462 0, volume, I_lin(0), 0, 2 * I_lin(1),
463 I_lin(1), I_lin(0), I_sqr(0) + I_sqr(1), 2 * I_mix(0), 2 * I_mix(0),
464 4 * I_lin(0), 2 * I_lin(1), 4 * I_mix(0), 6 * I_sqr(0), 2 * I_sqr(1),
465 2 * I_lin(0), 4 * I_lin(1), 4 * I_mix(0), 2 * I_sqr(0), 6 * I_sqr(1);
466 Eigen::FullPivLU<Eigen::Matrix<double, 5, 5>> lu(M);
467 assert(lu.isInvertible());
468
469 // show_matrix_stats(M);
470
471 // Compute L
472 L.resize(num_kernels + 1 + dim + dim * (dim + 1) / 2, num_kernels + 1);
473 L.setZero();
474 L.diagonal().setOnes();
475
476 L.block(num_kernels + 1, 0, dim, num_kernels) = -K_lin.transpose();
477 L.block(num_kernels + 1 + dim, 0, 1, num_kernels) = -K_mix.transpose().colwise().sum();
478 L.block(num_kernels + 1 + dim + 1, 0, dim, num_kernels) = -2.0 * (K_sqr.colwise() + K_cst).transpose();
479 L.bottomRightCorner(dim, 1).setConstant(-2.0 * volume);
480 // j \in [0, 4]
481 // i \in [0, num_kernels]
482 // ass_val = [q_10, q_01, q_11, q_20, q_02, psi_0, ..., psi_k]
483
484 // strong rows is the evaluation at quadrature points
485 // strong.col(0) = pde(q_10) (probably 0)
486 // strong.col(4) = pde(q_02) (it is 2 for laplacian)
487 // L.block(num_kernels + 1 + i, j) = +/- assembler.assemble(ass_val, j, 5 + i) +/- (strong.col(j).array() * ass_val.basis_values[5+i].val.array() * quadr.weights.array()).sum();
488
489 L.block(num_kernels + 1, 0, 5, num_kernels + 1) = lu.solve(L.block(num_kernels + 1, 0, 5, num_kernels + 1));
490
491 // Compute t
492 t.resize(L.rows(), num_bases);
493 t.setZero();
494 t.bottomRows(5) = local_basis_integral.transpose();
495 t.bottomRows(5) = lu.solve(weights_.bottomRows(5));
496}
497
499 const LinearAssembler &assembler,
500 const int num_bases,
501 const Quadrature &quadr,
502 const Eigen::MatrixXd &local_basis_integral,
503 Eigen::MatrixXd &L,
504 Eigen::MatrixXd &t) const
505{
506 // TODO
507 const double time = 0;
508 const int num_kernels = centers_.rows();
509 const int space_dim = centers_.cols();
510 const int assembler_dim = assembler.is_tensor() ? 2 : 1;
511 assert(space_dim == 2);
512
513 std::array<Eigen::MatrixXd, 5> strong;
514
515 // ass_val = [q_10, q_01, q_11, q_20, q_02, psi_0, ..., psi_k]
516 ElementAssemblyValues ass_val;
517 ass_val.has_parameterization = false;
518 ass_val.basis_values.resize(5 + num_kernels);
519
520 // evaluating monomial and grad of monomials at quad points
521 setup_monomials_vals_2d(0, quadr.points, ass_val);
522 setup_monomials_strong_2d(assembler_dim, assembler, quadr.points, quadr.weights.array(), strong);
523
524 // evaluating psi and grad psi at quadr points
525 for (int j = 0; j < num_kernels; ++j)
526 {
527 ass_val.basis_values[5 + j].val = Eigen::MatrixXd(quadr.points.rows(), 1);
528 ass_val.basis_values[5 + j].grad = Eigen::MatrixXd(quadr.points.rows(), quadr.points.cols());
529
530 for (int q = 0; q < quadr.points.rows(); ++q)
531 {
532 const RowVectorNd p = quadr.points.row(q) - centers_.row(j);
533 const double r = p.norm();
534
535 ass_val.basis_values[5 + j].val(q) = kernel(is_volume(), r);
536 ass_val.basis_values[5 + j].grad.row(q) = p * kernel_prime(is_volume(), r) / r;
537 }
538 }
539
540 for (size_t i = 5; i < ass_val.basis_values.size(); ++i)
541 {
542 ass_val.basis_values[i].grad_t_m = ass_val.basis_values[i].grad;
543 }
544
545 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 10, 10> M(5 * assembler_dim, 5 * assembler_dim);
546 for (int i = 0; i < 5; ++i)
547 {
548 for (int j = 0; j < 5; ++j)
549 {
550 const auto tmp = assembler.assemble(LinearAssemblerData(ass_val, time, i, j, quadr.weights));
551
552 for (int d1 = 0; d1 < assembler_dim; ++d1)
553 {
554 for (int d2 = 0; d2 < assembler_dim; ++d2)
555 {
556 const int loc_index = d1 * assembler_dim + d2;
557 M(i * assembler_dim + d1, j * assembler_dim + d2) = tmp(loc_index) + (strong[i].row(loc_index).transpose().array() * ass_val.basis_values[j].val.array()).sum();
558 }
559 }
560 }
561 }
562
563 Eigen::FullPivLU<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 10, 10>> lu(M);
564 assert(lu.isInvertible());
565
566 // Compute L
567 L.resize((num_kernels + 1 + space_dim + space_dim * (space_dim + 1) / 2) * assembler_dim, (num_kernels + 1) * assembler_dim);
568 L.setZero();
569 L.diagonal().setOnes();
570
571 for (int i = 0; i < 5; ++i)
572 {
573 for (int j = 0; j < num_kernels; ++j)
574 {
575 const auto tmp = assembler.assemble(LinearAssemblerData(ass_val, time, i, 5 + j, quadr.weights));
576 for (int d1 = 0; d1 < assembler_dim; ++d1)
577 {
578 for (int d2 = 0; d2 < assembler_dim; ++d2)
579 {
580 const int loc_index = d1 * assembler_dim + d2;
581 L((num_kernels + 1 + i) * assembler_dim + d1, j * assembler_dim + d2) = -tmp(loc_index) - (strong[i].row(loc_index).transpose().array() * ass_val.basis_values[5 + j].val.array()).sum();
582 // L(num_kernels + 1 + i*assembler_dim + d1, j*assembler_dim + d2) = -assembler.local_assemble(ass_val, i, 5 + j, quadr.weights)(0) - (strong[i].transpose().array() * ass_val.basis_values[5+j].val.array()).sum();
583 }
584 }
585 }
586 for (int d1 = 0; d1 < assembler_dim; ++d1)
587 {
588 for (int d2 = 0; d2 < assembler_dim; ++d2)
589 {
590 const int loc_index = d1 * assembler_dim + d2;
591 L(num_kernels + 1 + i * assembler_dim + d1, assembler_dim * num_kernels + d2) = -strong[i].row(loc_index).sum();
592 }
593 }
594
595 // L(num_kernels + 1 + i*assembler_dim + d1, assembler_dim*num_kernels) = - strong[i].sum();
596 }
597
598 L.block((num_kernels + 1) * assembler_dim, 0, 5 * assembler_dim, (num_kernels + 1) * assembler_dim) = lu.solve(L.block((num_kernels + 1) * assembler_dim, 0, 5 * assembler_dim, (num_kernels + 1) * assembler_dim));
599
600 // Compute t
601 // t == weights_
602 t.resize(L.rows(), num_bases * assembler_dim);
603 t.setZero();
604 t.bottomRows(5 * assembler_dim) = lu.solve(local_basis_integral.transpose());
605}
606
607// -----------------------------------------------------------------------------
608
610 const LinearAssembler &assembler,
611 const int num_bases,
612 const Quadrature &quadr,
613 const Eigen::MatrixXd &local_basis_integral,
614 Eigen::MatrixXd &L,
615 Eigen::MatrixXd &t) const
616{
617 const int num_kernels = centers_.rows();
618 const int dim = centers_.cols();
619 assert(dim == 3);
620 assert(local_basis_integral.cols() == 9);
621
622 // K_cst = ∫ψ_k
623 // K_lin = ∫∇x(ψ_k), ∫∇y(ψ_k), ∫∇z(ψ_k)
624 // K_mix = ∫(y·∇x(ψ_k)+x·∇y(ψ_k)), ∫(z·∇y(ψ_k)+y·∇z(ψ_k)), ∫(x·∇z(ψ_k)+z·∇x(ψ_k))
625 // K_sqr = ∫x·∇x(ψ_k), ∫y·∇y(ψ_k), ∫z·∇z(ψ_k)
626 Eigen::VectorXd K_cst = Eigen::VectorXd::Zero(num_kernels);
627 Eigen::MatrixXd K_lin = Eigen::MatrixXd::Zero(num_kernels, dim);
628 Eigen::MatrixXd K_mix = Eigen::MatrixXd::Zero(num_kernels, dim);
629 Eigen::MatrixXd K_sqr = Eigen::MatrixXd::Zero(num_kernels, dim);
630 for (int j = 0; j < num_kernels; ++j)
631 {
632 // ∫∇x(ψ_k)(p) = Σ_q (xq - xk) * 1/r * h'(r) * wq
633 // - xq is the x coordinate of the q-th quadrature point
634 // - wq is the q-th quadrature weight
635 // - r is the distance from pq to the kernel center
636 // - h is the RBF kernel (scalar function)
637 for (int q = 0; q < quadr.points.rows(); ++q)
638 {
639 const RowVectorNd p = quadr.points.row(q) - centers_.row(j);
640 const double r = p.norm();
641 const RowVectorNd gradPhi = p * kernel_prime(is_volume(), r) / r * quadr.weights(q);
642 K_cst(j) += kernel(is_volume(), r) * quadr.weights(q);
643 K_lin.row(j) += gradPhi;
644 for (int d = 0; d < dim; ++d)
645 {
646 K_mix(j, d) += quadr.points(q, (d + 1) % dim) * gradPhi(d) + quadr.points(q, d) * gradPhi((d + 1) % dim);
647 }
648 K_sqr.row(j) += (quadr.points.row(q).array() * gradPhi.array()).matrix();
649 }
650 }
651
652 // I_lin = ∫x, ∫y, ∫z
653 // I_sqr = ∫x², ∫y², ∫z²
654 // I_mix = ∫xy, ∫yz, ∫zx
655 Eigen::RowVectorXd I_lin = (quadr.points.array().colwise() * quadr.weights.array()).colwise().sum();
656 Eigen::RowVectorXd I_sqr = (quadr.points.array().square().colwise() * quadr.weights.array()).colwise().sum();
657 Eigen::RowVectorXd I_mix(3);
658 I_mix(0) = (quadr.points.col(0).array() * quadr.points.col(1).array() * quadr.weights.array()).sum();
659 I_mix(1) = (quadr.points.col(1).array() * quadr.points.col(2).array() * quadr.weights.array()).sum();
660 I_mix(2) = (quadr.points.col(2).array() * quadr.points.col(0).array() * quadr.weights.array()).sum();
661 double volume = quadr.weights.sum();
662
663 // Compute M
664 Eigen::Matrix<double, 9, 9> M;
665 M << volume, 0, 0, I_lin(1), 0, I_lin(2), 2 * I_lin(0), 0, 0,
666 0, volume, 0, I_lin(0), I_lin(2), 0, 0, 2 * I_lin(1), 0,
667 0, 0, volume, 0, I_lin(1), I_lin(0), 0, 0, 2 * I_lin(2),
668 I_lin(1), I_lin(0), 0, I_sqr(0) + I_sqr(1), I_mix(2), I_mix(1), 2 * I_mix(0), 2 * I_mix(0), 0,
669 0, I_lin(2), I_lin(1), I_mix(2), I_sqr(1) + I_sqr(2), I_mix(0), 0, 2 * I_mix(1), 2 * I_mix(1),
670 I_lin(2), 0, I_lin(0), I_mix(1), I_mix(0), I_sqr(2) + I_sqr(0), 2 * I_mix(2), 0, 2 * I_mix(2),
671 2 * I_lin(0), 0, 0, 2 * I_mix(0), 0, 2 * I_mix(2), 4 * I_sqr(0), 0, 0,
672 0, 2 * I_lin(1), 0, 2 * I_mix(0), 2 * I_mix(1), 0, 0, 4 * I_sqr(1), 0,
673 0, 0, 2 * I_lin(2), 0, 2 * I_mix(1), 2 * I_mix(2), 0, 0, 4 * I_sqr(2);
674 Eigen::Matrix<double, 1, 9> M_rhs;
675 M_rhs.segment<3>(0) = I_lin;
676 M_rhs.segment<3>(3) = I_mix;
677 M_rhs.segment<3>(6) = I_sqr;
678 // M_rhs << I_lin, I_mix, I_sqr;
679 M.bottomRows(dim).rowwise() += 2.0 * M_rhs;
680 Eigen::FullPivLU<Eigen::Matrix<double, 9, 9>> lu(M);
681 assert(lu.isInvertible());
682
683 // show_matrix_stats(M);
684
685 // Compute L
686 L.resize(num_kernels + 1 + dim + dim * (dim + 1) / 2, num_kernels + 1);
687 L.setZero();
688 L.diagonal().setOnes();
689 L.block(num_kernels + 1, 0, dim, num_kernels) = -K_lin.transpose();
690 L.block(num_kernels + 1 + dim, 0, dim, num_kernels) = -K_mix.transpose();
691 L.block(num_kernels + 1 + dim + dim, 0, dim, num_kernels) = -2.0 * (K_sqr.colwise() + K_cst).transpose();
692 L.bottomRightCorner(dim, 1).setConstant(-2.0 * volume);
693 L.block(num_kernels + 1, 0, 9, num_kernels + 1) = lu.solve(L.block(num_kernels + 1, 0, 9, num_kernels + 1));
694
695 // Compute t
696 t.resize(L.rows(), num_bases);
697 t.setZero();
698 t.bottomRows(9) = local_basis_integral.transpose();
699 t.bottomRows(9) = lu.solve(weights_.bottomRows(9));
700}
701
702// -----------------------------------------------------------------------------
703
704void RBFWithQuadratic::compute_weights(const LinearAssembler &assembler, const Eigen::MatrixXd &samples,
705 const Eigen::MatrixXd &local_basis_integral, const Quadrature &quadr,
706 Eigen::MatrixXd &rhs, bool with_constraints)
707{
708#ifdef VERBOSE
709 logger().trace("#kernel centers: {}", centers_.rows());
710 logger().trace("#collocation points: {}", samples.rows());
711 logger().trace("#quadrature points: {}", quadr.weights.size());
712 logger().trace("#non-vanishing bases: {}", rhs.cols());
713#endif
714
715 if (!with_constraints)
716 {
717 // Compute A
718 Eigen::MatrixXd A;
719 compute_kernels_matrix(samples, A);
720
721 // Solve the system
722 const int num_kernels = centers_.rows();
723 logger().trace("-- Solving system of size {}x{}", num_kernels, num_kernels);
724 weights_ = (A.transpose() * A).ldlt().solve(A.transpose() * rhs);
725 logger().trace("-- Solved!");
726
727 return;
728 }
729
730 const int num_bases = rhs.cols();
731
732 // Compute A
733 Eigen::MatrixXd A;
734 compute_kernels_matrix(samples, A);
735
736 // Compute L and t
737 // Note that t is stored into `weights_` for memory efficiency reasons
738 Eigen::MatrixXd L;
739 if (is_volume())
740 {
741 compute_constraints_matrix_3d(assembler, num_bases, quadr, local_basis_integral, L, weights_);
742 }
743 else
744 {
745 compute_constraints_matrix_2d(assembler, num_bases, quadr, local_basis_integral, L, weights_);
746 }
747
748 // Compute b = rhs - A t
749 Eigen::MatrixXd b = rhs - A * weights_;
750
751// Solve the system
752#ifdef VERBOSE
753 logger().trace("-- Solving system of size {}x{}", L.cols(), L.cols());
754#endif
755 auto ldlt = (L.transpose() * A.transpose() * A * L).ldlt();
756 if (ldlt.info() == Eigen::NumericalIssue)
757 {
758 logger().error("-- WARNING: Numerical issues when solving the harmonic least square.");
759 }
760 weights_ += L * ldlt.solve(L.transpose() * A.transpose() * b);
761#ifdef VERBOSE
762 logger().trace("-- Solved!");
763#endif
764
765#ifdef VERBOSE
766 logger().trace("-- Mean residual: {}", (A * weights_ - rhs).array().abs().colwise().maxCoeff().mean());
767#endif
768}
double val
Definition Assembler.cpp:86
QuadratureVector da
Definition Assembler.cpp:23
ElementAssemblyValues vals
Definition Assembler.cpp:22
int x
virtual bool is_tensor() const
virtual VectorNd compute_rhs(const AutodiffHessianPt &pt) const
stores per element basis values at given quadrature points and geometric mapping
assemble matrix based on the local assembler local assembler is eg Laplace, LinearElasticity etc
void assemble(const bool is_volume, const int n_basis, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const AssemblyValsCache &cache, const double t, StiffnessMatrix &stiffness, const bool is_mass=false) const override
assembles the stiffness matrix for the given basis the bilinear form (local assembler) is encoded by ...
static void setup_monomials_strong_2d(const int dim, const assembler::LinearAssembler &assembler, const Eigen::MatrixXd &pts, const QuadratureVector &da, std::array< Eigen::MatrixXd, 5 > &strong)
void basis(const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) const
Evaluates one RBF function over a list of coordinates.
void compute_constraints_matrix_3d(const assembler::LinearAssembler &assembler, const int num_bases, const quadrature::Quadrature &quadr, const Eigen::MatrixXd &local_basis_integral, Eigen::MatrixXd &L, Eigen::MatrixXd &t) const
void bases_values(const Eigen::MatrixXd &samples, Eigen::MatrixXd &val) const
Batch evaluates the RBF + polynomials on a set of sample points.
void bases_grads(const int axis, const Eigen::MatrixXd &samples, Eigen::MatrixXd &val) const
Batch evaluates the gradient of the RBF + polynomials on a set of sample points.
void compute_constraints_matrix_2d_old(const int num_bases, const quadrature::Quadrature &quadr, const Eigen::MatrixXd &local_basis_integral, Eigen::MatrixXd &L, Eigen::MatrixXd &t) const
void compute_weights(const assembler::LinearAssembler &assembler, const Eigen::MatrixXd &collocation_points, const Eigen::MatrixXd &local_basis_integral, const quadrature::Quadrature &quadr, Eigen::MatrixXd &rhs, bool with_constraints)
static void setup_monomials_vals_2d(const int star_index, const Eigen::MatrixXd &pts, assembler::ElementAssemblyValues &vals)
void compute_constraints_matrix_2d(const assembler::LinearAssembler &assembler, const int num_bases, const quadrature::Quadrature &quadr, const Eigen::MatrixXd &local_basis_integral, Eigen::MatrixXd &L, Eigen::MatrixXd &t) const
RBFWithQuadratic(const assembler::LinearAssembler &assembler, const Eigen::MatrixXd &centers, const Eigen::MatrixXd &collocation_points, const Eigen::MatrixXd &local_basis_integral, const quadrature::Quadrature &quadr, Eigen::MatrixXd &rhs, bool with_constraints=true)
Initialize RBF functions over a polytope element.
void grad(const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) const
Evaluates the gradient of one RBF function over a list of coordinates.
void compute_kernels_matrix(const Eigen::MatrixXd &samples, Eigen::MatrixXd &A) const
Used for test only.
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44
Eigen::Matrix< AutodiffScalarHessian, Eigen::Dynamic, 1, 0, 3, 1 > AutodiffHessianPt
DScalar2< double, Eigen::Matrix< double, Eigen::Dynamic, 1, 0, 3, 1 >, Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > > AutodiffScalarHessian
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, MAX_QUAD_POINTS, 1 > QuadratureVector
Definition Types.hpp:17
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:13
static void setVariableCount(size_t value)
Set the independent variable count used by the automatic differentiation layer.
Definition autodiff.h:54