Commit 45de2e64 authored by Phillip Lino Rall's avatar Phillip Lino Rall
Browse files

writing single values to hexahedron possible now : set_data ...

parent f3c87616
...@@ -221,6 +221,13 @@ void Variable<DTyp>::Print_VTK(std::ostream& Datei, double (*convert)(DTyp x),do ...@@ -221,6 +221,13 @@ void Variable<DTyp>::Print_VTK(std::ostream& Datei, double (*convert)(DTyp x),do
template <class DTyp> template <class DTyp>
void Variable<DTyp>::QPrint_VTK(QString DateiName, double (*convert)(DTyp x),double stretch_z, QString title, void Variable<DTyp>::QPrint_VTK(QString DateiName, double (*convert)(DTyp x),double stretch_z, QString title,
Unstructured_grid_Marker * marker, bool cutSmallData) { Unstructured_grid_Marker * marker, bool cutSmallData) {
if (this==NULL)
{
std::cout << "could not write file, var == NULL" << std::endl;
return;
}
D3vector v; D3vector v;
int Nx, Ny, Nz; int Nx, Ny, Nz;
int num_total = 0; int num_total = 0;
......
...@@ -294,6 +294,9 @@ class Variable : public Expr<Variable< DTyp > >, public Object_based_on_ug { ...@@ -294,6 +294,9 @@ class Variable : public Expr<Variable< DTyp > >, public Object_based_on_ug {
template <elementTyp TYP_EL> template <elementTyp TYP_EL>
inline DTyp Give_data ( params_in ) const; inline DTyp Give_data ( params_in ) const;
template <elementTyp TYP_EL>
inline void Set_data ( params_in, DTyp value );
// fuer debug-zwecke mit positionsrueckgabe // fuer debug-zwecke mit positionsrueckgabe
...@@ -301,7 +304,7 @@ class Variable : public Expr<Variable< DTyp > >, public Object_based_on_ug { ...@@ -301,7 +304,7 @@ class Variable : public Expr<Variable< DTyp > >, public Object_based_on_ug {
inline DTyp Give_data ( params_in, D3vector &posVector ) const ; inline DTyp Give_data ( params_in, D3vector &posVector ) const ;
void UpdateHexahedra(); ///> Updated alle Randdaten in allen Hexahedra void UpdateHexahedra(); ///> Updated alle Randdaten in allen Hexahedra
void UpdateHexahedraBack();
/// todo update nur falls notwendig /// todo update nur falls notwendig
template <elementTyp TYP_EL> template <elementTyp TYP_EL>
void Update ( int id ) const; void Update ( int id ) const;
...@@ -411,6 +414,13 @@ void Variable<DTyp>::UpdateHexahedra() { ...@@ -411,6 +414,13 @@ void Variable<DTyp>::UpdateHexahedra() {
} }
} }
template <class DTyp>
void Variable<DTyp>::UpdateHexahedraBack() {
for(int id_hex=0;id_hex<ug->Give_number_hexahedra();++id_hex) {
Update_back(id_hex);
}
}
template <class DTyp> template <class DTyp>
int Variable<DTyp>::checkThreshold ( double threshold, int kz ) { int Variable<DTyp>::checkThreshold ( double threshold, int kz ) {
int id_hex = 1; int id_hex = 1;
...@@ -492,6 +502,12 @@ inline double Variable<double>::Give_data<hexahedronEl> ( params_in ) const { ...@@ -492,6 +502,12 @@ inline double Variable<double>::Give_data<hexahedronEl> ( params_in ) const {
return data_hexahedra[id][i+ ( Nx+1 ) * ( j+ ( Ny+1 ) *k ) ]; return data_hexahedra[id][i+ ( Nx+1 ) * ( j+ ( Ny+1 ) *k ) ];
} }
template <>
template <>
inline void Variable<double>::Set_data<hexahedronEl> ( params_in, double value ) {
data_hexahedra[id][i+ ( Nx+1 ) * ( j+ ( Ny+1 ) *k ) ] = value;
}
// fuer debug zwecke koordinatenrueckgabe // fuer debug zwecke koordinatenrueckgabe
template <> template <>
template <> template <>
......
...@@ -845,7 +845,7 @@ assert( rank_rec == 0); ...@@ -845,7 +845,7 @@ assert( rank_rec == 0);
// data_hexahedra[id_hex][i*Ni+j*Nj+Nconst] = // data_hexahedra[id_hex][i*Ni+j*Nj+Nconst] =
// data_quadrangles[id_quad][i+ ( Nxquad+1 ) *j]; // data_quadrangles[id_quad][i+ ( Nxquad+1 ) *j];
data_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] = data_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] =
data_hexahedra[id_hex][i*Ni+j*Nj+Nconst]; data_hexahedra[id_hex][i*Ni+j*Nj+Nconst];
} }
} }
else { else {
...@@ -890,7 +890,7 @@ assert( rank_rec == 0); ...@@ -890,7 +890,7 @@ assert( rank_rec == 0);
for(i=1; i<Nxedge; ++i ) { for(i=1; i<Nxedge; ++i ) {
// vorwaerts // vorwaerts
//data_hexahedra[id_hex][i*Ni+Nconst] = data_edges[id_edge][i]; //data_hexahedra[id_hex][i*Ni+Nconst] = data_edges[id_edge][i];
data_edges[id_edge][i] = data_hexahedra[id_hex][i*Ni+Nconst]; data_edges[id_edge][i] = data_hexahedra[id_hex][i*Ni+Nconst];
} }
} }
else { else {
...@@ -920,7 +920,9 @@ assert( rank_rec == 0); ...@@ -920,7 +920,9 @@ assert( rank_rec == 0);
if(rank_send == rank_rec ) // for parallel if(rank_send == rank_rec ) // for parallel
// vorwaerts // vorwaerts
//data_hexahedra[id_hex][Nconst] = data_points[id_point][0]; //data_hexahedra[id_hex][Nconst] = data_points[id_point][0];
data_points[id_point][0] = data_hexahedra[id_hex][Nconst]; {
data_points[id_point][0] = data_hexahedra[id_hex][Nconst];
}
else { else {
assert(false); assert(false);
} }
......
...@@ -42,25 +42,6 @@ ...@@ -42,25 +42,6 @@
int Blockgrid::id_count_grid = 0; int Blockgrid::id_count_grid = 0;
//remove later
void writePoint(D3vector v, std::string str, int Counter)
{
std::stringstream ss;
ss << Counter;
std::string num = ss.str();
std::string s = str + num;
ofstream myfile;
myfile.open ("C:/Users/rall/Documents/UG_Blocks_thermal/testergebnisse/" + s + ".vtk");
myfile << "# vtk DataFile Version 4.2\n";
myfile << "vtk output\n";
myfile << "ASCII\n";
myfile << "DATASET POLYDATA\n";
myfile << "POINTS 1 float\n";
myfile << v.x << " " << v.y << " " << v.z << "\n" ;
myfile << "VERTICES 1 2\n";
myfile << "1 0 \n";
myfile.close();
}
Blockgrid::Blockgrid() { Blockgrid::Blockgrid() {
...@@ -120,6 +101,7 @@ Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, int Nz ) { ...@@ -120,6 +101,7 @@ Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, int Nz ) {
Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, int Nz, int Nr ) { Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, int Nz, int Nr ) {
cout << "Nx " << Nx << "Ny " << Ny << "Nz " << Nz << "Nr " << Nr << endl;
id_of_grid = id_count_grid; ++id_count_grid; id_of_grid = id_count_grid; ++id_count_grid;
phase_shift = NULL; phase_shift = NULL;
bg_coord = NULL; bg_coord = NULL;
......
...@@ -1607,10 +1607,6 @@ void Interpolate_direct::init() ...@@ -1607,10 +1607,6 @@ void Interpolate_direct::init()
j = 0; j = 0;
for (i = 0 ; i < Nx+1; i++) for (k = 0 ; k < Nz+1; k++) for (i = 0 ; i < Nx+1; i++) for (k = 0 ; k < Nz+1; k++)
{ {
if (i == 0 && j == 0 && k == 4)
{
cout << "debuigf " << endl;
}
if (!indexAllFacesOutside.at(3).empty()) if (!indexAllFacesOutside.at(3).empty())
{calculateNeighbourIndex(correctRelation.at(3), indexAllFacesOutside.at(3).at(0).at(0), id_hex, i,j,k,Nx,Ny,Nz);} {calculateNeighbourIndex(correctRelation.at(3), indexAllFacesOutside.at(3).at(0).at(0), id_hex, i,j,k,Nx,Ny,Nz);}
correctRelation.at(3); correctRelation.at(3);
...@@ -3182,16 +3178,17 @@ int Interpolate_direct::checkEdge(int idHex, int i, int j, int k, D3vector v) ...@@ -3182,16 +3178,17 @@ int Interpolate_direct::checkEdge(int idHex, int i, int j, int k, D3vector v)
bool Interpolate_direct::checkOverlapOfBoxes(D3vector vWSD, D3vector vENT, D3vector wWSD, D3vector wENT) bool Interpolate_direct::checkOverlapOfBoxes(D3vector vWSD, D3vector vENT, D3vector wWSD, D3vector wENT)
{ {
if ( (( (vWSD.x <= wENT.x) && (vWSD.x >= wWSD.x) ) || ( (vENT.x <= wENT.x) && (vENT.x >= wWSD.x) ) ) double epsilon = 1e-10;
&& (( (vWSD.y <= wENT.y) && (vWSD.y >= wWSD.y) ) || ( (vENT.y <= wENT.y) && (vENT.y >= wWSD.y) ) ) if ( (( (vWSD.x <= wENT.x + epsilon) && (vWSD.x + epsilon >= wWSD.x) ) || ( (vENT.x <= wENT.x + epsilon) && (vENT.x + epsilon >= wWSD.x) ) )
&& (( (vWSD.z <= wENT.z) && (vWSD.z >= wWSD.z) ) || ( (vENT.z <= wENT.z) && (vENT.z >= wWSD.z) ) ) ) && (( (vWSD.y <= wENT.y + epsilon) && (vWSD.y + epsilon >= wWSD.y) ) || ( (vENT.y <= wENT.y + epsilon) && (vENT.y + epsilon >= wWSD.y) ) )
&& (( (vWSD.z <= wENT.z + epsilon) && (vWSD.z + epsilon >= wWSD.z) ) || ( (vENT.z <= wENT.z + epsilon) && (vENT.z + epsilon >= wWSD.z) ) ) )
{ {
return true; return true;
} }
if ( (( (wWSD.x <= vENT.x) && (wWSD.x >= vWSD.x) ) || ( (wENT.x <= vENT.x) && (wENT.x >= vWSD.x) ) ) if ( (( (wWSD.x <= vENT.x + epsilon) && (wWSD.x + epsilon >= vWSD.x) ) || ( (wENT.x <= vENT.x + epsilon) && (wENT.x + epsilon >= vWSD.x) ) )
&& (( (wWSD.y <= vENT.y) && (wWSD.y >= vWSD.y) ) || ( (wENT.y <= vENT.y) && (wENT.y >= vWSD.y) ) ) && (( (wWSD.y <= vENT.y + epsilon) && (wWSD.y + epsilon >= vWSD.y) ) || ( (wENT.y <= vENT.y + epsilon) && (wENT.y + epsilon >= vWSD.y) ) )
&& (( (wWSD.z <= vENT.z) && (wWSD.z >= vWSD.z) ) || ( (wENT.z <= vENT.z) && (wENT.z >= vWSD.z) ) ) ) && (( (wWSD.z <= vENT.z + epsilon) && (wWSD.z + epsilon >= vWSD.z) ) || ( (wENT.z <= vENT.z + epsilon) && (wENT.z + epsilon >= vWSD.z) ) ) )
{ {
return true; return true;
} }
......
...@@ -56,6 +56,10 @@ class D3vector { ...@@ -56,6 +56,10 @@ class D3vector {
void Print(); void Print();
void Print(std::ofstream *Datei); void Print(std::ofstream *Datei);
void operator=(const D3vector& v) { x=v.x; y=v.y; z=v.z; } void operator=(const D3vector& v) { x=v.x; y=v.y; z=v.z; }
void operator+=(const D3vector& v) { x+=v.x; y+=v.y; z+=v.z; }
void operator-=(const D3vector& v) { x-=v.x; y-=v.y; z-=v.z; }
void operator*=(const double v) { x*=v; y*=v; z*=v; }
void operator/=(const double v) { x/=v; y/=v; z/=v; }
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(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;} }
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;} }
...@@ -149,6 +153,7 @@ class D3matrix { ...@@ -149,6 +153,7 @@ class D3matrix {
double Determinante() { double Determinante() {
return (x1 * (y2*z3 - y3*z2) - y1 * (x2*z3 - x3*z2) + z1 * (x2*y3 - x3*y2)); return (x1 * (y2*z3 - y3*z2) - y1 * (x2*z3 - x3*z2) + z1 * (x2*y3 - x3*y2));
} }
/* /*
D3matrix(D3vector cx, D3vector cy, D3vector cz) { D3matrix(D3vector cx, D3vector cy, D3vector cz) {
x1 = cx.x; y1 = cy.x; z1 = cz.x; x1 = cx.x; y1 = cy.x; z1 = cz.x;
...@@ -178,7 +183,7 @@ class D3matrix { ...@@ -178,7 +183,7 @@ class D3matrix {
D3vector invert_apply(const D3vector& v) { D3vector invert_apply(const D3vector& v) {
double x,y,z; double x,y,z;
double det = Determinante(); //x1 * (y2*z3 - y3*z2) - y1 * (x2*z3 - x3*z2) + z1 * (x2*y3 - x3*y2); double det = Determinante(); //x1 * (y2*z3 - y3*z2) - y1 * (x2*z3 - x3*z2) + z1 * (x2*y3 - x3*y2);
// std::cout << "det "<< det << std::endl;
/* /*
x = ( (y2*z3 - y3*z2) * v.x - (x2*z3 - x3*z2) * v.y + (x2*y3 - x3*y2) * v.z ) / det; x = ( (y2*z3 - y3*z2) * v.x - (x2*z3 - x3*z2) * v.y + (x2*y3 - x3*y2) * v.z ) / det;
y = (- (y1*z3 - y3*z1) * v.x + (x1*z3 - x3*z1) * v.y - (x1*y3 - x3*y1) * v.z ) / det; y = (- (y1*z3 - y3*z1) * v.x + (x1*z3 - x3*z1) * v.y - (x1*y3 - x3*y1) * v.z ) / det;
...@@ -192,6 +197,176 @@ class D3matrix { ...@@ -192,6 +197,176 @@ class D3matrix {
return D3vector(x,y,z); return D3vector(x,y,z);
} }
inline void Print() {
std::cout << "Matrix: " << std::endl;
std::cout << x1 << ", " << y1 << ", " << z1 << ";" <<std::endl;
std::cout << x2 << ", " << y2 << ", " << z2 << ";" <<std::endl;
std::cout << x3 << ", " << y3 << ", " << z3 << ";" <<std::endl;
}
D3matrix matrixMultiply(D3matrix m)
{
D3vector cx = {x1 * m.x1 + y1 * m.x2 +z1 * m.x3,
x2 * m.x1 + y2 * m.x2 +z2 * m.x3,
x3 * m.x1 + y3 * m.x2 +z3 * m.x3};
D3vector cy = {x1 * m.y1 + y1 * m.y2 +z1 * m.y3,
x2 * m.y1 + y2 * m.y2 +z2 * m.y3,
x3 * m.y1 + y3 * m.y2 +z3 * m.y3};
D3vector cz = {x1 * m.z1 + y1 * m.z2 +z1 * m.z3,
x2 * m.z1 + y2 * m.z2 +z2 * m.z3,
x3 * m.z1 + y3 * m.z2 +z3 * m.z3};
return D3matrix(cx,cy,cz);
}
D3vector vectorMultiply(D3vector m)
{
// D3vector res = {x1 * m.x + x2 * m.x + x3 * m.x,
// y1 * m.y + y2 * m.y + y3 * m.y,
// z1 * m.z + z2 * m.z + z3 * m.z};
// D3vector res = {x1 * m.x + x2 * m.y + x3 * m.z,
// y1 * m.x + y2 * m.y + y3 * m.z,
// z1 * m.x + z2 * m.y + z3 * m.z};
D3vector res = {x1 * m.x + y1 * m.y + z1 * m.z,
x2 * m.x + y2 * m.y + z2 * m.z,
x3 * m.x + y3 * m.y + z3 * m.z};
return res;
}
void transpose()
{
double temp;
temp = x2; x2 = y1 ; y1 = temp;
temp = x3; x3 = z1 ; z1 = temp;
temp = y3; y3 = z2 ; z2 = temp;
}
void invert_gauss_elimination()
{
int n;
double a[3][6];
for (int i = 0 ; i<3;i++)
{
for (int j = 0 ; j<6;j++)
{
//std::cout << "i j " << i << " " << j << std::endl;
a[i][j] = 0;
}
}
int order = 3;
n = 3;
a[0][0] = x1;
a[1][0] = y1;
a[2][0] = z1;
a[0][1] = x2;
a[1][1] = y2;
a[2][1] = z2;
a[0][2] = x3;
a[1][2] = y3;
a[2][2] = z3;
// a[0][0] = x1;
// a[1][0] = x2;
// a[2][0] = x3;
// a[0][1] = y1;
// a[1][1] = y2;
// a[2][1] = y3;
// a[0][2] = z1;
// a[1][2] = z2;
// a[2][2] = z3;
double temp;
for (int i = 0; i < order; i++) {
for (int j = 0; j < 2 * order; j++) {
// Add '1' at the diagonal places of
// the matrix to create a identity matirx
if (j == (i + order))
a[i][j] = 1;
}
}
for (int i = order - 1; i > 0; i--) {
//Swapping each and every element of the two rows
if (a[i - 1][0] < a[i][0])
for (int j = 0; j < 2 * order; j++) {
// Swapping of the row, if above
// condition satisfied.
temp = a[i][j];
a[i][j] = a[i - 1][j];
a[i - 1][j] = temp;
}
}
// Replace a row by sum of itself and a
// constant multiple of another row of the matrix
for (int i = 0; i < order; i++) {
for (int j = 0; j < order; j++) {
if (j != i) {
temp = a[j][i] / a[i][i];
for (int k = 0; k < 2 * order; k++) {
a[j][k] -= a[i][k] * temp;
}
}
}
}
// Multiply each row by a nonzero integer.
// Divide row element by the diagonal element
for (int i = 0; i < order; i++) {
temp = a[i][i];
for (int j = 0; j < 2 * order; j++) {
a[i][j] = a[i][j] / temp;
}
}
// std::cout << "inverted matrix : " << std::endl;
// for (int i = 0; i < n; i++) {
// for (int j = order; j < 6; j++) {
// std::cout << a[i][j] << " ";
// }
// std::cout<<std::endl;
// }
x1 = a[0][0+order];
y1 = a[1][0+order];
z1 = a[2][0+order];
x2 = a[0][1+order];
y2 = a[1][1+order];
z2 = a[2][1+order];
x3 = a[0][2+order];
y3 = a[1][2+order];
z3 = a[2][2+order];
// x1 = a[0][0+order];
// x2 = a[1][0+order];
// x3 = a[2][0+order];
// y1 = a[0][1+order];
// y2 = a[1][1+order];
// y3 = a[2][1+order];
// z1 = a[0][2+order];
// z2 = a[1][2+order];
// z3 = a[2][2+order];
return;
}
}; };
......
...@@ -89,11 +89,12 @@ void Sten_matrix<TYPE2, STENCIL_TYP>::Allocate_data() ...@@ -89,11 +89,12 @@ void Sten_matrix<TYPE2, STENCIL_TYP>::Allocate_data()
if (sten_typ.Not_const_y(id)==false) Ny = 2; if (sten_typ.Not_const_y(id)==false) Ny = 2;
if (sten_typ.Not_const_z(id)==false) Nz = 2; if (sten_typ.Not_const_z(id)==false) Nz = 2;
//cout << "Nx: " << Nx << " Ny: " << Ny << " Nz: " << Nz << endl; //cout << "Hex: "<<id<< "Nx: " << Nx << " Ny: " << Ny << " Nz: " << Nz << endl;
N_total = (Nx-1) * (Ny-1) * (Nz-1) * 27; N_total = (Nx-1) * (Ny-1) * (Nz-1) * 27;
sten_m_hexahedra[id] = new TYPE2[N_total]; sten_m_hexahedra[id] = new TYPE2[N_total];
for (i=0;i<N_total;++i) sten_m_hexahedra[id][i] = (TYPE2)0; for (i=0;i<N_total;++i)
{sten_m_hexahedra[id][i] = (TYPE2)0;}
} else sten_m_hexahedra[id] = NULL; } else sten_m_hexahedra[id] = NULL;
} }
......
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