Homec4science
Diffusion ConSol (master) Browse

ConSol/ConSol/ConSol/Reconstr/NoReconstrHybrid_2D.cpp
c23c58adfe18master

Last Change
Raw File

NoReconstrHybrid_2D.cpp
View Options

//
// Created by Fabian Moenkeberg on 2020-05-26.
//
#include "NoReconstrHybrid_2D.hpp"
void NoReconstrHybrid_2D::initialize(ModelBase *model, InterpBase *interpBase, LimiterBase *limiter, MeshBase *mesh, Source* source, Configuration config){
this->model = model;
this->source = source;
this->limiter = limiter;
this->interpBase = interpBase;
this->order = ceil((config.getOrder()+1)/2.0)+1;
this->nThreads = config.getNthreads();
if (nThreads == 0){
nThreads = omp_get_max_threads();
}
std::cout << " Reconstruction: Running threads: " << nThreads << " out of " << omp_get_max_threads() << "\n";
xCoord = mesh->getXCoord();
ngc = mesh->getNgc();
nCells = mesh->getNcells();
nQuads = QWeights.size();
nDomainsQUAD = mesh->getNdomainsQUAD() + mesh->getNdomainsQUADQUAD() + mesh->getNdomainsRQ();
nDomainsTRIQUAD = mesh->getNdomainsTRIQUAD();
nDomainsTRI = mesh->getNdomainsTRI();
nDomainsRT = mesh->getNdomainsRT();
submeshCell0 = mesh->getSubmeshCell0();
submeshEdg0 = mesh->getSubmeshEdg0();
nX = mesh->getNx();
nY = mesh->getNy();
}
void NoReconstrHybrid_2D::reconstruct(std::vector<std::vector<std::vector<std::vector<double>>>> &UReconstr, std::vector<std::vector<double>> &U,
MeshBase *mesh, double t, std::vector<std::vector<double>> &Source){
int dimSyst = mesh->getNelem();
int nCells = mesh->getNcells();
UReconstr.resize(2,std::vector<std::vector<std::vector<double>>>(1,std::vector<std::vector<double>>(mesh->getNedges(),std::vector<double>(dimSyst))));
std::vector<double> source_loc;
std::vector<int> c2e_loc(4);
int nNodes;
for (int ii = 0; ii < nCells; ii++){
c2e_loc = mesh->getC2E_internal(ii);
nNodes = c2e_loc.size();
int internal_i = mesh->getInternal(ii);
for (int j = 0; j < nNodes; j++){
if (mesh->getRight(c2e_loc[j])== internal_i){
UReconstr[0][c2e_loc[j]][0] = U[internal_i];
}else{
UReconstr[1][c2e_loc[j]][0] = U[internal_i];
}
}
if (source->getIsUsed()) {
std::vector<std::vector<double>> xy = mesh->getXCoordOfCell(internal_i);
source_loc = source->source(U[internal_i], t, std::vector<double>{(xy[0][0] + xy[1][0] + xy[2][0]) / 3,
(xy[0][1] + xy[1][1] + xy[2][1]) / 3});
for (int k = 0; k < dimSyst; k++) {
Source[internal_i][k] = source_loc[k];
}
}
}
mesh->updateEdges(UReconstr, t);
for (int ii = 0; ii < nDomainsQUAD; ii++){
omp_set_num_threads(nThreads);
#pragma omp parallel for
for (int i = 0; i < nX[ii]; i++) {
std::vector<std::vector<double>> uRec(dimSyst, std::vector<double>(2));
std::vector<double> uloc(2*(order-1) + 1);
std::vector<double> rec(2);
for (int j = 0; j < nY[ii]-1; j++) {
for (int k = 0; k < dimSyst; k++) {
uRec[k][0] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
uRec[k][1] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
}
limiter->limit_WENO(uRec, 0, mesh, U[i+ngc+(0+j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
for (int k = 0; k < dimSyst; k++) {
UReconstr[0][i+(j)*nX[ii] + submeshEdg0[ii + 1]][0][k] = uRec[k][0];
UReconstr[1][i+(j+1)*nX[ii] + submeshEdg0[ii + 1]][0][k] = uRec[k][1];
}
}
int j = nY[ii]-1;
for (int k = 0; k < dimSyst; k++) {
uRec[k][0] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
uRec[k][1] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
}
limiter->limit_WENO(uRec, 0, mesh, U[i+ngc+(0+j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
for (int k = 0; k < dimSyst; k++) {
UReconstr[0][i+(j)*nX[ii] + submeshEdg0[ii + 1]][0][k] = uRec[k][0];
UReconstr[0][i+(j+1)*nX[ii] + submeshEdg0[ii + 1]][0][k] = uRec[k][1];
}
}
omp_set_num_threads(nThreads);
#pragma omp parallel for
for (int j = 0; j < nY[ii]; j++) {
std::vector<std::vector<double>> uRec(dimSyst, std::vector<double>(2));
std::vector<double> uloc(2*(order-1) + 1);
std::vector<double> rec(2);
for (int i = 0; i < nX[ii]-1; i++) {
for (int k = 0; k < dimSyst; k++) {
uRec[k][0] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
uRec[k][1] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
}
limiter->limit_WENO(uRec, 0, mesh, U[i + ngc + 0 + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
for (int k = 0; k < dimSyst; k++) {
UReconstr[0][i+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][0];
UReconstr[1][i+1+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][1];
}
}
int i = nX[ii] - 1;
for (int k = 0; k < dimSyst; k++) {
uRec[k][0] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
uRec[k][1] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
}
limiter->limit_WENO(uRec, 0, mesh, U[i + ngc + 0 + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
for (int k = 0; k < dimSyst; k++) {
UReconstr[0][i+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][0];
UReconstr[0][i+1+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][1];
}
}
}
mesh->updateEdges(UReconstr, t);
}
c4science · Help