Homec4science
Diffusion ConSol (master) Browse

ConSol/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO.cpp
c0485f8af28dmaster

Last Change
Raw File

Multiquadratics_2D_WENO.cpp
View Options

//
// Created by Fabian Moenkeberg on 2019-12-16.
//
#include "Multiquadratics_2D_WENO.hpp"
#include "cblas.h"
Multiquadratics_2D_WENO::Multiquadratics_2D_WENO(int order, double eps){
this->order = order;
this->order_scheme = 2*((order+1)/2 + (order+1)%2);
this->shapeParam = eps;
this->dim = 2;
this->degPoly = order -1;
this->name = MULTIQUADS;
setMulti_Index(degPoly+1);
}
Multiquadratics_2D_WENO::Multiquadratics_2D_WENO(int order, int degPoly, double eps){
this->order = 1;
// this->order_scheme = 2*((order+1)/2 + (order+1)%2);
this->order_scheme = order+1;
this->shapeParam = eps;
this->dim = 2;
this->degPoly = degPoly;
this->name = MULTIQUADS;
setMulti_Index(degPoly+1);
}
std::vector<double> Multiquadratics_2D_WENO::evaluate(std::vector<std::vector<double>> evPoints, MeshBase *mesh, std::vector<int> stencil, std::vector<double> values){
int nStencils = (int) stencil.size();
std::vector<std::vector<std::vector<double>>> xB(nStencils,std::vector<std::vector<double>>(NODE_NUM, std::vector<double>(2)));
std::vector<std::vector<double>> xi(NODE_NUM,std::vector<double>(2));
for (int i = 0; i < nStencils; i++){
xi = mesh->getXCoord(mesh->getCells(stencil[i]));
xB[i] = xi;
}
double rad = stableRad / mesh->getApproxStencilSize(stencil);
// double rad = stableRad / sqrt(mesh->getSize(cellIndex));
std::vector<double> ret = Multiquadratics_2D_WENO::evaluateDirect(evPoints, xB, values, 1, rad);
return ret;
}
std::vector<double> Multiquadratics_2D_WENO::evaluateDirect(std::vector<std::vector<double>> evPoints,
std::vector<std::vector<std::vector<double>>> xB,
std::vector<double> values, double eps, double rad){
double eps0 = 0.1/rad;
eps = eps * rad;
int nEv = (int) evPoints.size();
std::vector<double> smInd(3);
std::vector<std::vector<double>> ret(3, std::vector<double>(nEv,0));
std::vector<double> res(nEv,0);
int nRbf = (int) xB.size();
std::vector<double> coeff_old;
int nRbf0 = std::min((int) xB.size(),3);
int nRbf_old;
int nP_old;
std::vector<std::vector<std::vector<double>>> xB_;
// std::vector<std::vector<std::vector<double>>> xB_.reserve(nRbf,std::vector<std::vector<double>>(3, std::vector<double>(2)));
std::vector<double> values_;
std::vector<double> values_old;
int degPolyMax_ = std::max(getDegPoly(nRbf),0);
std::vector<double> values_tmp;
for (int degPoly_ = 0; degPoly_ <= degPolyMax_; degPoly_++) {
int nP = multi_index.size();
int nP_tmp = 0;
int sum_multi_index = multi_index[0][0] + multi_index[0][1];
for (int i = 0; i < nP; i++) {
if (multi_index[i][0] + multi_index[i][1] <= degPoly_) {
nP_tmp++;
}
}
nP = nP_tmp;
int nP2 = multi_index.size();
nP_tmp = 0;
for (int i = 0; i < nP; i++) {
if (multi_index[i][0] + multi_index[i][1] <= degPoly_-2) {
nP_tmp++;
}
}
nP2 = nP_tmp;
// values[2] = 14;
if (degPoly_ == 0) {
nRbf = 1;
xB_.resize(nRbf,std::vector<std::vector<double>>(3, std::vector<double>(2)));
xB_[0] = xB[0];
values_.resize(nRbf);
values_[0] = values[0];
}else if (degPoly_ == 1) {
nRbf = 4;
xB_.resize(nRbf,std::vector<std::vector<double>>(3, std::vector<double>(2)));
values_.resize(nRbf);
for(int i = 0; i < nRbf; i++){
xB_[i] = xB[i];
values_[i] = values[i];
}
// values_[nRbf-1] = 14;
}else{
nRbf = xB.size();
// nRbf = 6;
xB_.resize(nRbf,std::vector<std::vector<double>>(3, std::vector<double>(2)));
values_.resize(nRbf);
for(int i = 0; i < nRbf; i++){
xB_[i] = xB[i];
values_[i] = values[i];
}
// values_[nRbf-1] = 14;
}
int m = nRbf + nP;
int approxOrder = std::min(degPoly_ + 1, order);
// std::vector<double> values_tmp = values_;
getInterpCoeff(xB_, values_, eps, nP, approxOrder);
// getInterpCoeff(xB_, values_tmp, eps, nP, approxOrder);
// getInterpCoeff(xB_, values_tmp, eps, 1, 1);
values_tmp = values_;
if (nRbf >1){
for ( int i = 0; i < nRbf_old; i++){
values_tmp[i] -= values_old[i];
}
for ( int i = 0; i < nP_old; i++){
values_tmp[i + nRbf] -= values_old[i + nRbf_old];
}
}
values_old = values_;
std::vector<std::vector<double>> EvA0(nEv, std::vector<double>(nRbf));
MVEvaluationMatrix(evPoints, xB_, eps, EvA0, approxOrder);
std::vector<std::vector<double>> PA(nEv, std::vector<double>(nP));
MVPolynomialEvaluationMatrixV(evPoints, xB_, PA, eps, nP);
double *EvA = new double[nEv * m];
for (int i = 0; i < nEv; i++) {
for (int j = 0; j < nRbf; j++) {
EvA[i + j * nEv] = EvA0[i][j];
}
for (int j = 0; j < nP; j++) {
EvA[i + (j + nRbf) * nEv] = PA[i][j];
}
}
if (nRbf > 1){
// smInd[degPoly_] = 1.0/(std::pow(ptwiseSmInd(xB[0], xB_, values_tmp, eps),2) + eps0) + 1;
smInd[degPoly_] = my_pow(ptwiseSmInd(xB[0], xB_, values_tmp, eps),1);
}else{
// double a1,a2,a3,b1,b2,b3;
// if (degPolyMax_>0){
// for (int i = 0; i < 3; i++){
//
// }
// }
smInd[0] = 1;
// double d = std::sqrt(std::fabs(calcSize(xB[0])));
// double tmp = 0;
// double sum_tt = 0;
// std::vector<double>a(3);
// std::vector<double>tt(3);
// for (int jj = 1; jj <= nRbf0; jj++){
// a[jj-1] = (values[jj] - values[0])/d;
// }
// double tmp2 = (std::fabs(a[1]-a[0]) + std::fabs(a[2]-a[0]) + std::fabs(a[1]-a[2]))/3;
// for (int jj = 0; jj < nRbf0; jj++){
// tt[jj] = (1 + tmp2*tmp2/(1e-10 + a[jj]))/3;
// sum_tt += tt[jj];
// }
// for (int jj = 0; jj < nRbf0; jj++){
// tmp += my_pow(tt[jj]/sum_tt*a[jj],2);
// }
// tmp = tmp*std::sqrt(std::abs(calcSize(xB[0])));
//// std::cout << tmp << "\n";
// smInd[0] = tmp;
}
double *sVecPt = values_tmp.data();
double *y = new double[nEv];
cblas_dgemv(CblasColMajor, CBLAS_TRANSPOSE::CblasNoTrans, (int) nEv, (int) (nRbf + nP), 1.0, EvA, (int) nEv,
sVecPt, 1, 0.0, y, (int) 1);
for (int i = 0; i < nEv; i++){
ret[degPoly_][i] = y[i];
}
// if (nRbf == 6){
// getInterpCoeff(xB_, values, eps, nP, approxOrder);
// double *sVecPt = values.data();
// cblas_dgemv(CblasColMajor, CBLAS_TRANSPOSE::CblasNoTrans, (int) nEv, (int) (nRbf + nP), 1.0, EvA, (int) nEv,
// sVecPt, 1, 0.0, y, (int) 1);
//
// for (int i = 0; i < nEv; i++){
// for ( int j = 1; j < degPolyMax_; j++) {
// ret[0][i] += ret[j][i];
// }
// ret[1][i] = y[i];
// }
// }
// values_old = values_;
nP_old = nP;
nRbf_old = nRbf;
delete[] y;
delete[] EvA;
// delete[] sVecPt;
}
// smInd[0] = 10;
for (int degPoly_ = 1; degPoly_ <= degPolyMax_; degPoly_++){
smInd[degPoly_] = 1/my_pow(eps0 + my_pow(100,degPoly_)*smInd[degPoly_] + smInd[0],2);
}
smInd[0] = 1/my_pow(eps0 + smInd[0] ,2);
// double sum_smInd = 0;
// for ( int i = 0; i <= degPolyMax_; i++){
// sum_smInd += smInd[i];
// }
for (int degPoly_ = 0; degPoly_ <= degPolyMax_; degPoly_++){
smInd[degPolyMax_-degPoly_] /=(smInd[0]);
// smInd[degPolyMax_-degPoly_] = 1;
// if (smInd[degPolyMax_-degPoly_] < 0.1){
// std::cout<< smInd[degPolyMax_-degPoly_];
// double t = ptwiseSmInd(xB[0], xB_, values_tmp, eps);
// }
}
for ( int j = 0; j <= degPolyMax_; j++){
for (int i = 0; i < nEv; i++){
// res[i] += (ret[degPolyMax_-j][i] - ret[degPolyMax_-j-1][i])*smInd[degPolyMax_-j]/sum_smInd*(degPolyMax_+1);
res[i] += ret[j][i]*smInd[j];
// ret[degPolyMax_-j][i] += ret[degPolyMax_-j][i]*smInd[degPolyMax_-j]/sum_smInd;
}
}
// for (int i = 0; i < nEv; i++){
// res[i] += ret[0][i]*smInd[0]/sum_smInd*(degPolyMax_+1);
// }
// if (smInd[0] > smInd[1]&&smInd[0]> smInd[2]){
// res = ret[0];
//// std::cout << 0;
// }else if( smInd[1]> smInd[2]){
// res = ret[1];
//// std::cout << 1;
// }else{
// res = ret[2];
//// std::cout << 2;
// }
// std::cout << smInd[0] << ", " << smInd[1] << ", " << smInd[2] << "\n";
return res;
}
double Multiquadratics_2D_WENO::ptwiseSmInd(std::vector<std::vector<double>> xBi,
std::vector<std::vector<std::vector<double>>> xB, std::vector<double> values,
double eps){
int nRbf = (int) xB.size();
int degPoly_ = std::max(getDegPoly(nRbf),0);
double Ind = 0;
// double s = 0;
// double p = 1;
// double dist = 0;
std::vector<double> xM(2, 0.0);
xM[0] = (xBi[0][0] + xBi[1][0] + xBi[2][0])/NODE_NUM;
xM[1] = (xBi[0][1] + xBi[1][1] + xBi[2][1])/NODE_NUM;
// double delta;
// double d_0, d_1;
// for (int i = 0; i < nRbf;i++){
// Ind += std::pow(values[i],2);
// delta = eps*std::sqrt(calcSize(xB[i]));
// Ind += values[i]*values[i];
// p *= delta*delta;
// s += delta*delta;
// d_0 = xM[0]-calcMidPt(xB[i],0);
// d_1 = xM[1]-calcMidPt(xB[i],1);
// dist += eps*std::sqrt(d_0*d_0 + d_1*d_1);
// }
// Ind = Ind*p*p/(s*s);
// Ind = 0;
// double delta0 = std::sqrt(calcSize(xBi));
double delta0 = (calcSize(xBi));
Ind*= delta0*delta0*0;
for (int j = 1;j <= degPoly_+1; j++){
// double Ds = 0;
// for (int i = 0; i < nRbf;i++){
// Ds += values[i]*DF(std::vector<double>{eps*(xM[0]-calcMidPt(xB[i],0)),eps*(xM[1]-calcMidPt(xB[i],1))},j);
// }
// Ds += DPolynomial(xM, xB, values, j, values.size() - nRbf, eps);
// Ind += std::pow(Ds,2);
for (int j2 = 0; j2 <= j; j2++){
double Ds = 0;
double Ds_tmp = 0;
for (int i = 0; i < nRbf;i++){
Ds += values[i]*DF(std::vector<double>{(xM[0]-calcMidPt(xB[i],0)),(xM[1]-calcMidPt(xB[i],1))},j2, j-j2, eps);
// Ds += values[i]*DF(std::vector<double>{eps*(xM[0]),eps*(xM[1])},j2, j-j2, eps);
// Ds_tmp += values[i]*DF(std::vector<double>{eps*(xM[0]-calcMidPt(xB[i],0)),eps*(xM[1]-calcMidPt(xB[i],1))},j2, j-j2, eps);
}
// Ds += DPolynomial(xM, xB, values, j2, j-j2, values.size() - nRbf, eps);
double Dt = DPolynomial(xM, xB, values, j2, j-j2, values.size() - nRbf, eps);
// Ind += (std::pow(my_pow(eps,j)*Ds + Dt,2)) *my_pow(delta0,j);
Ind += (std::pow(my_pow(1,j)*Ds,2) + std::pow(Dt,2)) *my_pow(delta0,j);
}
}
// if (Ind < 0.1){
// std::cout<< Ind;
// }
// for (int j2 = 0; j2 <= degPoly_; j2++){
// double Ds = 0;
// for (int i = 0; i < nRbf;i++){
// Ds += values[i]*DF(std::vector<double>{eps*(xM[0]),eps*(xM[1])},j2, degPoly_-j2, eps);
// }
//// Ds += DPolynomial(xM, xB, values, j2, degPoly_-j2, values.size() - nRbf, eps);
//// double Dt = DPolynomial(xM, xB, values, j2, j-j2, values.size() - nRbf, eps);
// Ind += std::pow(Ds,2)*my_pow(delta0,degPoly_);
// }
// Ind *= dist;
return std::fabs(Ind);
}
int Multiquadratics_2D_WENO::getDegPoly(int n){
int degPoly_;
if (false){
degPoly_ = 1;
}else{
degPoly_ = (int)floor(-1.5+sqrt(1+8*(n))/2);
}
return degPoly_;
}
c4science · Help