diff --git a/program/source/extemp/variable.h b/program/source/extemp/variable.h
index d775f70bdfb1be32ba1fdac4d13c7b8df8ec7394..eeb6b68a3c9c1e0aac094b7dc7fa90b20e01cd88 100644
--- a/program/source/extemp/variable.h
+++ b/program/source/extemp/variable.h
@@ -1449,7 +1449,7 @@ void Variable<DTyp>::operator= ( const Expr<A>& a ) {
for(int j = 1;j < Ny;++j )
for (int i = 1;i < Nx;++i ) {
data_quadrangles[id][i+ ( Nx+1 ) *j] =
- ao.template Give_data<quadrangleEl> ( id, i, j, 0, Nx, 0 );
+ ao.template Give_data<quadrangleEl> ( id, i, j, 0, Nx, 0 );
}
}
Update<quadrangleEl> ( id );
diff --git a/program/source/grid/blockgrid.cc b/program/source/grid/blockgrid.cc
index 7c9d71377d000901af84e9eae14e58cbfdcb97e9..88a05cbfb9086b16ed4715c8e03c650f653b3dd6 100644
--- a/program/source/grid/blockgrid.cc
+++ b/program/source/grid/blockgrid.cc
@@ -766,6 +766,15 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const
{
+ if(bg_coord != NULL) if(bg_coord->blockgrid_quad_coordinates_calculated)
+ {
+ int Nx = Give_Nx_quadrangle( id ) +1;
+ //cout << "id_hex " << id_hex << " i " << i << " j " << j << " k " << k << endl;
+ //cout << "index " << i + Nx* (j + k * (Ny))<< endl;
+ return bg_coord->blockgrid_quad_coordinates.at(id).at(i + Nx* (j ) );
+ }
+
+
double eta, xi;
D3vector PSW, PSE, PNW, PNE, PW, PE, PN, PS;
D3vector PTransFace, PWithoutTrans;
@@ -968,6 +977,10 @@ D3vector Blockgrid::Give_coord_edge ( int id, int i ) const
D3vector Blockgrid::Give_coord_point ( int id_point ) const
{
+ if(bg_coord != NULL) if(bg_coord->blockgrid_point_coordinates_calculated)
+ {
+ return bg_coord->blockgrid_point_coordinates.at(id_point);
+ }
return ug->Give_point ( id_point )->Give_coordinate();
}
@@ -1028,6 +1041,8 @@ void Blockgrid_coordinates::init_blockgrid_coordinates()
{
blockgrid_edge_coordinates_calculated = false;
blockgrid_hexa_coordinates_calculated = false;
+ blockgrid_quad_coordinates_calculated = false;
+ blockgrid_point_coordinates_calculated = false;
blockgrid_hexa_boundaries_calculated = false;
blockgrid_hexa_coordinates.resize(bg->Give_unstructured_grid()->Give_number_hexahedra());
@@ -1067,8 +1082,40 @@ void Blockgrid_coordinates::init_blockgrid_coordinates()
}
}
+ blockgrid_quad_coordinates.resize(bg->Give_unstructured_grid()->Give_number_quadrangles());
+ for (int id_quad = 0 ; id_quad < bg->Give_unstructured_grid()->Give_number_quadrangles() ; id_quad++)
+ {
+ int Nx = bg->Give_Nx_quadrangle(id_quad)+1;
+ int Ny = bg->Give_Ny_quadrangle(id_quad)+1;
+ blockgrid_quad_coordinates.at(id_quad).resize(Nx*Ny);
+ }
+ for (int id_quad = 0 ; id_quad < bg->Give_unstructured_grid()->Give_number_quadrangles() ; id_quad++)
+ {
+ int Nx = bg->Give_Nx_quadrangle(id_quad)+1;
+ int Ny = bg->Give_Ny_quadrangle(id_quad)+1;
+
+ for(int i=0;i<Nx;++i) for(int j=0;j<Ny;++j)
+ {
+ blockgrid_quad_coordinates.at(id_quad).at(i + Nx * j) = bg->Give_coord_quadrangle(id_quad,i,j);
+ }
+ }
+
+ blockgrid_point_coordinates.resize(bg->Give_unstructured_grid()->Give_number_points());
+ for (int id_point = 0 ; id_point < bg->Give_unstructured_grid()->Give_number_points() ; id_point++)
+ {
+ int N = bg->Give_unstructured_grid()->Give_number_points();
+
+ for(int i=0;i<N;++i)
+ {
+ blockgrid_point_coordinates.at(i) = bg->Give_coord_point(i);
+ }
+ }
+
+
blockgrid_edge_coordinates_calculated = true;
blockgrid_hexa_coordinates_calculated = true;
+ blockgrid_quad_coordinates_calculated = true;
+ blockgrid_point_coordinates_calculated = true;
}
void Blockgrid_coordinates::init_blockgrid_coordinates_boundary()
diff --git a/program/source/grid/blockgrid.h b/program/source/grid/blockgrid.h
index ff53b4b4db12ddf5f531146e7c7ff090cf8ed438..8d3eefb2387604bcb89fe766366196167019da1f 100644
--- a/program/source/grid/blockgrid.h
+++ b/program/source/grid/blockgrid.h
@@ -56,18 +56,28 @@ class Blockgrid_coordinates {
std::vector< std::vector< D3vector > > * getBlockgridHexaCoordinates(){return &blockgrid_hexa_coordinates;}
std::vector< std::vector< D3vector > > * getBlockgridEdgeCoordinates(){return &blockgrid_edge_coordinates;}
+ std::vector< D3vector > * getBlockgridPointCoordinates(){return &blockgrid_point_coordinates;}
+ std::vector< std::vector< D3vector > > * getBlockgridQuadCoordinates(){return &blockgrid_quad_coordinates;}
+
+
void setBlockgridHexaCoordinates(std::vector< std::vector< D3vector > > coords){blockgrid_hexa_coordinates = coords;}
void setBlockgridEdgeCoordinates(std::vector< std::vector< D3vector > > coords){blockgrid_edge_coordinates = coords;}
+ void setBlockgridPointCoordinates(std::vector< D3vector > coords) {blockgrid_point_coordinates = coords;}
+ void setBlockgridQuadCoordinates(std::vector< std::vector< D3vector > > coords){blockgrid_quad_coordinates = coords;}
std::vector< std::vector< D3vector > > blockgrid_hexa_coordinates;
std::vector< std::vector< D3vector > > blockgrid_edge_coordinates;
+ std::vector< std::vector< D3vector > > blockgrid_quad_coordinates;
+ std::vector< D3vector > blockgrid_point_coordinates;
Blockgrid *bg;
bool blockgrid_edge_coordinates_calculated;
bool blockgrid_hexa_coordinates_calculated;
+ bool blockgrid_quad_coordinates_calculated;
+ bool blockgrid_point_coordinates_calculated;
bool blockgrid_hexa_boundaries_calculated;
std::vector< std::vector< bool > > blockgrid_edge_coordinates_set;
diff --git a/program/source/interpol/interpol.cc b/program/source/interpol/interpol.cc
index 204fef3dc976fbfb8f4733e34917167bf4904a4d..4d6108ce3ab1aa8179772a538a61104364945bc6 100644
--- a/program/source/interpol/interpol.cc
+++ b/program/source/interpol/interpol.cc
@@ -682,6 +682,7 @@ Interpolate_on_block_grid::Interpolate_on_block_grid(int nx_, int ny_, int nz_,
}
+
void Interpolate_on_block_grid::interpolate(Variable<double>* U_from, Variable<double>* U_to,
double defaultInterpolation) {
@@ -706,6 +707,8 @@ void Interpolate_on_block_grid::interpolate(Variable<double>* U_from, Variable<d
(*U_to) = myFunctor(Xc,Yc,Zc);
}
+
+
double Interpolate_on_block_grid::evaluate(double coord_x, double coord_y, double coord_z) {
if(coord_x > pENT.x) return 0.0;
@@ -723,8 +726,7 @@ double Interpolate_on_block_grid::evaluate(double coord_x, double coord_y, doubl
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)];
double uESD = data[(i+1)+nx*(j +ny* k)];
double uWND = data[i +nx*((j+1)+ny* k)];
@@ -874,6 +876,31 @@ PointInterpolator::PointInterpolator(Interpolate_on_structured_grid* intermediat
hz = (pENT.z - pWSD.z) / (nz-1);
intermediateGrid->interpolate<double>(*U_from,data,defaultInterpolation_);
+}
+
+PointInterpolator::PointInterpolator(Interpolate_on_structured_grid *intermediateGrid, double defaultInterpolation_)
+{
+ defaultInterpolation = defaultInterpolation_;
+
+ nx = intermediateGrid->nx;
+ ny = intermediateGrid->ny;
+ nz = intermediateGrid->nz;
+
+ data = new double[nx*ny*nz];
+ dataCounter = new int[nx*ny*nz];
+ for (int iter = 0 ; iter < nx*ny*nz;iter++)
+ {
+ data[iter]=0.0;
+ dataCounter[iter]=0;
+ }
+
+ pENT = intermediateGrid->pENT;
+ pWSD = intermediateGrid->pWSD;
+
+ hx = intermediateGrid->getHx();
+ hy = intermediateGrid->getHy();
+ hz = intermediateGrid->getHz();
+
}
/*
@@ -942,7 +969,6 @@ double PointInterpolator::evaluate(double coord_x, double coord_y, double coord_
//cout << "id: " << id << " jd: " << jd << " kd: " << kd << endl;
//cout << "i: " << i << " j: " << j << " k: " << k << endl;
//cout << "nx: " << nx << " ny: " << ny << " nz: " << nz << endl;
-
double uWSD = data[i +nx*(j +ny* k)];
@@ -1194,6 +1220,17 @@ void PointInterpolator::smoothGrid()
}
}
+void PointInterpolator::normToNumberOfWritings()
+{
+ for (int iter = 0 ; iter < nx*ny*nz;iter++)
+ {
+ if (dataCounter[iter] >0)
+ {
+ data[iter] = data[iter] / double(dataCounter[iter]);
+ }
+ }
+}
+
void PointInterpolator::writeOnInterpolatedGrid(double coord_x, double coord_y, double coord_z, double value)
{
coord_x-=shiftx;
@@ -1228,23 +1265,39 @@ void PointInterpolator::writeOnInterpolatedGrid(double coord_x, double coord_y,
- double uWSD = value * (1.0 - locX) * (1.0 - locY) * (1.0 - locZ);
- double uESD = value * locX * (1.0 - locY) * (1.0 - locZ);
- double uWND = value * (1.0 - locX) * locY * (1.0 - locZ);
- double uEND = value * locX * locY * (1.0 - locZ);
- double uWST = value * (1.0 - locX) * (1.0 - locY) * locZ ;
- double uEST = value * locX * (1.0 - locY) * locZ ;
- double uWNT = value * (1.0 - locX) * locY * locZ ;
- double uENT = value * locX * locY * locZ ;
-
- data[i +nx*(j +ny* k)] = uWSD;
- data[(i+1)+nx*(j +ny* k)] = uESD;
- data[i +nx*((j+1)+ny* k)] = uWND;
- data[(i+1)+nx*((j+1)+ny* k)] = uEND;
- data[i +nx*(j +ny*(k+1))] = uWST;
- data[(i+1)+nx*(j +ny*(k+1))] = uEST;
- data[i +nx*((j+1)+ny*(k+1))] = uWNT;
- data[(i+1)+nx*((j+1)+ny*(k+1))] = uENT;
+// double uWSD = value * (1.0 - locX) * (1.0 - locY) * (1.0 - locZ);
+// double uESD = value * locX * (1.0 - locY) * (1.0 - locZ);
+// double uWND = value * (1.0 - locX) * locY * (1.0 - locZ);
+// double uEND = value * locX * locY * (1.0 - locZ);
+// double uWST = value * (1.0 - locX) * (1.0 - locY) * locZ ;
+// double uEST = value * locX * (1.0 - locY) * locZ ;
+// double uWNT = value * (1.0 - locX) * locY * locZ ;
+// double uENT = value * locX * locY * locZ ;
+
+// data[i +nx*(j +ny* k)] = uWSD;
+// data[(i+1)+nx*(j +ny* k)] = uESD;
+// data[i +nx*((j+1)+ny* k)] = uWND;
+// data[(i+1)+nx*((j+1)+ny* k)] = uEND;
+// data[i +nx*(j +ny*(k+1))] = uWST;
+// data[(i+1)+nx*(j +ny*(k+1))] = uEST;
+// data[i +nx*((j+1)+ny*(k+1))] = uWNT;
+// data[(i+1)+nx*((j+1)+ny*(k+1))] = uENT;
+ data[i +nx*(j +ny* k)] += value;
+ data[(i+1)+nx*(j +ny* k)] += value;
+ data[i +nx*((j+1)+ny* k)] += value;
+ data[(i+1)+nx*((j+1)+ny* k)] += value;
+ data[i +nx*(j +ny*(k+1))] += value;
+ data[(i+1)+nx*(j +ny*(k+1))] += value;
+ data[i +nx*((j+1)+ny*(k+1))] += value;
+ data[(i+1)+nx*((j+1)+ny*(k+1))] += value;
+ dataCounter[i +nx*(j +ny* k)]++;
+ dataCounter[(i+1)+nx*(j +ny* k)]++;
+ dataCounter[i +nx*((j+1)+ny* k)]++;
+ dataCounter[(i+1)+nx*((j+1)+ny* k)]++;
+ dataCounter[i +nx*(j +ny*(k+1))]++;
+ dataCounter[(i+1)+nx*(j +ny*(k+1))]++;
+ dataCounter[i +nx*((j+1)+ny*(k+1))]++;
+ dataCounter[(i+1)+nx*((j+1)+ny*(k+1))]++;
return;
@@ -1339,7 +1392,7 @@ void PointInterpolator::scaleInterpolatedData(double scale, double zeroVal)
PointInterpolator::~PointInterpolator() {
- delete interpolatorStructured;
+ //delete interpolatorStructured;
delete[] data;
}
@@ -3838,3 +3891,28 @@ int Interpolate_direct::checkForHexaNeighbours(int idHex, int i, int j, int k, D
counterHexa--;
return typ;
}
+
+Interpolate_on_block_grid_from_pointinterpolator::Interpolate_on_block_grid_from_pointinterpolator(PointInterpolator *interp, Blockgrid *blockgrid_to_)
+{
+
+ blockgrid_to = blockgrid_to_;
+ pointInterpolator = interp;
+
+}
+
+Interpolate_on_block_grid_from_pointinterpolator::~Interpolate_on_block_grid_from_pointinterpolator()
+{
+ delete pointInterpolator;
+ delete[] data;
+}
+
+void Interpolate_on_block_grid_from_pointinterpolator::interpolate(Variable<double> *U_to, double defaultInterpolation)
+{
+ Functor3<double,double,PointInterpolator> myFunctor(pointInterpolator);
+
+ X_coordinate Xc(*blockgrid_to);
+ Y_coordinate Yc(*blockgrid_to);
+ Z_coordinate Zc(*blockgrid_to);
+
+ (*U_to) = myFunctor(Xc,Yc,Zc);
+}
diff --git a/program/source/interpol/interpol.h b/program/source/interpol/interpol.h
index f3e886a780e401103f5fd07b8535b56dd0c46287..0dd1f17009fc3e0b6e64d40a31b5acb36edb6cdb 100644
--- a/program/source/interpol/interpol.h
+++ b/program/source/interpol/interpol.h
@@ -76,7 +76,15 @@ class Interpolate_on_structured_grid {
**/
template <class DTyp>
void interpolate(Variable<DTyp>& u, DTyp* data, DTyp defaultInterpolation);
- public:
+
+ int getNx(){return nx;}
+ int getNy(){return ny;}
+ int getNz(){return nz;}
+
+ double getHx(){return hx;}
+ double getHy(){return hy;}
+ double getHz(){return hz;}
+public:
int nx, ny, nz;
D3vector pENT,pWSD;
private:
@@ -200,6 +208,34 @@ void Interpolate_on_structured_grid::interpolate(Variable<DTyp>& u, DTyp* data,
/**
* interpolates data from blockgrid_from to blockgrid_to_ using an internal rectangular grid if size nx,ny,nz
**/
+class PointInterpolator;
+
+class Interpolate_on_block_grid_from_pointinterpolator {
+public:
+/**
+* If a pointinterpolator is accesible, it can be used directly to interpolate on blockgrid_to
+@param blockgrid_ of unstructured grid
+**/
+
+ Interpolate_on_block_grid_from_pointinterpolator(PointInterpolator *interp, Blockgrid* blockgrid_to_ );
+ ~Interpolate_on_block_grid_from_pointinterpolator();
+ void interpolate(Variable<double>* U_to, double defaultInterpolation = 0.0);
+
+ double evaluate(double coord_x, double coord_y, double coord_z);
+
+
+private:
+ int nx, ny, nz;
+ double hx, hy, hz;
+ PointInterpolator *pointInterpolator;
+ double* data;
+ Blockgrid* blockgrid_to;
+
+ D3vector pWSD;
+ D3vector pENT;
+
+};
+
class Interpolate_on_block_grid {
public:
/**
@@ -210,7 +246,8 @@ class Interpolate_on_block_grid {
@param blockgrid_ of unstructured grid
**/
Interpolate_on_block_grid(int nx_, int ny_, int nz_,
- Blockgrid* blockgrid_from, Blockgrid* blockgrid_to_);
+ Blockgrid* blockgrid_from, Blockgrid* blockgrid_to_);
+ Interpolate_on_block_grid(PointInterpolator *interp, Blockgrid* blockgrid_to_ );
~Interpolate_on_block_grid();
/**
@@ -220,13 +257,15 @@ class Interpolate_on_block_grid {
* verwendet.
**/
void interpolate(Variable<double>* U_from, Variable<double>* U_to, double defaultInterpolation = 0.0);
-
+
double evaluate(double coord_x, double coord_y, double coord_z);
-
+
+
private:
int nx, ny, nz;
double hx, hy, hz;
Interpolate_on_structured_grid* interpolatorStructured;
+ PointInterpolator *pointInterpolator;
double* data;
Blockgrid* blockgrid_to;
@@ -280,7 +319,7 @@ class PointInterpolator {
PointInterpolator(Interpolate_on_structured_grid* intermediateGrid,
Variable<double>* U_from, double defaultInterpolation_ = 0.0);
-
+ PointInterpolator(Interpolate_on_structured_grid* intermediateGrid, double defaultInterpolation_ = -1.0);
@@ -306,20 +345,22 @@ class PointInterpolator {
std::vector<double> evaluateGradient(double coord_x, double coord_y, double coord_z);
void smoothGrid();
-
+ void normToNumberOfWritings();
void writeOnInterpolatedGrid(double coord_x, double coord_y, double coord_z, double value);
void subtractOnInterpolatedGrid(double coord_x, double coord_y, double coord_z, double value);
void shiftInterpolatedGrid(double coord_x, double coord_y, double coord_z);
void scaleInterpolatedData(double scale, double zeroVal = 0.0);
void QPrint_VTK(QString DateiName);
-
+
+ double shiftx, shifty, shiftz;
+ double rotationx, rotationy, rotationz;
private:
int nx, ny, nz;
double hx, hy, hz;
- double shiftx, shifty, shiftz;
Interpolate_on_structured_grid* interpolatorStructured;
double* data;
+ int* dataCounter;
D3vector pWSD;
D3vector pENT;