132 const int num_kernels =
centers_.rows();
136 Eigen::MatrixXd A_prime(samples.rows(), num_kernels + 1 + dim + dim * (dim + 1) / 2);
139 for (
int j = 0; j < num_kernels; ++j)
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();
145 A_prime.middleCols(num_kernels + 1 + axis, 1).setOnes();
149 A_prime.col(num_kernels + 1 + dim) = samples.col(1 - axis);
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);
157 A_prime.rightCols(dim).col(axis) = 2.0 * samples.col(axis);
168 const int num_kernels =
centers_.rows();
171 A.resize(samples.rows(), num_kernels + 1 + dim + dim * (dim + 1) / 2);
172 for (
int j = 0; j < num_kernels; ++j)
174 A.col(j) = (samples.rowwise() -
centers_.row(j)).rowwise().norm().unaryExpr([
this](
double x) {
return kernel(
is_volume(),
x); });
176 A.col(num_kernels).setOnes();
177 A.middleCols(num_kernels + 1, dim) = samples;
180 A.middleCols(num_kernels + dim + 1, 1) = samples.rowwise().prod();
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();
189 A.rightCols(dim) = samples.array().square();
195 const int num_bases,
const Quadrature &quadr, Eigen::MatrixXd &C)
const
197 const int num_kernels =
centers_.rows();
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)
216 for (
int q = 0; q < quadr.
points.rows(); ++q)
219 const double r = p.norm();
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();
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();
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);
249 C.resize(5, num_kernels + 1 + dim + dim * (dim + 1) / 2);
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;
259 const int num_bases,
const Quadrature &quadr, Eigen::MatrixXd &C)
const
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);
268 std::array<Eigen::MatrixXd, 5> strong;
280 for (
int j = 0; j < num_kernels; ++j)
285 for (
int q = 0; q < quadr.
points.rows(); ++q)
288 const double r = p.norm();
295 for (
size_t i = 5; i < ass_val.
basis_values.size(); ++i)
304 for (
int d = 0; d < 5; ++d)
307 for (
int i = 0; i < num_kernels; ++i)
310 for (
int alpha = 0; alpha < assembler_dim; ++alpha)
312 for (
int beta = 0; beta < assembler_dim; ++beta)
314 const int loc_index = alpha * assembler_dim + beta;
321 for (
int i = 0; i < 5; ++i)
324 for (
int alpha = 0; alpha < assembler_dim; ++alpha)
326 for (
int beta = 0; beta < assembler_dim; ++beta)
328 const int loc_index = alpha * assembler_dim + beta;
335 for (
int alpha = 0; alpha < assembler_dim; ++alpha)
337 for (
int beta = 0; beta < assembler_dim; ++beta)
339 const int loc_index = alpha * assembler_dim + beta;
365 const int num_bases,
const Quadrature &quadr, Eigen::MatrixXd &C)
const
367 const int num_kernels =
centers_.rows();
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)
386 for (
int q = 0; q < quadr.
points.rows(); ++q)
389 const double r = p.norm();
392 K_lin.row(j) += gradPhi;
393 for (
int d = 0; d < dim; ++d)
395 K_mix(j, d) += quadr.
points(q, (d + 1) % dim) * gradPhi(d) + quadr.
points(q, d) * gradPhi((d + 1) % dim);
397 K_sqr.row(j) += (quadr.
points.row(q).array() * gradPhi.array()).matrix();
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();
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;
428 M.bottomRows(dim).rowwise() += 2.0 * M_rhs;
435 C.resize(9, num_kernels + 1 + dim + dim * (dim + 1) / 2);
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;
447 const Eigen::MatrixXd &local_basis_integral,
const Quadrature &quadr,
448 Eigen::MatrixXd &b,
bool with_constraints)
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());
460 if (!with_constraints)
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!");
471 const int num_bases = b.cols();
473 const Eigen::MatrixXd At = A.transpose();
487 assert(
centers_.rows() + 1 + dim + dim * (dim + 1) / 2 > C.rows());
488 logger().trace(
"#constraints: {}", C.rows());
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();
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();
508 logger().trace(
"-- Solving system of size {}x{}", M.rows(), M.cols());
509 auto ldlt = M.ldlt();
510 if (ldlt.info() == Eigen::NumericalIssue)
512 logger().error(
"-- WARNING: Numerical issues when solving the harmonic least square.");
514 const auto tmp = ldlt.solve(rhs);
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());