Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
RBFWithQuadraticLagrange.cpp
Go to the documentation of this file.
1
4
8
9#include <igl/Timer.h>
10#include <Eigen/Dense>
11#include <iostream>
12#include <fstream>
13#include <array>
15
16using namespace polyfem;
17using namespace polyfem::assembler;
18using namespace polyfem::basis;
19using namespace polyfem::quadrature;
20using namespace polyfem::utils;
21
22namespace
23{
24
25 // Harmonic kernel
26 double kernel(const bool is_volume, const double r)
27 {
28 if (r < 1e-8)
29 {
30 return 0;
31 }
32
33 if (is_volume)
34 {
35 return 1 / r;
36 }
37 else
38 {
39 return log(r);
40 }
41 }
42
43 double kernel_prime(const bool is_volume, const double r)
44 {
45 if (r < 1e-8)
46 {
47 return 0;
48 }
49
50 if (is_volume)
51 {
52 return -1 / (r * r);
53 }
54 else
55 {
56 return 1 / r;
57 }
58 }
59
60 // Biharmonic kernel (2d only)
61 // double kernel(const bool is_volume, const double r) {
62 // assert(!is_volume);
63 // if (r < 1e-8) { return 0; }
64
65 // return r * r * (log(r)-1);
66 // }
67
68 // double kernel_prime(const bool is_volume, const double r) {
69 // assert(!is_volume);
70 // if (r < 1e-8) { return 0; }
71
72 // return r * ( 2 * log(r) - 1);
73 // }
74
75} // anonymous namespace
76
78
80 const LinearAssembler &assembler,
81 const Eigen::MatrixXd &centers,
82 const Eigen::MatrixXd &collocation_points,
83 const Eigen::MatrixXd &local_basis_integral,
84 const Quadrature &quadr,
85 Eigen::MatrixXd &rhs,
86 bool with_constraints)
87 : centers_(centers)
88{
89 // centers_.resize(0, centers.cols());
90 compute_weights(assembler, collocation_points, local_basis_integral, quadr, rhs, with_constraints);
91}
92
93// -----------------------------------------------------------------------------
94
95void RBFWithQuadraticLagrange::basis(const int local_index, const Eigen::MatrixXd &samples, Eigen::MatrixXd &val) const
96{
97 Eigen::MatrixXd tmp;
98 bases_values(samples, tmp);
99 val = tmp.col(local_index);
100}
101
102// -----------------------------------------------------------------------------
103
104void RBFWithQuadraticLagrange::grad(const int local_index, const Eigen::MatrixXd &samples, Eigen::MatrixXd &val) const
105{
106 Eigen::MatrixXd tmp;
107 const int dim = centers_.cols();
108 val.resize(samples.rows(), dim);
109 for (int d = 0; d < dim; ++d)
110 {
111 bases_grads(d, samples, tmp);
112 val.col(d) = tmp.col(local_index);
113 }
114}
115
117
118void RBFWithQuadraticLagrange::bases_values(const Eigen::MatrixXd &samples, Eigen::MatrixXd &val) const
119{
120 // Compute A
121 Eigen::MatrixXd A;
122 compute_kernels_matrix(samples, A);
123
124 // Multiply by the weights
125 val = A * weights_;
126}
127
128// -----------------------------------------------------------------------------
129
130void RBFWithQuadraticLagrange::bases_grads(const int axis, const Eigen::MatrixXd &samples, Eigen::MatrixXd &val) const
131{
132 const int num_kernels = centers_.rows();
133 const int dim = (is_volume() ? 3 : 2);
134
135 // Compute ∇xA
136 Eigen::MatrixXd A_prime(samples.rows(), num_kernels + 1 + dim + dim * (dim + 1) / 2);
137 A_prime.setZero();
138
139 for (int j = 0; j < num_kernels; ++j)
140 {
141 A_prime.col(j) = (samples.rowwise() - centers_.row(j)).rowwise().norm().unaryExpr([this](double x) { return kernel_prime(is_volume(), x) / x; });
142 A_prime.col(j) = (samples.col(axis).array() - centers_(j, axis)) * A_prime.col(j).array();
143 }
144 // Linear terms
145 A_prime.middleCols(num_kernels + 1 + axis, 1).setOnes();
146 // Mixed terms
147 if (dim == 2)
148 {
149 A_prime.col(num_kernels + 1 + dim) = samples.col(1 - axis);
150 }
151 else
152 {
153 A_prime.col(num_kernels + 1 + dim + axis) = samples.col((axis + 1) % dim);
154 A_prime.col(num_kernels + 1 + dim + (axis + 2) % dim) = samples.col((axis + 2) % dim);
155 }
156 // Quadratic terms
157 A_prime.rightCols(dim).col(axis) = 2.0 * samples.col(axis);
158
159 // Apply weights
160 val = A_prime * weights_;
161}
162
164
165void RBFWithQuadraticLagrange::compute_kernels_matrix(const Eigen::MatrixXd &samples, Eigen::MatrixXd &A) const
166{
167 // Compute A
168 const int num_kernels = centers_.rows();
169 const int dim = (is_volume() ? 3 : 2);
170
171 A.resize(samples.rows(), num_kernels + 1 + dim + dim * (dim + 1) / 2);
172 for (int j = 0; j < num_kernels; ++j)
173 {
174 A.col(j) = (samples.rowwise() - centers_.row(j)).rowwise().norm().unaryExpr([this](double x) { return kernel(is_volume(), x); });
175 }
176 A.col(num_kernels).setOnes(); // constant term
177 A.middleCols(num_kernels + 1, dim) = samples; // linear terms
178 if (dim == 2)
179 {
180 A.middleCols(num_kernels + dim + 1, 1) = samples.rowwise().prod(); // mixed terms
181 }
182 else if (dim == 3)
183 {
184 A.middleCols(num_kernels + dim + 1, 3) = samples;
185 A.middleCols(num_kernels + dim + 1 + 0, 1).array() *= samples.col(1).array();
186 A.middleCols(num_kernels + dim + 1 + 1, 1).array() *= samples.col(2).array();
187 A.middleCols(num_kernels + dim + 1 + 2, 1).array() *= samples.col(0).array();
188 }
189 A.rightCols(dim) = samples.array().square(); // quadratic terms
190}
191
192// -----------------------------------------------------------------------------
193
195 const int num_bases, const Quadrature &quadr, Eigen::MatrixXd &C) const
196{
197 const int num_kernels = centers_.rows();
198 const int dim = centers_.cols();
199 assert(dim == 2);
200
201 // K_cst = ∫φj
202 // K_lin = ∫∇x(φj), ∫∇y(φj)
203 // K_mix = ∫y·∇x(φj), ∫x·∇y(φj)
204 // K_sqr = ∫x·∇x(φj), ∫y·∇y(φj)
205 Eigen::VectorXd K_cst = Eigen::VectorXd::Zero(num_kernels);
206 Eigen::MatrixXd K_lin = Eigen::MatrixXd::Zero(num_kernels, dim);
207 Eigen::MatrixXd K_mix = Eigen::MatrixXd::Zero(num_kernels, dim);
208 Eigen::MatrixXd K_sqr = Eigen::MatrixXd::Zero(num_kernels, dim);
209 for (int j = 0; j < num_kernels; ++j)
210 {
211 // ∫∇x(φj)(p) = Σ_q (xq - xk) * 1/r * h'(r) * wq
212 // - xq is the x coordinate of the q-th quadrature point
213 // - wq is the q-th quadrature weight
214 // - r is the distance from pq to the kernel center
215 // - h is the RBF kernel (scalar function)
216 for (int q = 0; q < quadr.points.rows(); ++q)
217 {
218 const RowVectorNd p = quadr.points.row(q) - centers_.row(j);
219 const double r = p.norm();
220 const RowVectorNd gradPhi = p * kernel_prime(is_volume(), r) / r * quadr.weights(q);
221 K_cst(j) += kernel(is_volume(), r) * quadr.weights(q);
222 K_lin.row(j) += gradPhi;
223 K_mix(j, 0) += quadr.points(q, 1) * gradPhi(0);
224 K_mix(j, 1) += quadr.points(q, 0) * gradPhi(1);
225 K_sqr.row(j) += (quadr.points.row(q).array() * gradPhi.array()).matrix();
226 }
227 }
228
229 // I_lin = ∫x, ∫y
230 // I_mix = ∫xy
231 // I_sqr = ∫x², ∫y²
232 Eigen::RowVectorXd I_lin = (quadr.points.array().colwise() * quadr.weights.array()).colwise().sum();
233 Eigen::RowVectorXd I_mix = (quadr.points.rowwise().prod().array() * quadr.weights.array()).colwise().sum();
234 Eigen::RowVectorXd I_sqr = (quadr.points.array().square().colwise() * quadr.weights.array()).colwise().sum();
235 double volume = quadr.weights.sum();
236
237 // Compute M
238 Eigen::Matrix<double, 5, 5> M;
239 M << volume, 0, I_lin(1), 2 * I_lin(0), 0,
240 0, volume, I_lin(0), 0, 2 * I_lin(1),
241 I_lin(1), I_lin(0), I_sqr(0) + I_sqr(1), 2 * I_mix(0), 2 * I_mix(0),
242 4 * I_lin(0), 2 * I_lin(1), 4 * I_mix(0), 6 * I_sqr(0), 2 * I_sqr(1),
243 2 * I_lin(0), 4 * I_lin(1), 4 * I_mix(0), 2 * I_sqr(0), 6 * I_sqr(1);
244 Eigen::FullPivLU<Eigen::Matrix<double, 5, 5>> lu(M);
245
247
248 // Compute L
249 C.resize(5, num_kernels + 1 + dim + dim * (dim + 1) / 2);
250 C.setZero();
251 C.block(0, 0, dim, num_kernels) = K_lin.transpose();
252 C.block(dim, 0, 1, num_kernels) = K_mix.transpose().colwise().sum();
253 C.block(dim + 1, 0, dim, num_kernels) = 2.0 * (K_sqr.colwise() + K_cst).transpose();
254 C.block(dim + 1, num_kernels, dim, 1).setConstant(2.0 * volume);
255 C.bottomRightCorner(5, 5) = M;
256}
257
259 const int num_bases, const Quadrature &quadr, Eigen::MatrixXd &C) const
260{
261 // TODO
262 const double t = 0;
263 const int num_kernels = centers_.rows();
264 const int space_dim = centers_.cols();
265 const int assembler_dim = assembler.is_tensor() ? 2 : 1;
266 assert(space_dim == 2);
267
268 std::array<Eigen::MatrixXd, 5> strong;
269
270 // ass_val = [q_10, q_01, q_11, q_20, q_02, psi_0, ..., psi_k]
271 ElementAssemblyValues ass_val;
272 ass_val.has_parameterization = false;
273 ass_val.basis_values.resize(5 + num_kernels);
274
275 // evaluating monomial and grad of monomials at quad points
277 RBFWithQuadratic::setup_monomials_strong_2d(assembler_dim, assembler, quadr.points, quadr.weights.array(), strong);
278
279 // evaluating psi and grad psi at quadr points
280 for (int j = 0; j < num_kernels; ++j)
281 {
282 ass_val.basis_values[5 + j].val = Eigen::MatrixXd(quadr.points.rows(), 1);
283 ass_val.basis_values[5 + j].grad = Eigen::MatrixXd(quadr.points.rows(), quadr.points.cols());
284
285 for (int q = 0; q < quadr.points.rows(); ++q)
286 {
287 const RowVectorNd p = quadr.points.row(q) - centers_.row(j);
288 const double r = p.norm();
289
290 ass_val.basis_values[5 + j].val(q) = kernel(is_volume(), r);
291 ass_val.basis_values[5 + j].grad.row(q) = p * kernel_prime(is_volume(), r) / r;
292 }
293 }
294
295 for (size_t i = 5; i < ass_val.basis_values.size(); ++i)
296 {
297 ass_val.basis_values[i].grad_t_m = ass_val.basis_values[i].grad;
298 }
299
300 // Compute C
301 C.resize(RBFWithQuadratic::index_mapping(assembler_dim - 1, assembler_dim - 1, 4, assembler_dim) + 1, num_kernels + 1 + 5);
302 C.setZero();
303
304 for (int d = 0; d < 5; ++d)
305 {
306 // first num_kernels bases
307 for (int i = 0; i < num_kernels; ++i)
308 {
309 const auto tmp = assembler.assemble(LinearAssemblerData(ass_val, t, d, 5 + i, quadr.weights));
310 for (int alpha = 0; alpha < assembler_dim; ++alpha)
311 {
312 for (int beta = 0; beta < assembler_dim; ++beta)
313 {
314 const int loc_index = alpha * assembler_dim + beta;
315 C(RBFWithQuadratic::index_mapping(alpha, beta, d, assembler_dim), i) = tmp(loc_index) + (strong[d].row(loc_index).transpose().array() * ass_val.basis_values[5 + i].val.array()).sum();
316 }
317 }
318 }
319
320 // second the q_i
321 for (int i = 0; i < 5; ++i)
322 {
323 const auto tmp = assembler.assemble(LinearAssemblerData(ass_val, t, d, i, quadr.weights));
324 for (int alpha = 0; alpha < assembler_dim; ++alpha)
325 {
326 for (int beta = 0; beta < assembler_dim; ++beta)
327 {
328 const int loc_index = alpha * assembler_dim + beta;
329 C(RBFWithQuadratic::index_mapping(alpha, beta, d, assembler_dim), num_kernels + i + 1) = tmp(loc_index) + (strong[d].row(loc_index).transpose().array() * ass_val.basis_values[i].val.array()).sum();
330 }
331 }
332 }
333
334 // finally the constant
335 for (int alpha = 0; alpha < assembler_dim; ++alpha)
336 {
337 for (int beta = 0; beta < assembler_dim; ++beta)
338 {
339 const int loc_index = alpha * assembler_dim + beta;
340 C(RBFWithQuadratic::index_mapping(alpha, beta, d, assembler_dim), num_kernels) = strong[d].row(loc_index).sum();
341 }
342 }
343 }
344
345 // {
346 // std::ofstream file;
347 // file.open("C.txt");
348 // file << C;
349 // file.close();
350 // }
351
352 // Eigen::MatrixXd Cold;
353 // compute_constraints_matrix_2d_old(num_bases, quadr,Cold);
354 // {
355 // std::ofstream file;
356 // file.open("Cold.txt");
357 // file << Cold;
358 // file.close();
359 // }
360}
361
362// -----------------------------------------------------------------------------
363
365 const int num_bases, const Quadrature &quadr, Eigen::MatrixXd &C) const
366{
367 const int num_kernels = centers_.rows();
368 const int dim = centers_.cols();
369 assert(dim == 3);
370
371 // K_cst = ∫φj
372 // K_lin = ∫∇x(φj), ∫∇y(φj), ∫∇z(φj)
373 // K_mix = ∫(y·∇x(φj)+x·∇y(φj)), ∫(z·∇y(φj)+y·∇z(φj)), ∫(x·∇z(φj)+z·∇x(φj))
374 // K_sqr = ∫x·∇x(φj), ∫y·∇y(φj), ∫z·∇z(φj)
375 Eigen::VectorXd K_cst = Eigen::VectorXd::Zero(num_kernels);
376 Eigen::MatrixXd K_lin = Eigen::MatrixXd::Zero(num_kernels, dim);
377 Eigen::MatrixXd K_mix = Eigen::MatrixXd::Zero(num_kernels, dim);
378 Eigen::MatrixXd K_sqr = Eigen::MatrixXd::Zero(num_kernels, dim);
379 for (int j = 0; j < num_kernels; ++j)
380 {
381 // ∫∇x(φj)(p) = Σ_q (xq - xk) * 1/r * h'(r) * wq
382 // - xq is the x coordinate of the q-th quadrature point
383 // - wq is the q-th quadrature weight
384 // - r is the distance from pq to the kernel center
385 // - h is the RBF kernel (scalar function)
386 for (int q = 0; q < quadr.points.rows(); ++q)
387 {
388 const RowVectorNd p = quadr.points.row(q) - centers_.row(j);
389 const double r = p.norm();
390 const RowVectorNd gradPhi = p * kernel_prime(is_volume(), r) / r * quadr.weights(q);
391 K_cst(j) += kernel(is_volume(), r) * quadr.weights(q);
392 K_lin.row(j) += gradPhi;
393 for (int d = 0; d < dim; ++d)
394 {
395 K_mix(j, d) += quadr.points(q, (d + 1) % dim) * gradPhi(d) + quadr.points(q, d) * gradPhi((d + 1) % dim);
396 }
397 K_sqr.row(j) += (quadr.points.row(q).array() * gradPhi.array()).matrix();
398 }
399 }
400
401 // I_lin = ∫x, ∫y, ∫z
402 // I_sqr = ∫x², ∫y², ∫z²
403 // I_mix = ∫xy, ∫yz, ∫zx
404 Eigen::RowVectorXd I_lin = (quadr.points.array().colwise() * quadr.weights.array()).colwise().sum();
405 Eigen::RowVectorXd I_sqr = (quadr.points.array().square().colwise() * quadr.weights.array()).colwise().sum();
406 Eigen::RowVectorXd I_mix(3);
407 I_mix(0) = (quadr.points.col(0).array() * quadr.points.col(1).array() * quadr.weights.array()).sum();
408 I_mix(1) = (quadr.points.col(1).array() * quadr.points.col(2).array() * quadr.weights.array()).sum();
409 I_mix(2) = (quadr.points.col(2).array() * quadr.points.col(0).array() * quadr.weights.array()).sum();
410 double volume = quadr.weights.sum();
411
412 // Compute M
413 Eigen::Matrix<double, 9, 9> M;
414 M << volume, 0, 0, I_lin(1), 0, I_lin(2), 2 * I_lin(0), 0, 0,
415 0, volume, 0, I_lin(0), I_lin(2), 0, 0, 2 * I_lin(1), 0,
416 0, 0, volume, 0, I_lin(1), I_lin(0), 0, 0, 2 * I_lin(2),
417 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,
418 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),
419 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),
420 2 * I_lin(0), 0, 0, 2 * I_mix(0), 0, 2 * I_mix(2), 4 * I_sqr(0), 0, 0,
421 0, 2 * I_lin(1), 0, 2 * I_mix(0), 2 * I_mix(1), 0, 0, 4 * I_sqr(1), 0,
422 0, 0, 2 * I_lin(2), 0, 2 * I_mix(1), 2 * I_mix(2), 0, 0, 4 * I_sqr(2);
423 Eigen::Matrix<double, 1, 9> M_rhs;
424 M_rhs.segment<3>(0) = I_lin;
425 M_rhs.segment<3>(3) = I_mix;
426 M_rhs.segment<3>(6) = I_sqr;
427 // M_rhs << I_lin, I_mix, I_sqr;
428 M.bottomRows(dim).rowwise() += 2.0 * M_rhs;
429 // TODO
430 // assert(false);
431
433
434 // Compute L
435 C.resize(9, num_kernels + 1 + dim + dim * (dim + 1) / 2);
436 C.setZero();
437 C.block(0, 0, dim, num_kernels) = K_lin.transpose();
438 C.block(dim, 0, dim, num_kernels) = K_mix.transpose();
439 C.block(dim + dim, 0, dim, num_kernels) = 2.0 * (K_sqr.colwise() + K_cst).transpose();
440 C.block(dim + dim, num_kernels, dim, 1).setConstant(2.0 * volume);
441 C.bottomRightCorner(9, 9) = M;
442}
443
444// -----------------------------------------------------------------------------
445
446void RBFWithQuadraticLagrange::compute_weights(const LinearAssembler &assembler, const Eigen::MatrixXd &samples,
447 const Eigen::MatrixXd &local_basis_integral, const Quadrature &quadr,
448 Eigen::MatrixXd &b, bool with_constraints)
449{
450 logger().trace("#kernel centers: {}", centers_.rows());
451 logger().trace("#collocation points: {}", samples.rows());
452 logger().trace("#quadrature points: {}", quadr.weights.size());
453 logger().trace("#non-vanishing bases: {}", b.cols());
454 logger().trace("#constraints: {}", b.cols());
455
456 // Compute A
457 Eigen::MatrixXd A;
458 compute_kernels_matrix(samples, A);
459
460 if (!with_constraints)
461 {
462 // Solve the system
463 const int num_kernels = centers_.rows();
464 logger().trace("-- Solving system of size {}x{}", num_kernels, num_kernels);
465 weights_ = (A.transpose() * A).ldlt().solve(A.transpose() * b);
466 logger().trace("-- Solved!");
467
468 return;
469 }
470
471 const int num_bases = b.cols();
472
473 const Eigen::MatrixXd At = A.transpose();
474
475 // Compute C
476 Eigen::MatrixXd C;
477 if (is_volume())
478 {
479 compute_constraints_matrix_3d(num_bases, quadr, C);
480 }
481 else
482 {
483 compute_constraints_matrix_2d(assembler, num_bases, quadr, C);
484 }
485
486 const int dim = centers_.cols();
487 assert(centers_.rows() + 1 + dim + dim * (dim + 1) / 2 > C.rows());
488 logger().trace("#constraints: {}", C.rows());
489
490 // Compute rhs = [ A^T b; d ]
491 assert(local_basis_integral.cols() == C.rows());
492 assert(local_basis_integral.rows() == b.cols());
493 assert(A.rows() == b.rows());
494 Eigen::MatrixXd rhs(A.cols() + local_basis_integral.cols(), b.cols());
495 rhs.topRows(A.cols()) = At * b;
496 rhs.bottomRows(local_basis_integral.cols()) = local_basis_integral.transpose();
497
498 // Compute M = [ A^T A, C^T; C, 0]
499 assert(C.cols() == A.cols());
500 assert(A.rows() == b.rows());
501 Eigen::MatrixXd M(A.cols() + C.rows(), A.cols() + C.rows());
502 M.topLeftCorner(A.cols(), A.cols()) = At * A;
503 M.topRightCorner(A.cols(), C.rows()) = C.transpose();
504 M.bottomLeftCorner(C.rows(), A.cols()) = C;
505 M.bottomRightCorner(C.rows(), C.rows()).setZero();
506
507 // Solve the system
508 logger().trace("-- Solving system of size {}x{}", M.rows(), M.cols());
509 auto ldlt = M.ldlt();
510 if (ldlt.info() == Eigen::NumericalIssue)
511 {
512 logger().error("-- WARNING: Numerical issues when solving the harmonic least square.");
513 }
514 const auto tmp = ldlt.solve(rhs);
515 weights_ = tmp.topRows(A.cols());
516 logger().trace("-- Solved!");
517 logger().trace("-- Mean residual: {}", (A * weights_ - b).array().abs().colwise().maxCoeff().mean());
518 logger().trace("-- Max constraints error: {}", (C * weights_ - local_basis_integral.transpose()).array().abs().maxCoeff());
519
520 // {
521 // std::ofstream file;
522 // file.open("M.txt");
523 // file << M;
524 // file.close();
525 // }
526
527 // {
528 // std::ofstream file;
529 // file.open("C.txt");
530 // file << C;
531 // file.close();
532 // }
533
534 // {
535 // std::ofstream file;
536 // file.open("b.txt");
537 // file << local_basis_integral.transpose();
538 // file.close();
539 // }
540
541 // {
542 // std::ofstream file;
543 // file.open("rhs.txt");
544 // file << rhs;
545 // file.close();
546 // }
547}
double val
Definition Assembler.cpp:86
int x
virtual bool is_tensor() 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)
static void setup_monomials_vals_2d(const int star_index, const Eigen::MatrixXd &pts, assembler::ElementAssemblyValues &vals)
static int index_mapping(const int alpha, const int beta, const int d, const int ass_dim)
void compute_constraints_matrix_2d_old(const int num_bases, const quadrature::Quadrature &quadr, Eigen::MatrixXd &C) const
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(const assembler::LinearAssembler &assembler, const int num_bases, const quadrature::Quadrature &quadr, Eigen::MatrixXd &C) const
void bases_values(const Eigen::MatrixXd &samples, Eigen::MatrixXd &val) const
Batch evaluates the RBF + polynomials on a set of sample points.
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)
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_constraints_matrix_3d(const int num_bases, const quadrature::Quadrature &quadr, Eigen::MatrixXd &C) const
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_kernels_matrix(const Eigen::MatrixXd &samples, Eigen::MatrixXd &A) const
RBFWithQuadraticLagrange(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.
Used for test only.
void show_matrix_stats(const Eigen::MatrixXd &M)
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:13