Commit ed88d577 authored by Phillip Lino Rall's avatar Phillip Lino Rall
Browse files

interpolator merged with christoph's version.

parents da39bbd0 04ee3325
/**********************************************************************************
* Copyright 2010 Christoph Pflaum
* Copyright 2010 Christoph Pflaum
* Department Informatik Lehrstuhl 10 - Systemsimulation
* Friedrich-Alexander Universität Erlangen-Nürnberg
*
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
......@@ -98,23 +98,23 @@ Intermadiate_grid_for_PointInterpolator::Intermadiate_grid_for_PointInterpolator
nx = nx_;
ny = ny_;
nz = nz_;
if(nx<=2) nx = 3;
if(ny<=2) ny = 3;
if(nz<=2) nz = 3;
Blockgrid* blockgrid_from = U_from->Give_blockgrid();
//Variable<double> coordXYZ(*blockgrid);
X_coordinate Xc(*blockgrid_from);
Y_coordinate Yc(*blockgrid_from);
Z_coordinate Zc(*blockgrid_from);
pWSD.x = Minimum(Xc); pWSD.y = Minimum(Yc); pWSD.z = Minimum(Zc);
pENT.x = Maximum(Xc); pENT.y = Maximum(Yc); pENT.z = Maximum(Zc);
pENT.x = Maximum(Xc); pENT.y = Maximum(Yc); pENT.z = Maximum(Zc);
interpolatorStructured = new Interpolate_on_structured_grid(nx,ny,nz, pWSD, pENT, *blockgrid_from);
}
*/
Interpolate_on_structured_grid::~Interpolate_on_structured_grid() {
......@@ -828,16 +828,16 @@ Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_,
D3vector lam;
//Variable<double> coordXYZ(*blockgrid);
X_coordinate Xc(blockgrid_);
Y_coordinate Yc(blockgrid_);
Z_coordinate Zc(blockgrid_);
//D3vector pWSD, pENT;
pWSD.x = Minimum(Xc); pWSD.y = Minimum(Yc); pWSD.z = Minimum(Zc);
pENT.x = Maximum(Xc); pENT.y = Maximum(Yc); pENT.z = Maximum(Zc);
pENT.x = Maximum(Xc); pENT.y = Maximum(Yc); pENT.z = Maximum(Zc);
blockgrid = &blockgrid_;
ug = blockgrid->Give_unstructured_grid();
......@@ -847,15 +847,15 @@ Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_,
if(nx_>1)
hx = (pENT.x - pWSD.x) / (nx_-1);
else
else
hx = 1.0;
if(ny_>1)
hy = (pENT.y - pWSD.y) / (ny_-1);
else
else
hy = 1.0;
if(nz_>1)
hz = (pENT.z - pWSD.z) / (nz_-1);
else
else
hz = 1.0;
int num_total = nx * ny * nz;
......@@ -888,13 +888,13 @@ Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_,
Nz = blockgrid->Give_Nz_hexahedron(id_hex);
for(int k=0;k<Nz;++k)
for(int j=0;j<Ny;++j)
for(int i=0;i<Nx;++i) {
// corner points of general hex-cell
for(int j=0;j<Ny;++j)
for(int i=0;i<Nx;++i) {
// corner points of general hex-cell
findLambdaForInterpolation(id_hex,i,j,k);
}
}
}
......@@ -902,9 +902,9 @@ Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_,
if(ids_hex[i]==-1) {
// wir nehmen default value!!
/*
cout << i
<< " Error: Interpolate_on_structured_grid: I cannot interpolate all data!"
<< endl;
cout << i
<< " Error: Interpolate_on_structured_grid: I cannot interpolate all data!"
<< endl;
ids_hex[i] = 0;
*/
}
......@@ -920,87 +920,87 @@ Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_,
Interpolate_on_block_grid::Interpolate_on_block_grid(int nx_, int ny_, int nz_,
Blockgrid* blockgrid_from, Blockgrid* blockgrid_to_) {
Blockgrid* blockgrid_from, Blockgrid* blockgrid_to_) {
nx = nx_;
ny = ny_;
nz = nz_;
if(nx<=2) nx = 3;
if(ny<=2) ny = 3;
if(nz<=2) nz = 3;
blockgrid_to = blockgrid_to_;
//Variable<double> coordXYZ(*blockgrid);
X_coordinate Xc(*blockgrid_to);
Y_coordinate Yc(*blockgrid_to);
Z_coordinate Zc(*blockgrid_to);
pWSD.x = Minimum(Xc); pWSD.y = Minimum(Yc); pWSD.z = Minimum(Zc);
pENT.x = Maximum(Xc); pENT.y = Maximum(Yc); pENT.z = Maximum(Zc);
pENT.x = Maximum(Xc); pENT.y = Maximum(Yc); pENT.z = Maximum(Zc);
interpolatorStructured = new Interpolate_on_structured_grid(nx,ny,nz, pWSD, pENT, *blockgrid_from);
data = new double[nx*ny*nz];
hx = (pENT.x - pWSD.x) / (nx-1);
hy = (pENT.y - pWSD.y) / (ny-1);
hz = (pENT.z - pWSD.z) / (nz-1);
/*
// test GGGG
cout << "\n WSD: " ; pWSD.Print();
cout << "\n ENT: " ; pENT.Print();
cout << "nx: " << nx << " ny: " << ny << " nz: " << nz << endl;
*/
}
void Interpolate_on_block_grid::interpolate(Variable<double>* U_from, Variable<double>* U_to,
double defaultInterpolation) {
double defaultInterpolation) {
/*
//test GGGG
X_coordinate Xfrom(*U_from->Give_blockgrid());
(*U_from) = Xfrom;
*/
(*U_from) = Xfrom;
*/
interpolatorStructured->interpolate<double>(*U_from,data,defaultInterpolation);
/*
//test GGGG
for(int i=0;i<Nx;++i) for(int j=0;j<nz;++j) for(int k=0;k<nz;++k)
for(int i=0;i<Nx;++i) for(int j=0;j<nz;++j) for(int k=0;k<nz;++k)
data[i +nx*(j +ny* k)] = hx * i;
*/
Functor3<double,double,Interpolate_on_block_grid> myFunctor(this);
X_coordinate Xc(*blockgrid_to);
Y_coordinate Yc(*blockgrid_to);
Z_coordinate Zc(*blockgrid_to);
(*U_to) = myFunctor(Xc,Yc,Zc);
}
double Interpolate_on_block_grid::evaluate(double coord_x, double coord_y, double coord_z) {
double Interpolate_on_block_grid::evaluate(double coord_x, double coord_y, double coord_z) {
if(coord_x > pENT.x) return 0.0;
if(coord_x < pWSD.x) return 0.0;
if(coord_y > pENT.y) return 0.0;
if(coord_y < pWSD.y) return 0.0;
if(coord_z > pENT.z) return 0.0;
if(coord_z < pWSD.z) return 0.0;
int i = (coord_x - pWSD.x) / hx;
int j = (coord_y - pWSD.y) / hy;
int k = (coord_z - pWSD.z) / hz;
int i = (coord_x - pWSD.x) / hx;
int j = (coord_y - pWSD.y) / hy;
int k = (coord_z - pWSD.z) / hz;
if(i < 0) i=0; if(j <0 ) j=0; if(k<0) k=0;
if(i>=nx-1) i=nx-2; if(j>=ny-1) j=ny-2; if(k>=nz-1) k=nz-2;
//cout << "i: " << i << " j: " << j << " k: " << k << endl;
double uWSD = data[i +nx*(j +ny* k)];
......@@ -1011,14 +1011,14 @@ double Interpolate_on_block_grid::evaluate(double coord_x, double coord_y, doubl
double uEST = data[(i+1)+nx*(j +ny*(k+1))];
double uWNT = data[i +nx*((j+1)+ny*(k+1))];
double uENT = data[(i+1)+nx*((j+1)+ny*(k+1))];
// assert( (i+1)+nx*((j+1)+ny*(k+1)) < nx*ny*nz);
double locX = (coord_x - pWSD.x) / hx - i;
double locY = (coord_y - pWSD.y) / hy - j;
double locZ = (coord_z - pWSD.z) / hz - k;
return uWSD * (1.0 - locX) * (1.0 - locY) * (1.0 - locZ) +
uESD * locX * (1.0 - locY) * (1.0 - locZ) +
......@@ -1029,7 +1029,7 @@ double Interpolate_on_block_grid::evaluate(double coord_x, double coord_y, doubl
uWNT * (1.0 - locX) * locY * locZ +
uENT * locX * locY * locZ;
}
Interpolate_on_block_grid::~Interpolate_on_block_grid() {
delete interpolatorStructured;
delete[] data;
......@@ -1037,10 +1037,10 @@ Interpolate_on_block_grid::~Interpolate_on_block_grid() {
/////////////////////////////////////////////////////////////
// 3. Interpolate from Variable on a blockgrid to any point using structured intermediate grid
// 3. Interpolate from Variable on a blockgrid to any point using structured intermediate grid
/////////////////////////////////////////////////////////////
PointInterpolator::PointInterpolator(int nx_, int ny_, int nz_,
Variable<double>* U_from, double defaultInterpolation_, bool trilinearInterpolation_ ) {
......@@ -1051,80 +1051,80 @@ PointInterpolator::PointInterpolator(int nx_, int ny_, int nz_,
nx = nx_;
ny = ny_;
nz = nz_;
new_lam_better(D3vector(0.033,0.799,0.316),D3vector(0.125,0.6742,0.166));
if(nx<=2) nx = 3;
if(ny<=2) ny = 3;
if(nz<=2) nz = 3;
Blockgrid* blockgrid_from = U_from->Give_blockgrid();
//Variable<double> coordXYZ(*blockgrid);
X_coordinate Xc(*blockgrid_from);
Y_coordinate Yc(*blockgrid_from);
Z_coordinate Zc(*blockgrid_from);
pWSD.x = Minimum(Xc); pWSD.y = Minimum(Yc); pWSD.z = Minimum(Zc);
pENT.x = Maximum(Xc); pENT.y = Maximum(Yc); pENT.z = Maximum(Zc);
pENT.x = Maximum(Xc); pENT.y = Maximum(Yc); pENT.z = Maximum(Zc);
interpolatorStructured = new Interpolate_on_structured_grid(nx,ny,nz, pWSD, pENT, *blockgrid_from, trilinearInterpolation_);
data = new double[nx*ny*nz];
hx = (pENT.x - pWSD.x) / (nx-1);
hy = (pENT.y - pWSD.y) / (ny-1);
hz = (pENT.z - pWSD.z) / (nz-1);
interpolatorStructured->interpolate<double>(*U_from,data,defaultInterpolation_);
/*
// test GGGG
cout << "\n WSD: " ; pWSD.Print();
cout << "\n ENT: " ; pENT.Print();
cout << "nx: " << nx << " ny: " << ny << " nz: " << nz << endl;
*/
}
PointInterpolator::PointInterpolator(int nx_, int ny_, int nz_,
D3vector pWSD_, D3vector pENT_,
Variable<double>* U_from, double defaultInterpolation_) {
D3vector pWSD_, D3vector pENT_,
Variable<double>* U_from, double defaultInterpolation_) {
defaultInterpolation = defaultInterpolation_;
shiftx = 0.0;
shifty = 0.0;
shiftz = 0.0;
nx = nx_;
ny = ny_;
nz = nz_;
if(nx<=2) nx = 3;
if(ny<=2) ny = 3;
if(nz<=2) nz = 3;
Blockgrid* blockgrid_from = U_from->Give_blockgrid();
pWSD = pWSD_;
pENT = pENT_;
pENT = pENT_;
interpolatorStructured = new Interpolate_on_structured_grid(nx,ny,nz, pWSD, pENT, *blockgrid_from);
data = new double[nx*ny*nz];
hx = (pENT.x - pWSD.x) / (nx-1);
hy = (pENT.y - pWSD.y) / (ny-1);
hz = (pENT.z - pWSD.z) / (nz-1);
interpolatorStructured->interpolate<double>(*U_from,data,defaultInterpolation);
/*
// test GGGG
cout << "\n WSD: " ; pWSD.Print();
cout << "\n ENT: " ; pENT.Print();
......@@ -1135,16 +1135,16 @@ PointInterpolator::PointInterpolator(int nx_, int ny_, int nz_,
PointInterpolator::PointInterpolator(Interpolate_on_structured_grid* intermediateGrid, Variable<double>* U_from, double defaultInterpolation_)
{
defaultInterpolation = defaultInterpolation_;
nx = intermediateGrid->nx;
ny = intermediateGrid->ny;
nz = intermediateGrid->nz;
data = new double[nx*ny*nz];
pENT = intermediateGrid->pENT;
pWSD = intermediateGrid->pWSD;
shiftx = 0.0;
shifty = 0.0;
shiftz = 0.0;
......@@ -1188,37 +1188,37 @@ PointInterpolator::PointInterpolator(Interpolate_on_structured_grid *intermediat
interpolatorStructured = intermediateGrid;
}
/*
Interpolate_on_structured_grid* PointInterpolator::intermediateGrid(int nx_, int ny_, int nz_, Variable<double>* U_from)
{
nx = nx_;
ny = ny_;
nz = nz_;
if(nx<=2) nx = 3;
if(ny<=2) ny = 3;
if(nz<=2) nz = 3;
Blockgrid* blockgrid_from = U_from->Give_blockgrid();
//Variable<double> coordXYZ(*blockgrid);
X_coordinate Xc(*blockgrid_from);
Y_coordinate Yc(*blockgrid_from);
Z_coordinate Zc(*blockgrid_from);
pWSD.x = Minimum(Xc); pWSD.y = Minimum(Yc); pWSD.z = Minimum(Zc);
pENT.x = Maximum(Xc); pENT.y = Maximum(Yc); pENT.z = Maximum(Zc);
pENT.x = Maximum(Xc); pENT.y = Maximum(Yc); pENT.z = Maximum(Zc);
interpolatorStructured = new Interpolate_on_structured_grid(nx,ny,nz, pWSD, pENT, *blockgrid_from);
return interpolatorStructured;
}
*/
double PointInterpolator::evaluate(double coord_x, double coord_y, double coord_z) {
double PointInterpolator::evaluate(double coord_x, double coord_y, double coord_z) {
coord_x-=shiftx;
coord_y-=shifty;
......@@ -1230,7 +1230,7 @@ double PointInterpolator::evaluate(double coord_x, double coord_y, double coord_
if(coord_y < (pWSD.y-(1e-9))) return defaultInterpolation;
if(coord_z > (pENT.z+(1e-9))) return defaultInterpolation;
if(coord_z < (pWSD.z-(1e-9))) return defaultInterpolation;
//cout << "coord_z " << coord_z << " pWSD.z " << pWSD.z << endl;
/*
int i = (coord_x - pWSD.x) / hx;
......@@ -1284,7 +1284,7 @@ double PointInterpolator::evaluate(double coord_x, double coord_y, double coord_
//cout << "x-1, y-1 " <<data[i-1 +nx*(j-1 +ny* k)] << endl;
// assert( (i+1)+nx*((j+1)+ny*(k+1)) < nx*ny*nz);
double posX = (coord_x - pWSD.x) ;
double locX = posX / hx - i;
double posY = (coord_y - pWSD.y);
......@@ -1296,8 +1296,8 @@ double PointInterpolator::evaluate(double coord_x, double coord_y, double coord_
//cout << "locX, Y, Z: " << locX << " " << locY << " " << locZ << endl;
//return uWSD;
//cout << "uPOS : " << uWSD << " , " << uESD << " , " << uWND << " , " << uEND << " , " << uWST << " , " << uEST << " , " << uWNT << " , " << uENT << endl;
double uTOT(0);
double uET, uWT, uWD, uED;
......@@ -1306,27 +1306,27 @@ double PointInterpolator::evaluate(double coord_x, double coord_y, double coord_
if ( (uEST != defaultInterpolation) == (uENT != defaultInterpolation) ) { uET = uEST * (1.0 - locY) + uENT * locY ;}
else if ( (uEST != defaultInterpolation) && (uENT == defaultInterpolation) ) { uET = uEST;}
else { uET = uENT;}
if ( (uWST != defaultInterpolation) == (uWNT != defaultInterpolation) ) {uWT = uWST * (1.0 - locY) + uWNT * locY ;}
else if ( (uWST != defaultInterpolation) && (uWNT == defaultInterpolation) ) {uWT = uWST;}
else {uWT = uWNT;}
if ( (uESD != defaultInterpolation) == (uEND != defaultInterpolation) ) {uED = uESD * (1.0 - locY) + uEND * locY ;}
else if ( (uESD != defaultInterpolation) && (uEND == defaultInterpolation) ) {uED = uESD;}
else {uED = uEND;}
if ( (uWSD != defaultInterpolation) == (uWND != defaultInterpolation) ) {uWD = uWSD * (1.0 - locY) + uWND * locY ;}
else if ( (uWSD != defaultInterpolation) && (uWND == defaultInterpolation) ) {uWD = uWSD;}
else {uWD = uWND;}
if ( (uET != defaultInterpolation) == (uWT != defaultInterpolation) ) {uT = uWT * (1.0 - locX) + uET * locX ;}
else if ( (uET != defaultInterpolation) && (uWT == defaultInterpolation) ) {uT = uET;}
else {uT = uWT;}
if ( (uED != defaultInterpolation) == (uWD != defaultInterpolation) ) {uD = uWD * (1.0 - locX) + uED * locX ;}
else if ( (uED != defaultInterpolation) && (uWD == defaultInterpolation) ) {uD = uED;}
else {uD = uWD;}
if ( (uT != defaultInterpolation) == (uD != defaultInterpolation) ) {uTOT = uD * (1.0 - locZ) + uT * locZ ;}
else if ( (uT != defaultInterpolation) && (uD == defaultInterpolation) ) {uTOT = uT;}
else {uTOT = uD;}
......@@ -1364,7 +1364,7 @@ double PointInterpolator::evaluate(double coord_x, double coord_y, double coord_
return uTOT;
/*
return uWSD * (1.0 - locX) * (1.0 - locY) * (1.0 - locZ) +
......@@ -1722,7 +1722,7 @@ void PointInterpolator::scaleInterpolatedData(double scale, double zeroVal)
}
}
void PointInterpolator::updateVariable(Variable<double> *U_from)
......
......@@ -196,6 +196,7 @@ void Interpolate_on_structured_grid::interpolate(Variable<DTyp>& u, DTyp* data,
<< endl;
*/
if(typ==0) data[ind_global] = interpolate_in_tet(lambda[ind_global],
du.WND(),du.WNT(),du.WST(),du.EST());
if(typ==1) data[ind_global] = interpolate_in_tet(lambda[ind_global],
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment