You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
121 lines
3.5 KiB
121 lines
3.5 KiB
#include <medusa/Medusa.hpp>
|
|
#include <medusa/bits/domains/GeneralFill.hpp>
|
|
#include <Eigen/SparseCore>
|
|
#include <Eigen/SparseLU>
|
|
#include <iostream>
|
|
#include <sstream>
|
|
#include "gtest/gtest.h"
|
|
|
|
namespace mm {
|
|
constexpr int dim = 2; // Change to 1 or 3.
|
|
|
|
|
|
TEST(Approximations, printOperators) {
|
|
Der1s<3> der1;
|
|
Der2s<3> der2;
|
|
Lap<3> lap;
|
|
std::stringstream out;
|
|
out << der1;
|
|
EXPECT_EQ(out.str(), "Der1s<3>");
|
|
out.str("");
|
|
out.clear();
|
|
out << der2;
|
|
EXPECT_EQ(out.str(), "Der2s<3>");
|
|
out.str("");
|
|
out.clear();
|
|
out << lap;
|
|
EXPECT_EQ(out.str(), "Laplacian<3>");
|
|
}
|
|
|
|
/// [Biharmonic def]
|
|
/// Represents the Biharmonic operator
|
|
template <int dim>
|
|
struct Biharmonic : public Operator<Biharmonic<dim>> {
|
|
typedef Vec<double, dim> vec;
|
|
static std::string type_name() { return format("Biharmonic<%d>", dim); }
|
|
|
|
/// Implement only the functions you need.
|
|
template <int k>
|
|
double applyAt0(const RBFBasis<Polyharmonic<double, k>, vec>&, int index,
|
|
const std::vector<vec>& support, double scale) const {
|
|
static_assert(k != -1, "If dynamic k is desired it can be obtained from basis.rbf().");
|
|
double r = support[index].norm();
|
|
return k*(k-2)*(dim+k-2)*(dim+k-4)*ipow<k-4>(r) / ipow<4>(scale);
|
|
}
|
|
|
|
double applyAt0(const Monomials<vec>& mon, int idx, const std::vector<vec>& q, double s) const {
|
|
double result = 0;
|
|
std::array<int, dim> orders;
|
|
for (int d = 0; d < dim; ++d) { orders[d] = 0; }
|
|
for (int d = 0; d < dim; ++d) {
|
|
orders[d] = 4;
|
|
result += mon.evalOpAt0(idx, Derivative<dim>(orders), q);
|
|
orders[d] = 0;
|
|
|
|
for (int d2 = 0; d2 < d; ++d2) {
|
|
orders[d] = 2;
|
|
orders[d2] = 2;
|
|
result += 2*mon.evalOpAt0(idx, Derivative<dim>(orders), q);
|
|
orders[d] = 0;
|
|
orders[d2] = 0;
|
|
}
|
|
}
|
|
return result / ipow<4>(s);
|
|
}
|
|
};
|
|
/// [Biharmonic def]
|
|
|
|
|
|
template <int dim>
|
|
int find_idx(const Eigen::Matrix<int, dim, Eigen::Dynamic>& powers,
|
|
const Eigen::Matrix<int, dim, 1>& mon) {
|
|
int idx = -1;
|
|
for (int i = 0; i < powers.cols(); ++i) {
|
|
EXPECT_FALSE(powers.col(i) == mon && idx != -1);
|
|
if (powers.col(i) == mon) {
|
|
idx = i;
|
|
}
|
|
}
|
|
EXPECT_NE(idx, -1);
|
|
return idx;
|
|
}
|
|
|
|
TEST(Approximations, BiharmonicCostumOp) {
|
|
/// [Biharmonic usage]
|
|
typedef Vec<double, dim> vec;
|
|
BallShape<vec> b(0.0, 1.0);
|
|
double dx = 0.05;
|
|
DomainDiscretization<vec> domain = b.discretizeBoundaryWithStep(dx);
|
|
auto fn = [=](const vec&) { return dx; };
|
|
GeneralFill<vec> fill; fill.seed(1337);
|
|
fill(domain, fn);
|
|
|
|
Monomials<vec> mon(6);
|
|
Polyharmonic<double, 5> rbf;
|
|
int n = 2*mon.size();
|
|
domain.findSupport(FindClosest(n));
|
|
RBFFD<Polyharmonic<double, 5>, vec, ScaleToClosest> approx(rbf, mon);
|
|
RBFBasis<Polyharmonic<double, 5>, Vec2d> basis(n);
|
|
|
|
Biharmonic<2> bih;
|
|
double val1 = mon.evalOpAt0(0, bih);
|
|
double val2 = mon.evalOpAt0(mon.size()-1, bih);
|
|
double val3 = basis.evalOpAt0(0, bih, domain.supportNodes(0));
|
|
|
|
approx.compute(domain.pos(0), domain.supportNodes(0));
|
|
auto shape = approx.getShape(bih);
|
|
/// [Biharmonic usage]
|
|
|
|
EXPECT_EQ(0, val1);
|
|
EXPECT_EQ(0, val2);
|
|
EXPECT_EQ(225, val3);
|
|
(void) shape;
|
|
|
|
int x4 = find_idx(mon.powers(), {4, 0});
|
|
int x2y2 = find_idx(mon.powers(), {2, 2});
|
|
EXPECT_EQ(24, mon.evalOpAt0(x4, bih));
|
|
EXPECT_EQ(8, mon.evalOpAt0(x2y2, bih));
|
|
}
|
|
|
|
} // namespace mm
|
|
|
|
|