Commit 3448ec54 authored by Phillip Lino Rall's avatar Phillip Lino Rall
Browse files

added the direct interpolator class. Now it is possible to interpolate in the...

added the direct interpolator class. Now it is possible to interpolate in the unstructured grid directly, without crating an intermediate mesh. slow for random access, but for structures access (as in ray tracing for example) similar to the intermediate method. Benefit: more accurate!
parent 0f270678
This diff is collapsed.
......@@ -55,7 +55,7 @@ class Blockgrid {
void Set_grid_points(int d, int Nd);
Unstructured_grid* Give_unstructured_grid() const { return ug; };
D3vector Give_coord_hexahedron(int id, int i, int j, int k) const;
D3vector Give_coord_hexahedron(int id, int i, int j, int k) const ;
D3vector Give_coord_quadrangle(int id, int i, int j) const;
D3vector Give_coord_edge(int id, int i) const;
D3vector Give_coord_point(int id) const;
......@@ -77,6 +77,10 @@ class Blockgrid {
int Total_number_of_points() const;
void init_blockgrid_coordinates();
void init_blockgrid_coordinates_boundary();
void test_blockgrid_coordinates();
void delete_blockgrid_coordinates();
/*** Komische Konstruktion wegen Phasenshift von Matthias und Christoph
* sollte vielleicht mal weg!!!!
***/
......@@ -93,7 +97,10 @@ class Blockgrid {
int getNumberPointsDegree(int d) { return number_points[d]; }
std::vector<std::vector<std::vector<int> > > blockgrid_hexa_boundary;
std::vector< std::vector< D3vector > > blockgrid_hexa_coordinates;
protected:
Blockgrid();
int *number_points; // number_points[i]
......@@ -104,6 +111,16 @@ class Blockgrid {
private:
Variable<std::complex<double> > *phase_shift;
std::vector< std::vector< bool > > blockgrid_hexa_coordinates_set;
std::vector< std::vector< D3vector > > blockgrid_edge_coordinates;
bool blockgrid_edge_coordinates_calculated = false;
bool blockgrid_hexa_coordinates_calculated = false;
bool blockgrid_hexa_boundaries_calculated = false;
std::vector< std::vector< bool > > blockgrid_edge_coordinates_set;
bool transformFromQuadrangle;
bool variable_set; // set true falls eine Variable konstruiert
......
......@@ -321,7 +321,7 @@ class Hexahedron_el : public Geometric_Object {
D3vector (*transform_hex)(int, double,double,double); // transformation of volumne
bool orientation[12];
bool orientationQuad[6];
//bool orientationQuad[6];
// geometric middle of coordinates of corner points
D3vector Give_middle_coord();
......
This diff is collapsed.
......@@ -92,6 +92,8 @@ class Interpolate_on_structured_grid {
Unstructured_grid* ug;
};
template <class DTyp>
class Give_corner_data_of_cube {
public:
......@@ -138,22 +140,22 @@ void Interpolate_on_structured_grid::interpolate(Variable<DTyp>& u, DTyp* data,
for(int ks = 0; ks < nz;++ks)
for(int js = 0; js < ny;++js)
for(int is = 0; is < nx;++is) {
ind_global = is+nx*(js+ny*ks);
i = ids_i[ind_global];
j = ids_j[ind_global];
k = ids_k[ind_global];
id_hex = ids_hex[ind_global];
/*
// test : das geht nicht::
if(id_hex < 0) {
ind_global = is+nx*((ny-1-js)+ny*ks);
i = ids_i[ind_global];
j = ids_j[ind_global];
k = ids_k[ind_global];
id_hex = ids_hex[ind_global];
}
*/
ind_global = is+nx*(js+ny*ks);
i = ids_i[ind_global];
j = ids_j[ind_global];
k = ids_k[ind_global];
id_hex = ids_hex[ind_global];
/*
// test : das geht nicht::
if(id_hex < 0) {
ind_global = is+nx*((ny-1-js)+ny*ks);
i = ids_i[ind_global];
j = ids_j[ind_global];
k = ids_k[ind_global];
id_hex = ids_hex[ind_global];
}
*/
if(id_hex < 0) data[ind_global] = defaultInterpolation;
else {
Give_corner_data_of_cube<DTyp> du(u, id_hex, i,j,k);
......@@ -172,17 +174,17 @@ void Interpolate_on_structured_grid::interpolate(Variable<DTyp>& u, DTyp* data,
<< endl;
*/
if(typ==0) data[ind_global] = interpolate_in_tet(lambda[ind_global],
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],
if(typ==1) data[ind_global] = interpolate_in_tet(lambda[ind_global],
du.EST(),du.WND(),du.WST(),du.ESD());
if(typ==2) data[ind_global] = interpolate_in_tet(lambda[ind_global],
if(typ==2) data[ind_global] = interpolate_in_tet(lambda[ind_global],
du.WND(),du.WSD(),du.WST(),du.ESD());
if(typ==3) data[ind_global] = interpolate_in_tet(lambda[ind_global],
du.EST(),du.WND(),du.ESD(),du.END());
if(typ==4) data[ind_global] = interpolate_in_tet(lambda[ind_global],
if(typ==3) data[ind_global] = interpolate_in_tet(lambda[ind_global],
du.EST(),du.WND(),du.ESD(),du.END());
if(typ==4) data[ind_global] = interpolate_in_tet(lambda[ind_global],
du.ENT(),du.WNT(),du.EST(),du.END());
if(typ==5) data[ind_global] = interpolate_in_tet(lambda[ind_global],
if(typ==5) data[ind_global] = interpolate_in_tet(lambda[ind_global],
du.WNT(),du.WND(),du.EST(),du.END());
}
}
......@@ -319,4 +321,123 @@ private:
};
//template <class DTyp>
class Interpolate_direct {
public:
Interpolate_direct (Blockgrid* bg):blockgrid(bg)
,idHexPrev(-1)
,iPrev(-1)
,jPrev(-1)
,kPrev(-1)
,idHexPrevPrev(-1)
,iPrevPrev(-1)
,jPrevPrev(-1)
,kPrevPrev(-1)
,idHexPrevPrevPrev(-1)
,iPrevPrevPrev(-1)
,jPrevPrevPrev(-1)
,kPrevPrevPrev(-1)
,counterFast(0)
,counterFastest(0)
,counterSecondTry(0)
,counterSlow(0)
,counterHexa(0)
,checkCounter(0)
,counterEdge(0)
,counterCorner(0)
,typCounter0(0)
,typCounter1(0)
,typCounter2(0)
,typCounter3(0)
,typCounter4(0)
,typCounter5(0){}
void init();
void initBoundary();
void find_cell(D3vector v);
void setPrevIndex();
void setPrevPrevIndex();
void setPrevPrevPrevIndex();
int checkBox(int idHex, int i, int j, int k, D3vector v);
int checkBoxSurrounding(int idHex, int i, int j, int k, D3vector v);
int checkBoxSurroundingOptimized(int idHex, int i, int j, int k, D3vector v);
int checkCorner(int idHex, int i, int j, int k, D3vector v);
int checkEdge(int idHex, int i, int j, int k, D3vector v);
bool checkOverlapOfBoxes(D3vector vWSD, D3vector vENT, D3vector wWSD, D3vector wENT);
void writeBox(D3vector vWSD, D3vector vENT, std::string str);
void writePoint(D3vector v, std::string str);
int checkForHexaNeighbours(int idHex, int i, int j, int k, D3vector v);
// void interpolate(Variable<DTyp>& u, DTyp* data, DTyp defaultInterpolation);
template <class DTyp>
DTyp evaluate(Variable<DTyp>& u);
int counterFast, counterFastest, counterSecondTry, counterSlow, counterHexa, counterCorner, counterEdge;
int checkCounter;
int boxCounter = 0;
int typCounter0,typCounter1,typCounter2,typCounter3,typCounter4,typCounter5;
bool debugTest;
D3vector lambda;
D3vector vNow, vPrev, vPrevPrev, vPrevPrevPrev;
private:
int idHexPrevPrevPrev, iPrevPrevPrev, jPrevPrevPrev, kPrevPrevPrev;
int idHexPrevPrev, iPrevPrev, jPrevPrev, kPrevPrev;
int idHexPrev, iPrev, jPrev, kPrev;
int idHexNow, iNow, jNow, kNow;
int typ_tet;
Blockgrid *blockgrid;
std::vector<std::vector<std::vector<D3vector> > > arrayBoxWSDENT;
std::vector<std::vector<std::vector<int> > > array_box_boundary;
// template <class DTyp>
// Variable<DTyp>* u;
};
template<class DTyp>
DTyp Interpolate_direct::evaluate(Variable<DTyp> &u)
{
int i,j,k, id_hex, typ;
int ind_global;
DTyp returnVal;
//for(int id=0;id<ug->Give_number_hexahedra();++id)
// u.template Update<hexahedronEl>(id);
if(idHexNow < 0)
{returnVal = -1;}
else {
Give_corner_data_of_cube<DTyp> du(u, idHexNow , iNow,jNow,kNow);
int typ = typ_tet;
if (lambda == D3vector(0,0,0))
{
cout << "no point found " << endl;
}
//cout << "lambda : " ; lambda.Print();cout<<endl;
if(typ==0) returnVal = interpolate_in_tet(lambda,
du.WND(),du.WNT(),du.WST(),du.EST());
if(typ==1) returnVal = interpolate_in_tet(lambda,
du.EST(),du.WND(),du.WST(),du.ESD());
if(typ==2) returnVal = interpolate_in_tet(lambda,
du.WND(),du.WSD(),du.WST(),du.ESD());
if(typ==3) returnVal = interpolate_in_tet(lambda,
du.EST(),du.WND(),du.ESD(),du.END());
if(typ==4) returnVal = interpolate_in_tet(lambda,
du.ENT(),du.WNT(),du.EST(),du.END());
if(typ==5) returnVal = interpolate_in_tet(lambda,
du.WNT(),du.WND(),du.EST(),du.END());
}
return returnVal;
};
#endif // INTERPOL_H
......@@ -67,7 +67,7 @@ D3vector lambda_of_p_in_tet(D3vector p,
D3vector t;
t = M.invert_apply(p-cA);
//cout<<"calculating lambda = ";t.Print();cout<<endl;
return M.invert_apply(p-cA);
}
......
......@@ -52,7 +52,9 @@ class D3vector {
void Print();
void Print(ofstream *Datei);
void operator=(const D3vector& v) { x=v.x; y=v.y; z=v.z; }
bool operator==(const D3vector& v) {if(v.x == x && v.y == y && v.z == z){return true;} else {return false;} }
bool operator==(const D3vector& v) {if(fabs(v.x - x) <1e-10 && fabs(v.y - y) <1e-10 && fabs(v.z - z) <1e-10 ){return true;} else {return false;} }
bool operator<(const D3vector& v) {if(v.x > x && v.y > y && v.z > z){return true;} else {return false;} }
bool operator>(const D3vector& v) {if(v.x < x && v.y < y && v.z < z){return true;} else {return false;} }
double operator[](int i) {
if(i==0) return x;
if(i==1) return y;
......@@ -82,6 +84,10 @@ inline D3vector MIN(D3vector V1, D3vector V2) {
return D3vector(MIN(V1.x,V2.x),MIN(V1.y,V2.y),MIN(V1.z,V2.z));
}
inline double SUM(D3vector V) {
return V.x+V.y+V.z;
}
inline D3vector operator+(const D3vector& v,const D3vector& w) {
return D3vector(v.x+w.x,v.y+w.y,v.z+w.z);
}
......
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