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.
 
 
 
 
 
 

43 lines
1.7 KiB

#ifndef EXAMPLES_CAHNHILLIARD_EQUATION_BIHARMONICOP_HPP_
#define EXAMPLES_CAHNHILLIARD_EQUATION_BIHARMONICOP_HPP_
#include <medusa/Medusa.hpp>
/// Represents the Biharmonic operator (see custom operator example)
template <int dim>
struct Biharmonic : public mm::Operator<Biharmonic<dim>> {
typedef mm::Vec<double, dim> vec;
static std::string type_name() { return mm::format("Biharmonic<%d>", dim); }
static std::string name() { return type_name(); }
template <int k>
double applyAt0(const mm::RBFBasis<mm::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)*mm::ipow<k-4>(r) / mm::ipow<4>(scale);
}
double applyAt0(const mm::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, mm::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, mm::Derivative<dim>(orders), q);
orders[d] = 0;
orders[d2] = 0;
}
}
return result / mm::ipow<4>(s);
}
};
#endif // EXAMPLES_CAHNHILLIARD_EQUATION_BIHARMONICOP_HPP_