144 vals.
basis_values[star_index + 0].grad = Eigen::MatrixXd(pts.rows(), pts.cols());
150 vals.
basis_values[star_index + 1].grad = Eigen::MatrixXd(pts.rows(), pts.cols());
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());
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());
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());
172 for (
size_t i = star_index; i < star_index + 5; ++i)
386 const int num_kernels =
centers_.rows();
389 A.resize(samples.rows(), num_kernels + 1 + dim + dim * (dim + 1) / 2);
390 for (
int j = 0; j < num_kernels; ++j)
392 A.col(j) = (samples.rowwise() -
centers_.row(j)).rowwise().norm().unaryExpr([
this](
double x) {
return kernel(
is_volume(),
x); });
394 A.col(num_kernels).setOnes();
395 A.middleCols(num_kernels + 1, dim) = samples;
398 A.middleCols(num_kernels + dim + 1, 1) = samples.rowwise().prod();
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();
407 A.rightCols(dim) = samples.array().square();
415 const Eigen::MatrixXd &local_basis_integral,
417 Eigen::MatrixXd &t)
const
419 const int num_kernels =
centers_.rows();
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)
438 for (
int q = 0; q < quadr.
points.rows(); ++q)
441 const double r = p.norm();
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();
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();
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());
472 L.resize(num_kernels + 1 + dim + dim * (dim + 1) / 2, num_kernels + 1);
474 L.diagonal().setOnes();
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);
489 L.block(num_kernels + 1, 0, 5, num_kernels + 1) = lu.solve(L.block(num_kernels + 1, 0, 5, num_kernels + 1));
492 t.resize(L.rows(), num_bases);
494 t.bottomRows(5) = local_basis_integral.transpose();
495 t.bottomRows(5) = lu.solve(
weights_.bottomRows(5));
502 const Eigen::MatrixXd &local_basis_integral,
504 Eigen::MatrixXd &t)
const
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);
513 std::array<Eigen::MatrixXd, 5> strong;
525 for (
int j = 0; j < num_kernels; ++j)
530 for (
int q = 0; q < quadr.
points.rows(); ++q)
533 const double r = p.norm();
540 for (
size_t i = 5; i < ass_val.
basis_values.size(); ++i)
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)
548 for (
int j = 0; j < 5; ++j)
552 for (
int d1 = 0; d1 < assembler_dim; ++d1)
554 for (
int d2 = 0; d2 < assembler_dim; ++d2)
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();
563 Eigen::FullPivLU<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, 10, 10>> lu(M);
564 assert(lu.isInvertible());
567 L.resize((num_kernels + 1 + space_dim + space_dim * (space_dim + 1) / 2) * assembler_dim, (num_kernels + 1) * assembler_dim);
569 L.diagonal().setOnes();
571 for (
int i = 0; i < 5; ++i)
573 for (
int j = 0; j < num_kernels; ++j)
576 for (
int d1 = 0; d1 < assembler_dim; ++d1)
578 for (
int d2 = 0; d2 < assembler_dim; ++d2)
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();
586 for (
int d1 = 0; d1 < assembler_dim; ++d1)
588 for (
int d2 = 0; d2 < assembler_dim; ++d2)
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();
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));
602 t.resize(L.rows(), num_bases * assembler_dim);
604 t.bottomRows(5 * assembler_dim) = lu.solve(local_basis_integral.transpose());
613 const Eigen::MatrixXd &local_basis_integral,
615 Eigen::MatrixXd &t)
const
617 const int num_kernels =
centers_.rows();
620 assert(local_basis_integral.cols() == 9);
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)
637 for (
int q = 0; q < quadr.
points.rows(); ++q)
640 const double r = p.norm();
643 K_lin.row(j) += gradPhi;
644 for (
int d = 0; d < dim; ++d)
646 K_mix(j, d) += quadr.
points(q, (d + 1) % dim) * gradPhi(d) + quadr.
points(q, d) * gradPhi((d + 1) % dim);
648 K_sqr.row(j) += (quadr.
points.row(q).array() * gradPhi.array()).matrix();
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();
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;
679 M.bottomRows(dim).rowwise() += 2.0 * M_rhs;
680 Eigen::FullPivLU<Eigen::Matrix<double, 9, 9>> lu(M);
681 assert(lu.isInvertible());
686 L.resize(num_kernels + 1 + dim + dim * (dim + 1) / 2, num_kernels + 1);
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));
696 t.resize(L.rows(), num_bases);
698 t.bottomRows(9) = local_basis_integral.transpose();
699 t.bottomRows(9) = lu.solve(
weights_.bottomRows(9));
705 const Eigen::MatrixXd &local_basis_integral,
const Quadrature &quadr,
706 Eigen::MatrixXd &rhs,
bool with_constraints)
710 logger().trace(
"#collocation points: {}", samples.rows());
711 logger().trace(
"#quadrature points: {}", quadr.
weights.size());
712 logger().trace(
"#non-vanishing bases: {}", rhs.cols());
715 if (!with_constraints)
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!");
730 const int num_bases = rhs.cols();
749 Eigen::MatrixXd b = rhs - A *
weights_;
753 logger().trace(
"-- Solving system of size {}x{}", L.cols(), L.cols());
755 auto ldlt = (L.transpose() * A.transpose() * A * L).ldlt();
756 if (ldlt.info() == Eigen::NumericalIssue)
758 logger().error(
"-- WARNING: Numerical issues when solving the harmonic least square.");
760 weights_ += L * ldlt.solve(L.transpose() * A.transpose() * b);
762 logger().trace(
"-- Solved!");
766 logger().trace(
"-- Mean residual: {}", (A *
weights_ - rhs).array().abs().colwise().maxCoeff().mean());