#ifndef _SCTL_CHEB_UTILS_HPP_
#define _SCTL_CHEB_UTILS_HPP_
#include <algorithm> // for max, sort
#include <cassert> // for assert
#include <functional> // for function
#include "sctl/common.hpp" // for Integer, Long, SCTL_ASSERT, SCTL_NA...
#include "sctl/iterator.hpp" // for ConstIterator, Iterator
#include "sctl/math_utils.hpp" // for fabs, const_pi, cos
#include "sctl/math_utils.txx" // for pow, machine_eps
#include "sctl/matrix.hpp" // for Matrix
#include "sctl/quadrule.hpp" // for LegQuadRule
#include "sctl/static-array.hpp" // for StaticArray
#include "sctl/static-array.txx" // for StaticArray::operator[]
#include "sctl/vector.hpp" // for Vector
namespace sctl {
template <class ValueType, class Derived> class BasisInterface {
public:
template <Integer DIM> static void Nodes(Integer order, Vector<ValueType>& nodes) {
if (DIM == 1) {
Derived::Nodes1D(order, nodes);
return;
}
Vector<ValueType> nodes1d;
Derived::Nodes1D(order, nodes1d);
Integer order_DIM = pow<Integer>(order, DIM);
if (nodes.Dim() != order_DIM * DIM) nodes.ReInit(order_DIM * DIM);
StaticArray<Integer, DIM> idx;
for (Integer i = 0; i < DIM; i++) idx[i] = 0;
Integer itr = 0;
for (Integer j = 0; j < order_DIM; j++) {
for (Integer i = 0; i < DIM; i++) {
if (idx[i] == order) idx[i] = 0;
nodes[itr + i] = nodes1d[idx[i]];
}
itr += DIM;
idx[0]++;
for (Integer i = 1; i < DIM; i++) {
if (idx[i - 1] == order) idx[i]++;
}
}
}
/**
* Computes approximation from function values at node points.
* \param[in] fn_v Function values at node points (dof x order^DIM).
* \param[out] coeff Coefficient values (dof x Ncoeff).
*/
template <Integer DIM> static void Approx(Integer order, const Vector<ValueType>& fn_v, Vector<ValueType>& coeff) {
Matrix<ValueType> Mp;
{ // Precompute
static Vector<Matrix<ValueType>> precomp(1000);
SCTL_ASSERT(order < precomp.Dim());
if (precomp[order].Dim(0) * precomp[order].Dim(1) == 0) {
#pragma omp critical(SCTL_BASIS_APPROX)
if (precomp[order].Dim(0) * precomp[order].Dim(1) == 0) {
Vector<ValueType> x, p;
Derived::Nodes1D(order, x);
Derived::EvalBasis1D(order, x, p);
Matrix<ValueType> Mp1(order, order, p.begin(), false);
Mp1.pinv(machine_eps<ValueType>()).Swap(precomp[order]);
}
}
Mp.ReInit(precomp[order].Dim(0), precomp[order].Dim(1), precomp[order].begin(), false);
}
Integer order_DIM = pow<Integer>(order, DIM);
Integer order_DIM_ = pow<Integer>(order, DIM - 1);
Long dof = fn_v.Dim() / order_DIM;
SCTL_ASSERT(fn_v.Dim() == dof * order_DIM);
// Create work buffers
Long buff_size = dof * order_DIM;
Vector<ValueType> buff(2 * buff_size);
Iterator<ValueType> buff1 = buff.begin() + buff_size * 0;
Iterator<ValueType> buff2 = buff.begin() + buff_size * 1;
Vector<ValueType> fn(order_DIM * dof, (Iterator<ValueType>)fn_v.begin(), false);
for (Integer k = 0; k < DIM; k++) { // Apply Mp along k-dimension
Matrix<ValueType> Mi(dof * order_DIM_, order, fn.begin(), false);
Matrix<ValueType> Mo(dof * order_DIM_, order, buff2, false);
Matrix<ValueType>::GEMM(Mo, Mi, Mp);
Matrix<ValueType> Mo_t(order, dof * order_DIM_, buff1, false);
for (Long i = 0; i < Mo.Dim(0); i++) {
for (Long j = 0; j < Mo.Dim(1); j++) {
Mo_t[j][i] = Mo[i][j];
}
}
fn.ReInit(order_DIM * dof, buff1, false);
}
{ // Rearrange and write to coeff
Vector<ValueType> tensor(order_DIM * dof, buff1, false);
tensor2coeff<DIM>(order, tensor, coeff);
}
}
template <Integer DIM> static void Approx_(Integer order, const Vector<ValueType>& fn_v, Vector<ValueType>& coeff, ValueType scale) {
Matrix<ValueType> Mp;
{ // Precompute
static Vector<Matrix<ValueType>> precomp(1000);
SCTL_ASSERT(order < precomp.Dim());
if (precomp[order].Dim(0) * precomp[order].Dim(1) == 0) {
#pragma omp critical(SCTL_BASIS_APPROX)
if (precomp[order].Dim(0) * precomp[order].Dim(1) == 0) {
Vector<ValueType> x, p;
Derived::Nodes1D(order, x);
for (Integer i = 0; i < order; i++) x[i] = (x[i] - 0.5) * scale + 0.5;
Derived::EvalBasis1D(order, x, p);
Matrix<ValueType> Mp1(order, order, p.begin(), false);
Mp1.pinv(machine_eps<ValueType>()).Swap(precomp[order]);
}
}
Mp.ReInit(precomp[order].Dim(0), precomp[order].Dim(1), precomp[order].begin(), false);
}
Integer order_DIM = pow<Integer>(order, DIM);
Integer order_DIM_ = pow<Integer>(order, DIM - 1);
Long dof = fn_v.Dim() / order_DIM;
SCTL_ASSERT(fn_v.Dim() == dof * order_DIM);
// Create work buffers
Long buff_size = dof * order_DIM;
Vector<ValueType> buff(2 * buff_size);
Iterator<ValueType> buff1 = buff.begin() + buff_size * 0;
Iterator<ValueType> buff2 = buff.begin() + buff_size * 1;
Vector<ValueType> fn(order_DIM * dof, (Iterator<ValueType>)fn_v.begin(), false);
for (Integer k = 0; k < DIM; k++) { // Apply Mp along k-dimension
Matrix<ValueType> Mi(dof * order_DIM_, order, fn.begin(), false);
Matrix<ValueType> Mo(dof * order_DIM_, order, buff2, false);
Matrix<ValueType>::GEMM(Mo, Mi, Mp);
Matrix<ValueType> Mo_t(order, dof * order_DIM_, buff1, false);
for (Long i = 0; i < Mo.Dim(0); i++) {
for (Long j = 0; j < Mo.Dim(1); j++) {
Mo_t[j][i] = Mo[i][j];
}
}
fn.ReInit(order_DIM * dof, buff1, false);
}
{ // Rearrange and write to coeff
Vector<ValueType> tensor(order_DIM * dof, buff1, false);
tensor2coeff<DIM>(order, tensor, coeff);
}
}
/**
* Evaluates values from input coefficients at points on a regular
* grid defined by in_x, in_y, in_z the values in the input vector.
* \param[in] coeff Coefficient values (dof x Ncoeff).
* \param[out] out Values at node points (in_x[DIM-1].Dim() x ... x in_x[0].Dim() x dof).
*/
template <Integer DIM> static void Eval(Integer order, const Vector<ValueType>& coeff, ConstIterator<Vector<ValueType>> in_x, Vector<ValueType>& out) {
Integer Ncoeff = 1;
for (Integer i = 0; i < DIM; i++) {
Ncoeff = (Ncoeff * (order + i)) / (i + 1);
}
Long dof = coeff.Dim() / Ncoeff;
SCTL_ASSERT(coeff.Dim() == Ncoeff * dof);
// Precomputation
Long buff_size = dof;
StaticArray<Matrix<ValueType>, DIM> Mp;
for (Integer i = 0; i < DIM; i++) {
Integer n = in_x[i].Dim();
if (!n) return;
Mp[i].ReInit(order, n);
Vector<ValueType> p(order * n, Mp[i].begin(), false);
Derived::EvalBasis1D(order, in_x[i], p);
buff_size *= std::max(order, n);
}
// Create work buffers
Vector<ValueType> buff(2 * buff_size);
Iterator<ValueType> buff1 = buff.begin() + buff_size * 0;
Iterator<ValueType> buff2 = buff.begin() + buff_size * 1;
{ // Rearrange coefficients into a tensor.
Vector<ValueType> tensor(dof * pow<Integer>(order, DIM), buff1, false);
coeff2tensor<DIM>(order, coeff, tensor);
}
{ // ReInit out
Long len = dof;
for (Integer i = 0; i < DIM; i++) len *= in_x[i].Dim();
if (out.Dim() != len) out.ReInit(len);
}
for (Integer k = 0; k < DIM; k++) { // Apply Mp along k-dimension
Integer order_DIM = pow<Integer>(order, DIM - k - 1);
for (Integer i = 0; i < k; i++) order_DIM *= in_x[i].Dim();
Matrix<ValueType> Mi(dof * order_DIM, order, buff1, false);
Matrix<ValueType> Mo(dof * order_DIM, in_x[k].Dim(), buff2, false);
Matrix<ValueType>::GEMM(Mo, Mi, Mp[k]);
Matrix<ValueType> Mo_t(in_x[k].Dim(), dof * order_DIM, buff1, false);
if (k == DIM - 1) Mo_t.ReInit(in_x[k].Dim(), dof * order_DIM, out.begin(), false);
for (Long i = 0; i < Mo.Dim(0); i++) {
for (Long j = 0; j < Mo.Dim(1); j++) {
Mo_t[j][i] = Mo[i][j];
}
}
}
}
/**
* Returns the sum of the absolute value of coefficients of the
* highest order terms as an estimate of truncation error.
* \param[in] coeff Coefficient values (dof x Ncoeff).
*/
template <Integer DIM> static ValueType TruncErr(Integer order, const Vector<ValueType>& coeff) {
Integer Ncoeff = 1;
{ // Set Ncoeff
for (Integer i = 0; i < DIM; i++) Ncoeff = (Ncoeff * (order + i)) / (i + 1);
}
Long dof = coeff.Dim() / Ncoeff;
SCTL_ASSERT(coeff.Dim() == Ncoeff * dof);
ValueType err = 0;
for (Long l = 0; l < dof; l++) { // TODO: optimize this
Long offset0 = l * Ncoeff;
Integer indx0 = 0;
//Integer indx1 = 0;
StaticArray<Integer, DIM + 1> i0;
for (Integer i = 0; i <= DIM; i++) i0[i] = 0;
Integer sum = 0;
while (1) {
if (sum < order) {
if (sum == order - 1) err += fabs<ValueType>(coeff[offset0 + indx0]);
indx0++;
}
//indx1++;
sum++;
i0[0]++;
for (Integer j = 0; j < DIM && i0[j] == order; j++) {
i0[j] = 0;
i0[j + 1]++;
sum = sum + 1 - order;
}
if (i0[DIM]) break;
}
}
return err;
}
/**
* Compute gradient.
* \param[in] coeff_in Input coefficients (dof x Ncoeff)
* \param[out] coeff_out Output coefficients (dof x DIM x Ncoeff)
*/
template <Integer DIM> static void Grad(Integer order, const Vector<ValueType>& coeff_in, Vector<ValueType>* coeff_out) {
Integer Ncoeff = 1;
for (Integer i = 0; i < DIM; i++) {
Ncoeff = (Ncoeff * (order + i)) / (i + 1);
}
Long dof = coeff_in.Dim() / Ncoeff;
SCTL_ASSERT(coeff_in.Dim() == Ncoeff * dof);
Matrix<ValueType> Mdiff;
{ // Precompute
static Vector<Matrix<ValueType>> precomp(1000);
SCTL_ASSERT(order < precomp.Dim());
if (precomp[order].Dim(0) * precomp[order].Dim(1) == 0) {
#pragma omp critical(SCTL_BASIS_GRAD)
if (precomp[order].Dim(0) * precomp[order].Dim(1) == 0) {
Matrix<ValueType> M;
diff_1d(order, &M);
M.Swap(precomp[order]);
}
}
Mdiff.ReInit(precomp[order].Dim(0), precomp[order].Dim(1), precomp[order].begin(), false);
}
// Create work buffers
Long buff_size = dof * pow<Integer>(order, DIM);
Vector<ValueType> buff((3 + DIM) * buff_size);
Vector<ValueType> buff0(buff_size * 1, buff.begin() + buff_size * 0, false);
Vector<ValueType> buff1(buff_size * 1, buff.begin() + buff_size * 1, false);
Vector<ValueType> buff2(buff_size * 1, buff.begin() + buff_size * 2, false);
Vector<ValueType> buff3(buff_size * DIM, buff.begin() + buff_size * 3, false);
{ // buff0 <-- coeff2tensor(coeff_in);
coeff2tensor<DIM>(order, coeff_in, buff0);
}
for (Integer k = 0; k < DIM; k++) { // buff2 <-- Grad(buff0)
Long N0 = pow<Integer>(order, k);
Long N1 = order;
Long N2 = pow<Integer>(order, DIM - k - 1);
for (Long i3 = 0; i3 < dof; i3++) { // buff1 <-- transpose(buff0)
for (Long i2 = 0; i2 < N2; i2++) {
for (Long i1 = 0; i1 < N1; i1++) {
for (Long i0 = 0; i0 < N0; i0++) {
buff1[((i3 * N2 + i2) * N0 + i0) * N1 + i1] = buff0[((i3 * N2 + i2) * N1 + i1) * N0 + i0];
}
}
}
}
{ // buff2 <-- buff1 * Mdiff
Matrix<ValueType> Mi(dof * N0 * N2, N1, buff1.begin(), false);
Matrix<ValueType> Mo(dof * N0 * N2, N1, buff2.begin(), false);
Matrix<ValueType>::GEMM(Mo, Mi, Mdiff);
}
for (Long i3 = 0; i3 < dof; i3++) { // buff3 <-- transpose(buff2)
for (Long i2 = 0; i2 < N2; i2++) {
for (Long i1 = 0; i1 < N1; i1++) {
for (Long i0 = 0; i0 < N0; i0++) {
buff3[(((i2 * N1 + i1) * N0 + i0) * dof + i3) * DIM + k] = buff2[((i3 * N2 + i2) * N0 + i0) * N1 + i1];
}
}
}
}
}
{ // coeff_out <-- tensor2coeff(buff2);
tensor2coeff<DIM>(order, buff3, *coeff_out);
}
}
template <Integer DIM, Integer SUBDIM, class Kernel> static void Integ(Matrix<ValueType>& Mcoeff, Integer order, ConstIterator<ValueType> trg_, ValueType side, Integer src_face, const Kernel& ker, ValueType tol = -1, Integer Nq = 0) {
if (!Nq) Nq = order;
Integ_<DIM, SUBDIM>(Mcoeff, order, trg_, side, src_face, ker, Nq);
if (tol < 0) tol = machine_eps<ValueType>() * 256;
ValueType err = tol + 1;
Matrix<ValueType> Mtmp;
while (err > tol) {
err = 0;
ValueType max_val = pow<SUBDIM>(side);
Nq = std::max((Integer)(Nq * 1.26), Nq + 1);
Integ_<DIM, SUBDIM>(Mtmp, order, trg_, side, src_face, ker, Nq);
for (Integer i = 0; i < Mtmp.Dim(0) * Mtmp.Dim(1); i++) {
err = std::max(err, fabs<ValueType>(Mtmp[0][i] - Mcoeff[0][i]));
max_val = std::max(max_val, fabs<ValueType>(Mtmp[0][i]));
}
err /= max_val;
Mcoeff = Mtmp;
if (Nq>200) {
SCTL_WARN("Failed to converge, error = "<<err);
break;
}
}
Mcoeff = Mcoeff.Transpose();
}
template <Integer DIM> static void tensor2coeff(Integer order, const Vector<ValueType>& tensor, Vector<ValueType>& coeff) {
Integer Ncoeff = 1, Ntensor = pow<Integer>(order, DIM);
for (Integer i = 0; i < DIM; i++) Ncoeff = (Ncoeff * (order + i)) / (i + 1);
Long dof = tensor.Dim() / Ntensor;
SCTL_ASSERT(tensor.Dim() == Ntensor * dof);
if (coeff.Dim() != Ncoeff * dof) coeff.ReInit(Ncoeff * dof);
for (Long l = 0; l < dof; l++) { // TODO: optimize this
Long offset0 = l * Ncoeff;
Integer indx0 = 0;
Integer indx1 = 0;
StaticArray<Integer, DIM + 1> i0;
for (Integer i = 0; i <= DIM; i++) i0[i] = 0;
Integer sum = 0;
while (1) {
if (sum < order) {
coeff[offset0 + indx0] = tensor[l + indx1 * dof];
indx0++;
}
indx1++;
sum++;
i0[0]++;
for (Integer j = 0; j < DIM && i0[j] == order; j++) {
i0[j] = 0;
i0[j + 1]++;
sum = sum + 1 - order;
}
if (i0[DIM]) break;
}
}
}
template <Integer DIM> static void coeff2tensor(Integer order, const Vector<ValueType>& coeff, Vector<ValueType>& tensor) {
Integer Ncoeff = 1, Ntensor = pow<Integer>(order, DIM);
for (Integer i = 0; i < DIM; i++) Ncoeff = (Ncoeff * (order + i)) / (i + 1);
Long dof = coeff.Dim() / Ncoeff;
SCTL_ASSERT(coeff.Dim() == Ncoeff * dof);
if (tensor.Dim() != Ntensor * dof) tensor.ReInit(Ntensor * dof);
for (Long l = 0; l < dof; l++) { // TODO: optimize this
Long offset0 = l * Ncoeff;
Long offset1 = l * Ntensor;
Integer indx0 = 0;
Integer indx1 = 0;
StaticArray<Integer, DIM + 1> i0;
for (Integer i = 0; i <= DIM; i++) i0[i] = 0;
Integer sum = 0;
while (1) {
if (sum < order) {
tensor[offset1 + indx1] = coeff[offset0 + indx0];
indx0++;
} else {
tensor[offset1 + indx1] = 0;
}
indx1++;
sum++;
i0[0]++;
for (Integer j = 0; j < DIM && i0[j] == order; j++) {
i0[j] = 0;
i0[j + 1]++;
sum = sum + 1 - order;
}
if (i0[DIM]) break;
}
}
}
template <Integer DIM> static void Truncate(Vector<ValueType> &coeff0, Integer order0, Integer order1) {
SCTL_ASSERT(order1 <= order0);
Integer Ncoeff0 = 1, Ncoeff1 = 1;
for (Integer i = 0; i < DIM; i++) Ncoeff0 = (Ncoeff0 * (order0 + i)) / (i + 1);
for (Integer i = 0; i < DIM; i++) Ncoeff1 = (Ncoeff1 * (order1 + i)) / (i + 1);
Long dof = coeff0.Dim() / Ncoeff0;
SCTL_ASSERT(coeff0.Dim() == Ncoeff0 * dof);
Vector<ValueType> coeff1(dof * Ncoeff1);
coeff1.SetZero();
for (Long l = 0; l < dof; l++) { // TODO: optimize this
Long offset0 = l * Ncoeff0;
Long offset1 = l * Ncoeff1;
Integer indx0 = 0;
Integer indx1 = 0;
StaticArray<Integer, DIM + 1> i0;
for (Integer i = 0; i <= DIM; i++) i0[i] = 0;
Integer sum = 0;
while (1) {
if (sum < order1) coeff1[offset1 + indx1] = coeff0[offset0 + indx0];
if (sum < order0) indx0++;
if (sum < order1) indx1++;
sum++;
i0[0]++;
for (Integer j = 0; j < DIM && i0[j] == order0; j++) {
i0[j] = 0;
i0[j + 1]++;
sum = sum + 1 - order0;
}
if (i0[DIM]) break;
}
}
coeff0 = coeff1;
}
template <Integer DIM> static void Reflect(Vector<ValueType> &coeff, Integer order, Integer dir) {
SCTL_ASSERT(dir < DIM);
Integer Ncoeff = 1;
for (Integer i = 0; i < DIM; i++) Ncoeff = (Ncoeff * (order + i)) / (i + 1);
Long dof = coeff.Dim() / Ncoeff;
SCTL_ASSERT(coeff.Dim() == Ncoeff * dof);
for (Long l = 0; l < dof; l++) { // TODO: optimize this
Long offset = l * Ncoeff;
Integer indx = 0;
StaticArray<Integer, DIM + 1> i0;
for (Integer i = 0; i <= DIM; i++) i0[i] = 0;
Integer sum = 0;
while (1) {
if (sum < order) coeff[offset + indx] = coeff[offset + indx] * (i0[dir] % 2 ? -1 : 1) * (1);
if (sum < order) indx++;
sum++;
i0[0]++;
for (Integer j = 0; j < DIM && i0[j] == order; j++) {
i0[j] = 0;
i0[j + 1]++;
sum = sum + 1 - order;
}
if (i0[DIM]) break;
}
}
}
template <Integer DIM, Integer CONTINUITY> static void MakeContinuous(Vector<ValueType> &coeff0, Vector<ValueType> &coeff1, Integer order, Integer dir0, Integer dir1) {
if (dir0>=2*DIM || dir1>=2*DIM) return;
Integer Ncoeff = 1;
for (Integer i = 0; i < DIM; i++) Ncoeff = (Ncoeff * (order + i)) / (i + 1);
Long dof = coeff0.Dim() / Ncoeff;
SCTL_ASSERT(coeff0.Dim() == Ncoeff * dof);
SCTL_ASSERT(coeff1.Dim() == Ncoeff * dof);
static Matrix<Matrix<ValueType>> M(2*DIM, 2*DIM);
if (M[dir0][dir1].Dim(0) != 2 * Ncoeff) {
Integer Ngrid = pow<Integer>(order, DIM - 1);
Vector<ValueType> nodes;
Nodes<1>(order, nodes);
Matrix<ValueType> M_diff(2*Ncoeff, Ngrid);
{ // Set M_diff
M_diff.SetZero();
StaticArray<Vector<ValueType>, DIM> nodes_;
for (Integer i = 0; i < DIM; i++) { // Set nodes_
nodes_[i].ReInit(nodes.Dim(), nodes.begin(), false);
}
Vector<ValueType> nodes0, nodes1;
nodes0.PushBack(0);
nodes1.PushBack(1);
Vector<ValueType> value;
Vector<ValueType> coeff(Ncoeff);
coeff.SetZero();
for (Integer i = 0; i < Ncoeff; i++) {
coeff[i]=0.5;
value.ReInit(Ngrid, M_diff[i + Ncoeff * 0], false);
nodes_[dir0/2].ReInit(1, (dir0 & 1 ? nodes1.begin() : nodes0.begin()), false);
Eval<DIM>(order, coeff, nodes_, value);
nodes_[dir0/2].ReInit(nodes.Dim(), nodes.begin(), false);
coeff[i]=-0.5;
value.ReInit(Ngrid, M_diff[i + Ncoeff * 1], false);
nodes_[dir1/2].ReInit(1, (dir1 & 1 ? nodes1.begin() : nodes0.begin()), false);
Eval<DIM>(order, coeff, nodes_, value);
nodes_[dir1/2].ReInit(nodes.Dim(), nodes.begin(), false);
coeff[i]=0;
}
}
Matrix<ValueType> M_grad(2 * Ncoeff, 2 * Ncoeff);
{ // Set M_grad
M_grad.SetZero();
Vector<ValueType> coeff(Ncoeff * Ncoeff), coeff_grad;
coeff.SetZero();
for(Integer i = 0; i < Ncoeff; i++) coeff[i * Ncoeff + i] = 1;
Grad<DIM>(order, coeff, &coeff_grad);
for (Integer i = 0; i < Ncoeff; i++){
for (Integer j = 0; j < Ncoeff; j++){
M_grad[i + Ncoeff * 0][j + Ncoeff * 0] = coeff_grad[j + (i * DIM + dir0/2) * Ncoeff];
M_grad[i + Ncoeff * 1][j + Ncoeff * 1] = coeff_grad[j + (i * DIM + dir1/2) * Ncoeff];
}
}
}
auto fn_perturb = [&](std::function<ValueType(ValueType)> fn, bool even) { // Set M0
Matrix<ValueType> M0(Ngrid, 2 * Ncoeff);
M0.SetZero();
{ // dir0
Integer N0=pow<Integer>(order, dir0/2);
Integer N1=pow<Integer>(order, 1);
Integer N2=pow<Integer>(order, DIM - dir0/2 - 1);
SCTL_ASSERT(N0 * N2 == Ngrid);
Vector<ValueType> val(Ngrid * order), coeff;
val.SetZero();
for (Integer i0=0;i0<N0;i0++){
for (Integer i2=0;i2<N2;i2++){
for (Integer i1=0;i1<N1;i1++){
val[(i2*N1+i1)*N0+i0] = (dir0 & 1 ? fn(nodes[i1]) : fn(1.0 - nodes[i1])) * (even ? 1.0 : -1.0);
}
coeff.ReInit(Ncoeff, M0[i2 * N0 + i0] + Ncoeff * 0, false);
Approx<DIM>(order, val, coeff);
for (Integer i1=0;i1<N1;i1++){
val[(i2*N1+i1)*N0+i0] = 0;
}
}
}
}
{ // dir1
Integer N0=pow<Integer>(order, dir1/2);
Integer N1=pow<Integer>(order, 1);
Integer N2=pow<Integer>(order, DIM - dir1/2 - 1);
SCTL_ASSERT(N0 * N2 == Ngrid);
Vector<ValueType> val(Ngrid * order), coeff;
val.SetZero();
for (Integer i0=0;i0<N0;i0++){
for (Integer i2=0;i2<N2;i2++){
for (Integer i1=0;i1<N1;i1++){
val[(i2*N1+i1)*N0+i0] = (dir1 & 1 ? fn(nodes[i1]) : fn(1.0 - nodes[i1]));
}
coeff.ReInit(Ncoeff, M0[i2 * N0 + i0] + Ncoeff * 1, false);
Approx<DIM>(order, val, coeff);
for (Integer i1=0;i1<N1;i1++){
val[(i2*N1+i1)*N0+i0] = 0;
}
}
}
}
return M0;
};
if (CONTINUITY == 0) {
auto fn0 = [](ValueType x) {return x;};
Matrix<ValueType> M0 = fn_perturb(fn0, 0);
M[dir0][dir1] = M_diff * M0;
} else if (CONTINUITY == 1) {
auto fn0 = [](ValueType x) {return (-2*x + 3) * x * x;};
auto fn1 = [](ValueType x) {return (1 - x) * x * x;};
Matrix<ValueType> M0 = fn_perturb(fn0, 0);
Matrix<ValueType> M1 = fn_perturb(fn1, 1);
M[dir0][dir1] = M_diff * M0 + M_grad * M_diff * M1;
} else if (CONTINUITY == 2) {
auto fn0 = [](ValueType x) {return x*x*x*(6*x*x-15*x+10);};
auto fn1 = [](ValueType x) {return x*x*x*(-3*x*x+7*x-4);};
auto fn2 = [](ValueType x) {return x*x*x*(0.5*x*x-1*x+0.5);};
Matrix<ValueType> M0 = fn_perturb(fn0, 0);
Matrix<ValueType> M1 = fn_perturb(fn1, 1);
Matrix<ValueType> M2 = fn_perturb(fn2, 0);
M[dir0][dir1] = M_diff * M0 + M_grad * M_diff * M1 + M_grad * M_grad * M_diff * M2;
}
for (Integer i=0;i<2*Ncoeff;i++){
M[dir0][dir1][i][i]+=1.0;
}
if(0){ //// Alternate approach // DOESN'T WORK
//Matrix<ValueType> Mgrid2coeff;
//{ // Set Mgrid2coeff
// Integer Ngrid = pow<Integer>(order, DIM);
// Matrix<ValueType> M(Ngrid, Ncoeff);
// Vector<ValueType> val(Ngrid);
// val.SetZero();
// for (Integer i=0;i<Ngrid;i++) {
// val[i]=1.0;
// Vector<ValueType> coeff(Ncoeff, M[i], false);
// Approx<DIM>(order, val, coeff);
// val[i]=0.0;
// }
// Mgrid2coeff.ReInit(2*Ngrid, 2*Ncoeff);
// Mgrid2coeff.SetZero();
// for(Integer i=0;i<Ngrid;i++){
// for(Integer j=0;j<Ncoeff;j++){
// Mgrid2coeff[i+Ngrid*0][j+Ncoeff*0]=M[i][j];
// Mgrid2coeff[i+Ngrid*1][j+Ncoeff*1]=M[i][j];
// }
// }
//}
//Matrix<ValueType> Mcoeff2grid;
//{ // Set Mgrid2coeff
// StaticArray<Vector<ValueType>, DIM> nodes_;
// for (Integer i = 0; i < DIM; i++) { // Set nodes_
// nodes_[i].ReInit(nodes.Dim(), nodes.begin(), false);
// }
// Integer Ngrid = pow<Integer>(order, DIM);
// Matrix<ValueType> M(Ncoeff, Ngrid);
// Vector<ValueType> coeff(Ncoeff);
// coeff.SetZero();
// for (Integer i=0;i<Ncoeff;i++) {
// coeff[i]=1.0;
// Vector<ValueType> val(Ngrid, M[i], false);
// Eval<DIM>(order, coeff, nodes_, val);
// coeff[i]=0.0;
// }
// Mcoeff2grid.ReInit(2*Ncoeff, 2*Ngrid);
// Mcoeff2grid.SetZero();
// for(Integer i=0;i<Ncoeff;i++){
// for(Integer j=0;j<Ngrid;j++){
// Mcoeff2grid[i+Ncoeff*0][j+Ngrid*0]=M[i][j];
// Mcoeff2grid[i+Ncoeff*1][j+Ngrid*1]=M[i][j];
// }
// }
//}
//if(0){
// Integer Ngrid0 = Ngrid*order;
// Matrix<ValueType> MM(2*Ngrid0 + 2*Ngrid, 2*Ngrid0);
// MM.SetZero();
// for (Integer i=0;i<2*Ngrid0;i++) MM[i][i]=1;
// Matrix<ValueType> M0_(Ngrid, 2 * Ngrid0, MM[2 * Ngrid0 + Ngrid * 0], false); M0_ = (Mgrid2coeff * M_diff).Transpose();
// Matrix<ValueType> M1_(Ngrid, 2 * Ngrid0, MM[2 * Ngrid0 + Ngrid * 1], false); M1_ = (Mgrid2coeff * M_grad * M_diff).Transpose();
// for (Long i=0;i<2*Ngrid*2*Ngrid0;i++) MM[0][2*Ngrid0*2*Ngrid0 +i] *= 10000;
// MM = MM.Transpose().pinv();
// M[dir].ReInit(2 * Ngrid0, 2 * Ngrid0, MM.begin());
// M[dir] = Mcoeff2grid * M[dir] * Mgrid2coeff;
//} else {
// SCTL_ASSERT(DIM==2);
// Vector<ValueType> coeff_weight;
// for (Integer i=0;i<order;i++) {
// for (Integer j=0;j<order;j++) {
// if(i+j<order) coeff_weight.PushBack(pow<ValueType>(1.5, i+j)*1e-4);
// }
// }
// SCTL_ASSERT(coeff_weight.Dim()==Ncoeff);
// auto M0_ = M_diff.Transpose();
// auto M1_ = (M_grad * M_diff).Transpose();
// Matrix<ValueType> MM(2*Ncoeff + 6*Ngrid, 2*Ncoeff);
// MM.SetZero();
// for (Integer i=0;i<Ncoeff;i++) {
// MM[i+Ncoeff*0][i+Ncoeff*0]=coeff_weight[i];
// MM[i+Ncoeff*1][i+Ncoeff*1]=coeff_weight[i];
// }
// for (Integer i=0;i<Ngrid;i++) {
// for (Integer j=0;j<Ncoeff;j++) {
// MM[2 * Ncoeff + 0 * Ngrid +i][0 * Ncoeff + j] = M0_[0 * Ngrid + i][0 * Ncoeff + j];
// MM[2 * Ncoeff + 0 * Ngrid +i][1 * Ncoeff + j] = M0_[0 * Ngrid + i][1 * Ncoeff + j];
// MM[2 * Ncoeff + 1 * Ngrid +i][0 * Ncoeff + j] = M1_[0 * Ngrid + i][0 * Ncoeff + j];
// MM[2 * Ncoeff + 1 * Ngrid +i][1 * Ncoeff + j] = M1_[0 * Ngrid + i][1 * Ncoeff + j];
// MM[2 * Ncoeff + 2 * Ngrid +i][0 * Ncoeff + j] = M0_[0 * Ngrid + i][1 * Ncoeff + j];
// MM[2 * Ncoeff + 3 * Ngrid +i][1 * Ncoeff + j] = M0_[0 * Ngrid + i][0 * Ncoeff + j];
// MM[2 * Ncoeff + 4 * Ngrid +i][0 * Ncoeff + j] = M1_[0 * Ngrid + i][1 * Ncoeff + j];
// MM[2 * Ncoeff + 5 * Ngrid +i][1 * Ncoeff + j] = M1_[0 * Ngrid + i][0 * Ncoeff + j];
// }
// }
// Matrix<ValueType> MMM(2*Ncoeff + 6*Ngrid, 2*Ncoeff);
// MMM.SetZero();
// for (Integer i=0;i<Ncoeff;i++) {
// MMM[i+Ncoeff*0][i+Ncoeff*0]=coeff_weight[i];
// MMM[i+Ncoeff*1][i+Ncoeff*1]=coeff_weight[i];
// }
// for (Integer i=0;i<Ngrid;i++) {
// for (Integer j=0;j<Ncoeff;j++) {
// // MMM[2 * Ncoeff + 0 * Ngrid +i][0 * Ncoeff + j] = M0_[0 * Ngrid + i][0 * Ncoeff + j];
// // MMM[2 * Ncoeff + 0 * Ngrid +i][1 * Ncoeff + j] = M0_[0 * Ngrid + i][1 * Ncoeff + j];
// // MMM[2 * Ncoeff + 1 * Ngrid +i][0 * Ncoeff + j] = M1_[0 * Ngrid + i][0 * Ncoeff + j];
// // MMM[2 * Ncoeff + 1 * Ngrid +i][1 * Ncoeff + j] = M1_[0 * Ngrid + i][1 * Ncoeff + j];
// MMM[2 * Ncoeff + 2 * Ngrid +i][0 * Ncoeff + j] = M0_[0 * Ngrid + i][1 * Ncoeff + j];
// MMM[2 * Ncoeff + 3 * Ngrid +i][1 * Ncoeff + j] = M0_[0 * Ngrid + i][0 * Ncoeff + j];
// MMM[2 * Ncoeff + 4 * Ngrid +i][0 * Ncoeff + j] = M1_[0 * Ngrid + i][1 * Ncoeff + j];
// MMM[2 * Ncoeff + 5 * Ngrid +i][1 * Ncoeff + j] = M1_[0 * Ngrid + i][0 * Ncoeff + j];
// }
// }
// M[dir] = (MM.pinv(1e-10) * MMM).Transpose();
//}
//M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
////M[dir]=M[dir]*M[dir];
}
}
Matrix<ValueType> x(dof, 2 * Ncoeff), y(dof, 2 * Ncoeff);
for (Long i = 0; i < dof; i++) {
for (Integer j = 0; j < Ncoeff; j++) {
x[i][Ncoeff * 0 + j] = coeff0[i * Ncoeff + j];
x[i][Ncoeff * 1 + j] = coeff1[i * Ncoeff + j];
}
}
Matrix<ValueType>::GEMM(y, x, M[dir0][dir1]);
for (Long i = 0; i < dof; i++) {
for (Integer j = 0; j < Ncoeff; j++) {
coeff0[i * Ncoeff + j] = y[i][Ncoeff * 0 + j];
coeff1[i * Ncoeff + j] = y[i][Ncoeff * 1 + j];
}
}
}
template <Integer DIM, Integer CONTINUITY> static void MakeContinuousEdge(Vector<ValueType> &coeff0, Vector<ValueType> &coeff1, Integer order, Integer dir0, Integer dir1, Integer norm0, Integer norm1) {
SCTL_ASSERT(DIM==2);
if (dir0>=2*DIM || dir1>=2*DIM) return;
Integer Ncoeff = 1;
for (Integer i = 0; i < DIM; i++) Ncoeff = (Ncoeff * (order + i)) / (i + 1);
Long dof = coeff0.Dim() / Ncoeff;
SCTL_ASSERT(coeff0.Dim() == Ncoeff * dof);
SCTL_ASSERT(coeff1.Dim() == Ncoeff * dof);
static Matrix<Matrix<ValueType>> M(2*DIM, 2*DIM);
static Matrix<Matrix<ValueType>> MM(2*DIM, 2*DIM);
if (M[dir0][dir1].Dim(0) != 2 * Ncoeff) {
Integer Ngrid = pow<Integer>(order, DIM - 1);
Vector<ValueType> nodes;
Nodes<1>(order, nodes);
Matrix<ValueType> Mtrunc(2*Ncoeff, 2*Ncoeff);
{ // Set Mtrunc
Vector<ValueType> w;
w.SetZero();
for (Integer i=0;i<order;i++){
for (Integer j=0;j<order;j++){
if (i+j<order) {
w.PushBack(i<order-CONTINUITY*2-1 && j<order-CONTINUITY*2-1);
}
}
}
Mtrunc.SetZero();
for (Integer i=0;i<Ncoeff;i++){
Mtrunc[i + Ncoeff * 0][i + Ncoeff * 0] = w[i];
Mtrunc[i + Ncoeff * 1][i + Ncoeff * 1] = w[i];
}
}
Matrix<ValueType> M_diff(2*Ncoeff, Ngrid);
{ // Set M_diff
M_diff.SetZero();
StaticArray<Vector<ValueType>, DIM> nodes_;
for (Integer i = 0; i < DIM; i++) { // Set nodes_
nodes_[i].ReInit(nodes.Dim(), nodes.begin(), false);
}
Vector<ValueType> nodes0, nodes1;
nodes0.PushBack(0);
nodes1.PushBack(1);
Vector<ValueType> value;
Vector<ValueType> coeff(Ncoeff);
coeff.SetZero();
for (Integer i = 0; i < Ncoeff; i++) {
coeff[i]=0.5;
value.ReInit(Ngrid, M_diff[i + Ncoeff * 0], false);
nodes_[dir0/2].ReInit(1, (dir0 & 1 ? nodes1.begin() : nodes0.begin()), false);
Eval<DIM>(order, coeff, nodes_, value);
nodes_[dir0/2].ReInit(nodes.Dim(), nodes.begin(), false);
coeff[i]=-0.5;
value.ReInit(Ngrid, M_diff[i + Ncoeff * 1], false);
nodes_[dir1/2].ReInit(1, (dir1 & 1 ? nodes1.begin() : nodes0.begin()), false);
Eval<DIM>(order, coeff, nodes_, value);
nodes_[dir1/2].ReInit(nodes.Dim(), nodes.begin(), false);
coeff[i]=0;
}
}
Matrix<ValueType> M_grad(2 * Ncoeff, 2 * Ncoeff);
{ // Set M_grad
M_grad.SetZero();
Vector<ValueType> coeff(Ncoeff * Ncoeff), coeff_grad;
coeff.SetZero();
for(Integer i = 0; i < Ncoeff; i++) coeff[i * Ncoeff + i] = 1;
Grad<DIM>(order, coeff, &coeff_grad);
for (Integer i = 0; i < Ncoeff; i++){
for (Integer j = 0; j < Ncoeff; j++){
M_grad[i + Ncoeff * 0][j + Ncoeff * 0] = coeff_grad[j + (i * DIM + dir0/2) * Ncoeff];
M_grad[i + Ncoeff * 1][j + Ncoeff * 1] = coeff_grad[j + (i * DIM + dir1/2) * Ncoeff];
}
}
}
auto fn_perturb = [&](std::function<ValueType(ValueType)> fn, bool even) { // Set M0
Matrix<ValueType> M0(Ngrid, 2 * Ncoeff);
M0.SetZero();
{ // dir0
Integer N0=pow<Integer>(order, dir0/2);
Integer N1=pow<Integer>(order, 1);
Integer N2=pow<Integer>(order, DIM - dir0/2 - 1);
SCTL_ASSERT(N0 * N2 == Ngrid);
Vector<ValueType> val(Ngrid * order), coeff;
val.SetZero();
for (Integer i0=0;i0<N0;i0++){
for (Integer i2=0;i2<N2;i2++){
for (Integer i1=0;i1<N1;i1++){
val[(i2*N1+i1)*N0+i0] = (dir0 & 1 ? fn(nodes[i1]) : fn(1.0 - nodes[i1])) * (even ? 1.0 : -1.0);
}
coeff.ReInit(Ncoeff, M0[i2 * N0 + i0] + Ncoeff * 0, false);
Approx<DIM>(order, val, coeff);
for (Integer i1=0;i1<N1;i1++){
val[(i2*N1+i1)*N0+i0] = 0;
}
}
}
}
{ // dir1
Integer N0=pow<Integer>(order, dir1/2);
Integer N1=pow<Integer>(order, 1);
Integer N2=pow<Integer>(order, DIM - dir1/2 - 1);
SCTL_ASSERT(N0 * N2 == Ngrid);
Vector<ValueType> val(Ngrid * order), coeff;
val.SetZero();
for (Integer i0=0;i0<N0;i0++){
for (Integer i2=0;i2<N2;i2++){
for (Integer i1=0;i1<N1;i1++){
val[(i2*N1+i1)*N0+i0] = (dir1 & 1 ? fn(nodes[i1]) : fn(1.0 - nodes[i1]));
}
coeff.ReInit(Ncoeff, M0[i2 * N0 + i0] + Ncoeff * 1, false);
Approx<DIM>(order, val, coeff);
for (Integer i1=0;i1<N1;i1++){
val[(i2*N1+i1)*N0+i0] = 0;
}
}
}
}
return M0;
};
Matrix<ValueType> Mfilter[2];
{ // Set Mfilter
Mfilter[0].ReInit(2*Ncoeff, 2*Ncoeff);
Mfilter[1].ReInit(2*Ncoeff, 2*Ncoeff);
Mfilter[0].SetZero();
Mfilter[1].SetZero();
for (Integer i=0;i<Ncoeff;i++) {
Mfilter[0][i + Ncoeff * 0][i + Ncoeff * 0] = 1;
Mfilter[1][i + Ncoeff * 1][i + Ncoeff * 1] = 1;
}
}
if (CONTINUITY == 0) {
auto fn0 = [](ValueType x) {return x;};
Matrix<ValueType> M0 = fn_perturb(fn0, 0);
M[dir0][dir1] = M_diff * M0;
} else if (CONTINUITY == 1) {
auto fn0 = [](ValueType x) {return (-2*x + 3) * x * x;};
auto fn1 = [](ValueType x) {return (1 - x) * x * x;};
Matrix<ValueType> M0 = fn_perturb(fn0, 0);
Matrix<ValueType> M1 = fn_perturb(fn1, 1);
M[dir0][dir1] = M_diff * M0;
if (dir0 & 1) M[dir0][dir1] += (M_grad+M_grad) * (Mfilter[0] * M_diff * M1 * Mfilter[0]);
else M[dir0][dir1] -= (M_grad+M_grad) * (Mfilter[0] * M_diff * M1 * Mfilter[0]);
if (dir1 & 1) M[dir0][dir1] -= (M_grad+M_grad) * (Mfilter[1] * M_diff * M1 * Mfilter[1]);
else M[dir0][dir1] += (M_grad+M_grad) * (Mfilter[1] * M_diff * M1 * Mfilter[1]);
}
if (CONTINUITY == 1) {
auto fn1 = [](ValueType x) {return (1 - x) * x * x;};
Matrix<ValueType> M1 = fn_perturb(fn1, 1);
MM[dir0][dir1] = M_grad * M_diff * M1;
for (Integer i=0;i<Ncoeff;i++){
for (Integer j=0;j<2*Ncoeff;j++){
MM[dir0][dir1][i + Ncoeff*1][j]*=-1;
}
}
for (Integer i=0;i<2*Ncoeff;i++){
for (Integer j=0;j<Ncoeff;j++){
MM[dir0][dir1][i][j + Ncoeff*0]*=(dir0 & 1 ? 1.0 : -1.0);
MM[dir0][dir1][i][j + Ncoeff*1]*=(dir1 & 1 ? 1.0 : -1.0);
}
}
}
for (Integer i=0;i<2*Ncoeff;i++){
M[dir0][dir1][i][i]+=1.0;
MM[dir0][dir1][i][i]+=1.0;
}
M[dir0][dir1] = Mtrunc * M[dir0][dir1] * Mtrunc;
MM[dir0][dir1] = Mtrunc * MM[dir0][dir1] * Mtrunc;
for (Integer i=0;i<10;i++) {
M[dir0][dir1]=M[dir0][dir1]*M[dir0][dir1];
}
}
Matrix<ValueType> x(dof, 2 * Ncoeff), y(dof, 2 * Ncoeff);
for (Long i = 0; i < dof; i++) {
for (Integer j = 0; j < Ncoeff; j++) {
x[i][Ncoeff * 0 + j] = coeff0[i * Ncoeff + j];
x[i][Ncoeff * 1 + j] = coeff1[i * Ncoeff + j];
}
}
Matrix<ValueType>::GEMM(y, x, M[dir0][dir1]);
{ ////
Matrix<ValueType> xx(1, 2*Ncoeff), yy(1, 2*Ncoeff);
for (Integer j = 0; j < Ncoeff; j++) {
xx[0][Ncoeff * 0 + j] = coeff0[norm0 * Ncoeff + j];
xx[0][Ncoeff * 1 + j] = coeff1[norm1 * Ncoeff + j];
}
Matrix<ValueType>::GEMM(yy, xx, MM[dir0][dir1]);
for (Integer j = 0; j < Ncoeff; j++) {
y[norm0][Ncoeff * 0 + j] = yy[0][Ncoeff * 0 + j];
y[norm1][Ncoeff * 1 + j] = yy[0][Ncoeff * 1 + j];
}
}
for (Long i = 0; i < dof; i++) {
for (Integer j = 0; j < Ncoeff; j++) {
coeff0[i * Ncoeff + j] = y[i][Ncoeff * 0 + j];
coeff1[i * Ncoeff + j] = y[i][Ncoeff * 1 + j];
}
}
}
static void quad_rule(Integer order, Vector<ValueType>& x, Vector<ValueType>& w) {
static Vector<Vector<ValueType>> x_lst(10000);
static Vector<Vector<ValueType>> w_lst(x_lst.Dim());
SCTL_ASSERT(order < x_lst.Dim());
if (x.Dim() != order) x.ReInit(order);
if (w.Dim() != order) w.ReInit(order);
if (!order) return;
bool done = false;
#pragma omp critical(SCTL_QUAD_RULE)
if (x_lst[order].Dim()) {
Vector<ValueType>& x_ = x_lst[order];
Vector<ValueType>& w_ = w_lst[order];
for (Integer i = 0; i < order; i++) {
x[i] = x_[i];
w[i] = w_[i];
}
done = true;
}
if (done) return;
Vector<ValueType> x_, w_;
LegQuadRule<ValueType>::ComputeNdsWts(&x_, &w_, order);
#pragma omp critical(SCTL_QUAD_RULE)
if (!x_lst[order].Dim()) { // Set x_lst, w_lst
x_lst[order].Swap(x_);
w_lst[order].Swap(w_);
}
quad_rule(order, x, w);
}
private:
BasisInterface() {
void (*EvalBasis1D)(Integer, const Vector<ValueType>&, Vector<ValueType>&) = Derived::EvalBasis1D;
void (*Nodes1D)(Integer, Vector<ValueType>&) = Derived::Nodes1D;
}
static void cheb_nodes_1d(Integer order, Vector<ValueType>& nodes) {
if (nodes.Dim() != order) nodes.ReInit(order);
for (Integer i = 0; i < order; i++) {
nodes[i] = -cos<ValueType>((i + (ValueType)0.5) * const_pi<ValueType>() / order) * (ValueType)0.5 + (ValueType)0.5;
}
}
static void cheb_basis_1d(Integer order, const Vector<ValueType>& x, Vector<ValueType>& y) {
Integer n = x.Dim();
if (y.Dim() != order * n) y.ReInit(order * n);
if (order > 0) {
for (Long i = 0; i < n; i++) {
y[i] = 1.0;
}
}
if (order > 1) {
for (Long i = 0; i < n; i++) {
y[i + n] = x[i] * 2 - 1;
}
}
for (Integer i = 2; i < order; i++) {
for (Long j = 0; j < n; j++) {
y[i * n + j] = 2 * y[n + j] * y[i * n - 1 * n + j] - y[i * n - 2 * n + j];
}
}
}
template <Integer DIM, Integer SUBDIM, class Kernel> static void Integ_(Matrix<ValueType>& Mcoeff, Integer order, ConstIterator<ValueType> trg_, ValueType side, Integer src_face, const Kernel& ker, Integer Nq = 0) {
constexpr ValueType eps = machine_eps<ValueType>() * 64;
ValueType side_inv = 1.0 / side;
if (!Nq) Nq = order;
Vector<ValueType> qp, qw;
quad_rule(Nq, qp, qw);
Integer Ncoeff;
{ // Set Ncoeff
Ncoeff = 1;
for (Integer i = 0; i < SUBDIM; i++) Ncoeff = (Ncoeff * (order + i)) / (i + 1);
}
StaticArray<Integer, 2> kdim;
kdim[0] = ker.Dim(0);
kdim[1] = ker.Dim(1);
StaticArray<Integer, DIM> perm0;
StaticArray<ValueType, DIM> trg; // target after rotation
{ // Set perm0
SCTL_ASSERT(0 <= src_face && src_face < 2 * DIM);
if (SUBDIM == DIM - 1) {
for (Integer i = 0; i < DIM; i++) {
perm0[i] = (i + (src_face >> 1) + 1) % DIM;
}
} else {
for (Integer i = 0; i < DIM; i++) {
perm0[i] = i;
}
}
for (Integer i = 0; i < DIM; i++) trg[i] = trg_[perm0[i]];
if (SUBDIM == DIM - 1) trg[DIM - 1] -= side * (src_face & 1);
}
Vector<ValueType> r;
{ // Set r
Vector<ValueType> r_;
r_.PushBack(0);
for (Integer i = 0; i < SUBDIM; i++) {
r_.PushBack(fabs(trg[i] - 0.0));
r_.PushBack(fabs(trg[i] - side));
}
std::sort(r_.begin(), r_.begin() + r_.Dim());
ValueType r0, r1 = r_[r_.Dim() - 1];
r0 = (r1 > side ? r1 - side : 0.0);
for (Integer i = SUBDIM; i < DIM; i++) r0 = std::max(r0, fabs(trg[i]));
if (r0 > eps) r.PushBack(-r0);
r.PushBack(r0);
for (Integer i = 0; i < r_.Dim(); i++) {
if (r_[i] > r0) {
while (r[r.Dim() - 1] > 0.0 && 3.0 * r[r.Dim() - 1] < r_[i]) r.PushBack(3.0 * r[r.Dim() - 1]);
r.PushBack(r_[i]);
}
}
}
// Work vectors
StaticArray<Vector<ValueType>, SUBDIM> eval_mesh;
StaticArray<Vector<ValueType>, SUBDIM> eval_poly;
Vector<ValueType> eval_coord_tmp;
Vector<ValueType> eval_coord;
Vector<ValueType> kern_value;
// Temporary vectors
Vector<ValueType> r_src, n_src, v_src;
{ // Init r_src, n_src, v_src
r_src.ReInit(DIM);
for (Integer k = 0; k < DIM; k++) r_src[k] = 0.0;
if (SUBDIM == DIM - 1) {
n_src.ReInit(DIM);
for (Integer k = 0; k < DIM; k++) n_src[k] = 0.0;
n_src[src_face >> 1] = (src_face & 1 ? -1.0 : 1.0);
}
v_src.ReInit(kdim[0]);
}
Vector<ValueType> v0;
Vector<ValueType> v1;
Matrix<ValueType> Mtensor(kdim[1] * kdim[0], pow<Integer>(order, SUBDIM));
Mtensor.SetZero();
for (Integer i0 = 0; i0 < r.Dim() - 1; i0++) { // for each layer
for (Integer i1 = 0; i1 < 2 * SUBDIM; i1++) { // for each direction
StaticArray<ValueType, 2 * SUBDIM> range0;
StaticArray<ValueType, 2 * SUBDIM> range1;
{ // Set range0, range1
for (Integer k = 0; k < SUBDIM; k++) {
if (i1 >> 1 == k) {
ValueType s = (i1 & 1 ? 1.0 : -1.0);
range0[k * 2 + 0] = trg[k] + s * r[i0 + 0];
range0[k * 2 + 1] = trg[k] + s * r[i0 + 0];
range1[k * 2 + 0] = trg[k] + s * r[i0 + 1];
range1[k * 2 + 1] = trg[k] + s * r[i0 + 1];
} else {
range0[k * 2 + 0] = trg[k] - fabs(r[i0 + 0]);
range0[k * 2 + 1] = trg[k] + fabs(r[i0 + 0]);
range1[k * 2 + 0] = trg[k] - fabs(r[i0 + 1]);
range1[k * 2 + 1] = trg[k] + fabs(r[i0 + 1]);
}
}
for (Integer k = 0; k < 2 * SUBDIM; k++) {
if (range0[k] > side) range0[k] = side;
if (range0[k] < 0.0) range0[k] = 0.0;
if (range1[k] > side) range1[k] = side;
if (range1[k] < 0.0) range1[k] = 0.0;
}
bool continue_flag = false;
for (Integer k = 0; k < SUBDIM; k++) { // continue if volume if 0
if (i1 >> 1 == k) {
if (fabs(range0[2 * k + 0] - range1[2 * k + 0]) < eps && fabs(range0[2 * k + 1] - range1[2 * k + 1]) < eps) {
continue_flag = true;
break;
}
} else {
if (fabs(range0[2 * k + 0] - range0[2 * k + 1]) < eps && fabs(range1[2 * k + 0] - range1[2 * k + 1]) < eps) {
continue_flag = true;
break;
}
}
}
if (continue_flag) continue;
}
for (Integer i2 = 0; i2 < Nq; i2++) { // for each quadrature point
StaticArray<ValueType, 2 * SUBDIM> range;
for (Integer k = 0; k < 2 * SUBDIM; k++) { // Set range
range[k] = range0[k] + (range1[k] - range0[k]) * qp[i2];
}
for (Integer k = 0; k < SUBDIM; k++) { // Set eval_mesh
if (k == (i1 >> 1)) {
eval_mesh[k].ReInit(1);
eval_mesh[k][0] = range[2 * k];
} else {
eval_mesh[k].ReInit(Nq);
for (Integer l = 0; l < Nq; l++) eval_mesh[k][l] = range[2 * k + 0] + (range[2 * k + 1] - range[2 * k + 0]) * qp[l];
}
}
{ // Set eval_coord
Integer N = 1;
eval_coord.ReInit(0);
for (Integer k = 0; k < SUBDIM; k++) {
Integer Nk = eval_mesh[k].Dim();
eval_coord_tmp.Swap(eval_coord);
eval_coord.ReInit(Nk * N * DIM);
for (Integer l0 = 0; l0 < Nk; l0++) {
for (Integer l1 = 0; l1 < N; l1++) {
for (Integer l2 = 0; l2 < k; l2++) {
eval_coord[DIM * (N * l0 + l1) + l2] = eval_coord_tmp[DIM * l1 + l2];
}
eval_coord[DIM * (N * l0 + l1) + k] = trg[k] - eval_mesh[k][l0];
}
}
N *= Nk;
}
StaticArray<ValueType, DIM> c;
for (Integer k = 0; k < N; k++) { // Rotate
for (Integer l = 0; l < SUBDIM; l++) c[l] = eval_coord[k * DIM + l];
for (Integer l = SUBDIM; l < DIM; l++) c[l] = trg[l];
for (Integer l = 0; l < DIM; l++) eval_coord[k * DIM + perm0[l]] = c[l];
}
}
for (Integer k = 0; k < SUBDIM; k++) { // Set eval_poly
Integer N = eval_mesh[k].Dim();
for (Integer l = 0; l < eval_mesh[k].Dim(); l++) { // Scale eval_mesh to [0, 1]
eval_mesh[k][l] *= side_inv;
}
Derived::EvalBasis1D(order, eval_mesh[k], eval_poly[k]);
if (k == (i1 >> 1)) {
assert(N == 1);
ValueType qscal = fabs(range1[i1] - range0[i1]);
for (Integer l0 = 0; l0 < order; l0++) {
eval_poly[k][l0] *= qscal * qw[i2];
}
} else {
assert(N == Nq);
ValueType qscal = (range[2 * k + 1] - range[2 * k + 0]);
for (Integer l0 = 0; l0 < order; l0++) {
for (Integer l1 = 0; l1 < N; l1++) {
eval_poly[k][N * l0 + l1] *= qscal * qw[l1];
}
}
}
}
{ // Set kern_value
Integer N = eval_coord.Dim() / DIM;
kern_value.ReInit(kdim[0] * N * kdim[1]);
kern_value.SetZero();
for (Integer j = 0; j < kdim[0]; j++) { // Evaluate ker
for (Integer k = 0; k < kdim[0]; k++) v_src[k] = 0.0;
v_src[j] = 1.0;
Vector<ValueType> ker_value(N * kdim[1], kern_value.begin() + j * N * kdim[1], false);
ker(r_src, n_src, v_src, eval_coord, ker_value);
}
{ // Transpose
v0.ReInit(kern_value.Dim());
for (Integer k = 0; k < v0.Dim(); k++) v0[k] = kern_value[k];
Matrix<ValueType> M0(kdim[0], N * kdim[1], v0.begin(), false);
Matrix<ValueType> M1(N * kdim[1], kdim[0], kern_value.begin(), false);
for (Integer l0 = 0; l0 < M1.Dim(0); l0++) { // Transpose
for (Integer l1 = 0; l1 < M1.Dim(1); l1++) {
M1[l0][l1] = M0[l1][l0];
}
}
}
}
{ // Set Update M
Matrix<ValueType> Mkern(eval_mesh[SUBDIM - 1].Dim(), kern_value.Dim() / eval_mesh[SUBDIM - 1].Dim(), kern_value.begin(), false);
for (Integer k = SUBDIM - 1; k >= 0; k--) { // Compute v0
Matrix<ValueType> Mpoly(order, eval_mesh[k].Dim(), eval_poly[k].begin(), false);
v1.ReInit(Mpoly.Dim(0) * Mkern.Dim(1));
Matrix<ValueType> Mcoef(Mpoly.Dim(0), Mkern.Dim(1), v1.begin(), false);
Matrix<ValueType>::GEMM(Mcoef, Mpoly, Mkern);
v0.ReInit(v1.Dim());
Matrix<ValueType> Mt(Mkern.Dim(1), Mpoly.Dim(0), v0.begin(), false);
for (Integer l0 = 0; l0 < Mt.Dim(0); l0++) { // Transpose
for (Integer l1 = 0; l1 < Mt.Dim(1); l1++) {
Mt[l0][l1] = Mcoef[l1][l0];
}
}
if (k > 0) { // Reinit Mkern
Mkern.ReInit(eval_mesh[k - 1].Dim(), v0.Dim() / eval_mesh[k - 1].Dim(), v0.begin(), false);
}
}
assert(v0.Dim() == Mtensor.Dim(0) * Mtensor.Dim(1));
for (Integer k = 0; k < v0.Dim(); k++) { // Update M
Mtensor[0][k] += v0[k];
}
}
}
if (r[i0] < 0.0) break;
}
}
Mtensor = Mtensor.Transpose();
{ // Set Mcoeff
if (Mcoeff.Dim(0) != kdim[1] || Mcoeff.Dim(1) != kdim[0] * Ncoeff) {
Mcoeff.ReInit(kdim[1], kdim[0] * Ncoeff);
}
Vector<ValueType> Mtensor_(Mtensor.Dim(0) * Mtensor.Dim(1), Mtensor.begin(), false);
Vector<ValueType> Mcoeff_(Mcoeff.Dim(0) * Mcoeff.Dim(1), Mcoeff.begin(), false);
tensor2coeff<SUBDIM>(order, Mtensor_, Mcoeff_);
}
}
static void diff_1d(Integer order, Matrix<ValueType>* M) {
Vector<ValueType> nodes;
Nodes<1>(order, nodes);
Integer N = nodes.Dim();
Matrix<ValueType> M0(N, N);
for (Integer i = 0; i < N; i++) {
for (Integer j = 0; j < N; j++) {
M0[i][j] = 0;
for (Integer l = 0; l < N; l++) {
if (l != i) {
ValueType Mij = 1;
for (Integer k = 0; k < N; k++) {
if (k != i) {
if (l == k) {
Mij *= 1 / (nodes[i] - nodes[k]);
} else {
Mij *= (nodes[j] - nodes[k]) / (nodes[i] - nodes[k]);
}
}
}
M0[i][j] += Mij;
}
}
}
}
Vector<ValueType> p;
Derived::EvalBasis1D(order, nodes, p);
Matrix<ValueType> Mp(order, N, p.begin(), false);
M0 = Mp * M0;
Vector<ValueType> coeff;
Approx<1>(order, Vector<ValueType>(M0.Dim(0) * M0.Dim(1), M0.begin(), false), coeff);
(*M) = Matrix<ValueType>(M0.Dim(0), coeff.Dim() / M0.Dim(0), coeff.begin(), false);
}
friend Derived;
};
template <class ValueType> class ChebBasis : public BasisInterface<ValueType, ChebBasis<ValueType>> {
private:
ChebBasis();
static void Nodes1D(Integer order, Vector<ValueType>& nodes) { BasisInterface<ValueType, ChebBasis<ValueType>>::cheb_nodes_1d(order, nodes); }
/**
* Returns the values of all Chebyshev polynomials up to degree d,
* evaluated at points in the input vector. Output format:
* { T0[x[0]], ..., T0[x[n-1]], T1[x[0]], ..., Td[x[n-1]] }
*/
static void EvalBasis1D(Integer order, const Vector<ValueType>& x, Vector<ValueType>& y) { BasisInterface<ValueType, ChebBasis<ValueType>>::cheb_basis_1d(order, x, y); }
friend BasisInterface<ValueType, ChebBasis<ValueType>>;
};
} // end namespace
#endif // _SCTL_CHEB_UTILS_HPP_