Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62618624
Multiquadratics_2D_WENO.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, May 14, 09:29
Size
12 KB
Mime Type
text/x-c
Expires
Thu, May 16, 09:29 (2 d)
Engine
blob
Format
Raw Data
Handle
17668859
Attached To
R3721 ConSol
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_;
}
Event Timeline
Log In to Comment