Commit 51060ea8 authored by Phillip Lino Rall's avatar Phillip Lino Rall
Browse files

core aufgeräumt, updatefunktionen hinzugefügt, minimalbeispiel für jens auch

parents e7b82f34 bb52a776
...@@ -372,6 +372,11 @@ inline std::vector<int> find_p_quad ( dir3D quad, Hexahedron_el *hex ) ...@@ -372,6 +372,11 @@ inline std::vector<int> find_p_quad ( dir3D quad, Hexahedron_el *hex )
// WDed, EDed, WTed, ETed, // WDed, EDed, WTed, ETed,
// SWed, SEed, NWed, NEed }; // SWed, SEed, NWed, NEed };
std::vector<int> ijk = {1 , 1 , 1}; std::vector<int> ijk = {1 , 1 , 1};
// for (int iter = 0 ; iter < 12; iter++)
// {
// std::cout << hex->orientation[iter] << "\n" ;
// }
// std::cout << std::endl;
if (quad == Wdir3D) if (quad == Wdir3D)
{ {
...@@ -550,14 +555,45 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex, ...@@ -550,14 +555,45 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
int quadWaddK = (p_quad_Wdir3D[2] == 1) ? 0 : nEdgeZ; int quadWaddK = (p_quad_Wdir3D[2] == 1) ? 0 : nEdgeZ;
// std::cout <<"Nx " << Give_Nx_quadrangle ( idW ) << " Ny " << Give_Ny_quadrangle ( idW ) << "\n";
// std::cout <<"Nx " << Give_Nx_quadrangle ( idE ) << " Ny " << Give_Ny_quadrangle ( idE ) << "\n";
// std::cout <<"Nx " << Give_Nx_quadrangle ( idS ) << " Ny " << Give_Ny_quadrangle ( idS ) << "\n";
// std::cout <<"Nx " << Give_Nx_quadrangle ( idN ) << " Ny " << Give_Ny_quadrangle ( idN ) << "\n";
// std::cout <<"Nx " << Give_Nx_quadrangle ( idD ) << " Ny " << Give_Ny_quadrangle ( idD ) << "\n";
// std::cout <<"Nx " << Give_Nx_quadrangle ( idT ) << " Ny " << Give_Ny_quadrangle ( idT ) << std::endl;
// std::cout << std::endl;
// std::cout << " transface coordiantes: "<< std::endl;
//with modification //with modification
QPsi[0] = Give_coord_quadrangle( idW , quadWaddJ + j*p_quad_Wdir3D[1] , quadWaddK + k*p_quad_Wdir3D[2]);
QPsi[1] = Give_coord_quadrangle( idE , quadEaddJ + j*p_quad_Edir3D[1] , quadEaddK + k*p_quad_Edir3D[2]);
QPsi[2] = Give_coord_quadrangle( idS , quadSaddI + i*p_quad_Sdir3D[0] , quadSaddK + k*p_quad_Sdir3D[2]);
QPsi[3] = Give_coord_quadrangle( idN , quadNaddI + i*p_quad_Ndir3D[0] , quadNaddK + k*p_quad_Ndir3D[2]);
QPsi[4] = Give_coord_quadrangle( idD , quadDaddI + i*p_quad_Ddir3D[0] , quadDaddJ + j*p_quad_Ddir3D[1]);
QPsi[5] = Give_coord_quadrangle( idT , quadTaddI + i*p_quad_Tdir3D[0] , quadTaddJ + j*p_quad_Tdir3D[1]);
bool invertTD = false;
bool invertNS = false;
bool invertEW = false;
if (!(nEdgeY == Give_Nx_quadrangle ( idW ) && nEdgeZ == Give_Ny_quadrangle ( idW )))
{
invertEW = true;
}
if (!(nEdgeX == Give_Nx_quadrangle ( idS ) && nEdgeZ == Give_Ny_quadrangle ( idS ) ))
{
invertNS = true;
}
if (!(nEdgeX == Give_Nx_quadrangle ( idD ) && nEdgeY == Give_Ny_quadrangle ( idD )))
{
invertTD = true;
}
QPsi[0] = Give_coord_quadrangle( idW , quadWaddJ + j*p_quad_Wdir3D[1] , quadWaddK + k*p_quad_Wdir3D[2],invertEW);
QPsi[1] = Give_coord_quadrangle( idE , quadEaddJ + j*p_quad_Edir3D[1] , quadEaddK + k*p_quad_Edir3D[2],invertEW);
QPsi[2] = Give_coord_quadrangle( idS , quadSaddI + i*p_quad_Sdir3D[0] , quadSaddK + k*p_quad_Sdir3D[2],invertNS);
QPsi[3] = Give_coord_quadrangle( idN , quadNaddI + i*p_quad_Ndir3D[0] , quadNaddK + k*p_quad_Ndir3D[2],invertNS);
QPsi[4] = Give_coord_quadrangle( idD , quadDaddI + i*p_quad_Ddir3D[0] , quadDaddJ + j*p_quad_Ddir3D[1],invertTD);
QPsi[5] = Give_coord_quadrangle( idT , quadTaddI + i*p_quad_Tdir3D[0] , quadTaddJ + j*p_quad_Tdir3D[1],invertTD);
//delete again!!!
// QPsi[5] = Give_coord_quadrangle( idT , quadTaddJ + j*p_quad_Tdir3D[1] , quadTaddI + i*p_quad_Tdir3D[0]);
//std::cout << std::endl;
// Version: 1 x Flaechen + 2 x Kanten + 3 x Ecken // Version: 1 x Flaechen + 2 x Kanten + 3 x Ecken
// ///////////////////////////////////////////////////////// // /////////////////////////////////////////////////////////
Pres = (QPsi[0] + eta * (QPsi[1]-QPsi[0]) + QPsi[2] + xi * (QPsi[3]-QPsi[2]) + QPsi[4] + phi * (QPsi[5]-QPsi[4])); Pres = (QPsi[0] + eta * (QPsi[1]-QPsi[0]) + QPsi[2] + xi * (QPsi[3]-QPsi[2]) + QPsi[4] + phi * (QPsi[5]-QPsi[4]));
...@@ -597,21 +633,21 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex, ...@@ -597,21 +633,21 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
if ( md==0 ) // EW if ( md==0 ) // EW
{ {
PSW = Give_coord_quadrangle(idS,quadSaddI+i*p_quad_Sdir3D[0] ,quadSaddK); PSW = Give_coord_quadrangle(idS,quadSaddI+i*p_quad_Sdir3D[0] ,quadSaddK, invertNS);
if(coordianteFromCorners) if(coordianteFromCorners)
PSW = Psi[SDed]; PSW = Psi[SDed];
PSE = Give_coord_quadrangle(idN,quadNaddI+i * p_quad_Ndir3D[0],quadNaddK); PSE = Give_coord_quadrangle(idN,quadNaddI+i * p_quad_Ndir3D[0],quadNaddK, invertNS);
if(coordianteFromCorners) if(coordianteFromCorners)
PSE = Psi[NDed]; PSE = Psi[NDed];
PNW = Give_coord_quadrangle(idS,quadSaddI+i*p_quad_Sdir3D[0] ,nEdgeZ-quadSaddK); PNW = Give_coord_quadrangle(idS,quadSaddI+i*p_quad_Sdir3D[0] ,nEdgeZ-quadSaddK, invertNS);
if(coordianteFromCorners) if(coordianteFromCorners)
PNW = Psi[STed]; PNW = Psi[STed];
PNE = Give_coord_quadrangle(idN,quadNaddI+i * p_quad_Ndir3D[0],nEdgeZ-quadNaddK); PNE = Give_coord_quadrangle(idN,quadNaddI+i * p_quad_Ndir3D[0],nEdgeZ-quadNaddK, invertNS);
if(coordianteFromCorners) if(coordianteFromCorners)
PNE = Psi[NTed]; PNE = Psi[NTed];
...@@ -625,21 +661,21 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex, ...@@ -625,21 +661,21 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
{ {
PSW = Give_coord_quadrangle(idD,quadDaddI,quadDaddJ+j*p_quad_Ddir3D[1]); PSW = Give_coord_quadrangle(idD,quadDaddI,quadDaddJ+j*p_quad_Ddir3D[1], invertTD);
if(coordianteFromCorners) if(coordianteFromCorners)
PSW = Psi[WDed]; PSW = Psi[WDed];
PSE = Give_coord_quadrangle(idD,nEdgeX-quadDaddI,quadDaddJ+j*p_quad_Ddir3D[1]); PSE = Give_coord_quadrangle(idD,nEdgeX-quadDaddI,quadDaddJ+j*p_quad_Ddir3D[1], invertTD);
if(coordianteFromCorners) if(coordianteFromCorners)
PSE = Psi[EDed]; PSE = Psi[EDed];
PNW = Give_coord_quadrangle(idT,quadTaddI,quadTaddJ+j*p_quad_Tdir3D[1]); PNW = Give_coord_quadrangle(idT,quadTaddI,quadTaddJ+j*p_quad_Tdir3D[1], invertTD);
if(coordianteFromCorners) if(coordianteFromCorners)
PNW = Psi[WTed]; PNW = Psi[WTed];
PNE = Give_coord_quadrangle(idT,nEdgeX-quadTaddI,quadTaddJ+j*p_quad_Tdir3D[1]); PNE = Give_coord_quadrangle(idT,nEdgeX-quadTaddI,quadTaddJ+j*p_quad_Tdir3D[1], invertTD);
if(coordianteFromCorners) if(coordianteFromCorners)
PNE = Psi[ETed]; PNE = Psi[ETed];
...@@ -654,20 +690,20 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex, ...@@ -654,20 +690,20 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
{ {
PSW = Give_coord_quadrangle(idW,quadWaddJ,quadWaddK+k*p_quad_Wdir3D[2]); PSW = Give_coord_quadrangle(idW,quadWaddJ,quadWaddK+k*p_quad_Wdir3D[2], invertEW);
if(coordianteFromCorners) if(coordianteFromCorners)
PSW = Psi[SWed]; PSW = Psi[SWed];
PSE = Give_coord_quadrangle(idE,quadEaddJ,quadEaddK+k*p_quad_Edir3D[2]); PSE = Give_coord_quadrangle(idE,quadEaddJ,quadEaddK+k*p_quad_Edir3D[2], invertEW);
if(coordianteFromCorners) if(coordianteFromCorners)
PSE = Psi[SEed]; PSE = Psi[SEed];
PNW = Give_coord_quadrangle(idW,nEdgeY-quadWaddJ,quadWaddK+k*p_quad_Wdir3D[2]); PNW = Give_coord_quadrangle(idW,nEdgeY-quadWaddJ,quadWaddK+k*p_quad_Wdir3D[2], invertEW);
if(coordianteFromCorners) if(coordianteFromCorners)
PNW = Psi[NWed]; PNW = Psi[NWed];
PNE = Give_coord_quadrangle(idE,nEdgeY-quadEaddJ,quadEaddK+k*p_quad_Edir3D[2]); PNE = Give_coord_quadrangle(idE,nEdgeY-quadEaddJ,quadEaddK+k*p_quad_Edir3D[2], invertEW);
if(coordianteFromCorners) if(coordianteFromCorners)
PNE = Psi[NEed]; PNE = Psi[NEed];
...@@ -764,7 +800,7 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex, ...@@ -764,7 +800,7 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
return Pres + hex->transform_hex ( id_hex,eta,xi,phi ); return Pres + hex->transform_hex ( id_hex,eta,xi,phi );
} }
D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const D3vector Blockgrid::Give_coord_quadrangle (int id, int i, int j, bool invert) const
{ {
if(bg_coord != NULL) if(bg_coord->blockgrid_quad_coordinates_calculated) if(bg_coord != NULL) if(bg_coord->blockgrid_quad_coordinates_calculated)
{ {
...@@ -775,6 +811,8 @@ D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const ...@@ -775,6 +811,8 @@ D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const
} }
double eta, xi; double eta, xi;
D3vector PSW, PSE, PNW, PNE, PW, PE, PN, PS; D3vector PSW, PSE, PNW, PNE, PW, PE, PN, PS;
D3vector PTransFace, PWithoutTrans; D3vector PTransFace, PWithoutTrans;
...@@ -783,9 +821,30 @@ D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const ...@@ -783,9 +821,30 @@ D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const
quad = this->ug->Give_quadrangle ( id ); quad = this->ug->Give_quadrangle ( id );
double Nx = (double) Give_Nx_quadrangle ( id );
double Ny = ( double ) Give_Ny_quadrangle ( id );
// std::cout << "nx ny " << Nx << " " << Ny << std::endl;
// if (id == 23 && Nx != Ny)
// {
// std::cout << "break";
// }
// std::cout << "nx ny " << Nx << " " << Ny << std::endl;
if (invert)
{
double temp = Nx; Nx = Ny; Ny = temp;
}
//std::cout << "nx ny " << Nx << " " << Ny << std::endl;
eta = ( double ) i/ Nx;
xi = ( double ) j/ Ny;
if (eta > 1 || xi > 1)
{
std::cout << "debug : flag is " << invert << ", but should be " << !invert << std::endl;
std::cout << std::endl;
}
eta = ( double ) i/ ( double ) Give_Nx_quadrangle ( id );
xi = ( double ) j/ ( double ) Give_Ny_quadrangle ( id );
// @todo : so sollte es sein // @todo : so sollte es sein
PSW = quad->Give_coord ( SWdir2D ); PSW = quad->Give_coord ( SWdir2D );
...@@ -814,7 +873,9 @@ D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const ...@@ -814,7 +873,9 @@ D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const
if (quad->transform != NULL) if (quad->transform != NULL)
{ {
PTransFace = quad->transform ( eta, xi, quad->shiftPointer(ug->Get_pointer_global_data())); PTransFace = quad->transform ( eta, xi, quad->shiftPointer(ug->Get_pointer_global_data()));
PTransFace = PSW + (PSE-PSW) * eta + (PNW-PSW)*xi + (PNE-PSE-PNW+PSW)*xi*eta + PTransFace; //std::cout << "PTransFace ";PTransFace.Print();
D3vector TEMP = PSW + (PSE-PSW) * eta + (PNW-PSW)*xi + (PNE-PSE-PNW+PSW)*xi*eta ;
PTransFace = TEMP + PTransFace;// + PTransFace;
return PTransFace; return PTransFace;
} }
...@@ -849,7 +910,6 @@ D3vector Blockgrid::Give_coord_edge ( int id, int i ) const ...@@ -849,7 +910,6 @@ D3vector Blockgrid::Give_coord_edge ( int id, int i ) const
if(bg_coord != NULL) if(bg_coord->blockgrid_edge_coordinates_calculated) if(bg_coord != NULL) if(bg_coord->blockgrid_edge_coordinates_calculated)
{ {
return bg_coord->blockgrid_edge_coordinates.at(id).at(i) ; return bg_coord->blockgrid_edge_coordinates.at(id).at(i) ;
} }
...@@ -903,49 +963,85 @@ D3vector Blockgrid::Give_coord_edge ( int id, int i ) const ...@@ -903,49 +963,85 @@ D3vector Blockgrid::Give_coord_edge ( int id, int i ) const
} }
} }
} }
ug->Set_edge_to_hex_id(id, idHex);
ug->Set_edge_to_quad_id(id,quadId); ug->Set_edge_to_quad_id(id,quadId);
ug->Set_edge_to_quad_dir(id,(int)dirForQuad); ug->Set_edge_to_quad_dir(id,(int)dirForQuad);
} }
else else
{ {
idHex = ug->Give_edge_to_hex_id(id);
quadId = ug->Give_edge_to_quad_id(id); quadId = ug->Give_edge_to_quad_id(id);
dirForQuad = (dir3D)ug->Give_edge_to_quad_dir(id); dirForQuad = (dir3D)ug->Give_edge_to_quad_dir(id);
quad = ug->Give_quadrangle(quadId); quad = ug->Give_quadrangle(quadId);
} }
int Nx = Give_Nx_quadrangle(dirForQuad); int Nx = Give_Nx_quadrangle(quadId);
int Ny = Give_Ny_quadrangle(dirForQuad); int Ny = Give_Ny_quadrangle(quadId);
int nEdgeX = Give_Nx_hexahedron ( idHex );
int nEdgeY = Give_Ny_hexahedron ( idHex );
int nEdgeZ = Give_Nz_hexahedron ( idHex );
bool invert = false;
// std::cout << "dirForQuad " << dirForQuad << std::endl;
// std::cout << "Nx " << Nx << " Ny " << Ny << std::endl;
// if (!(Nx == Give_Nx_quadrangle ( dirForQuad ) && Ny == Give_Ny_quadrangle ( dirForQuad ) ))
// {
// invert = true;
// }
// if (dirForQuad == Edir3D || dirForQuad == Wdir3D)
// {
// if (!(nEdgeY == Give_Nx_quadrangle ( dirForQuad ) && nEdgeZ == Give_Ny_quadrangle ( dirForQuad ) ))
// {
// invert = true;
// }
// }
// else if (dirForQuad == Ndir3D || dirForQuad == Sdir3D)
// {
// if (!(nEdgeX == Give_Nx_quadrangle ( dirForQuad ) && nEdgeZ == Give_Ny_quadrangle ( dirForQuad ) ))
// {
// invert = true;
// }
// }
// else if (dirForQuad == Tdir3D || dirForQuad == Ddir3D)
// {
// if (!(nEdgeX == Give_Nx_quadrangle ( dirForQuad ) && nEdgeY == Give_Ny_quadrangle ( dirForQuad )))
// {
// invert = true;
// }
// }
// else
// {
// std::cout << "should not happen, dirForQuad only has 6 possibilities\n";
// }
if (Nx == nEdgeX && Ny == nEdgeY || Nx == nEdgeX && Ny == nEdgeZ || Nx == nEdgeY && Ny == nEdgeZ )
{
invert = false;
}
else
{
invert = true;
}
invert = false;
if (id == quad->Give_id_edge(Wdir2D)) //Wdir2D if (id == quad->Give_id_edge(Wdir2D)) //Wdir2D
{ {
Psi = Give_coord_quadrangle(quadId,0,i,invert);
Psi = Give_coord_quadrangle(quadId,0,i);
} }
if (id == quad->Give_id_edge(Edir2D)) //Edir2D if (id == quad->Give_id_edge(Edir2D)) //Edir2D
{ {
Psi = Give_coord_quadrangle(quadId,Nx,i,invert);
Psi = Give_coord_quadrangle(quadId,Nx,i);
} }
if (id == quad->Give_id_edge(Sdir2D)) //Sdir2D if (id == quad->Give_id_edge(Sdir2D)) //Sdir2D
{ {
Psi = Give_coord_quadrangle(quadId,i,0,invert);
Psi = Give_coord_quadrangle(quadId,i,0);
} }
if (id == quad->Give_id_edge(Ndir2D)) //Ndir2D if (id == quad->Give_id_edge(Ndir2D)) //Ndir2D
{ {
Psi = Give_coord_quadrangle(quadId,i,Ny,invert);
Psi = Give_coord_quadrangle(quadId,i,Ny);
} }
return Psi; return Psi;
...@@ -1066,6 +1162,8 @@ void Blockgrid_coordinates::init_blockgrid_coordinates() ...@@ -1066,6 +1162,8 @@ void Blockgrid_coordinates::init_blockgrid_coordinates()
} }
} }
//test:: ??? Nx = 3????
blockgrid_edge_coordinates.resize(bg->Give_unstructured_grid()->Give_number_edges()); blockgrid_edge_coordinates.resize(bg->Give_unstructured_grid()->Give_number_edges());
for (int id_edge = 0 ; id_edge < bg->Give_unstructured_grid()->Give_number_edges() ; id_edge++) for (int id_edge = 0 ; id_edge < bg->Give_unstructured_grid()->Give_number_edges() ; id_edge++)
{ {
......
...@@ -114,7 +114,7 @@ class Blockgrid { ...@@ -114,7 +114,7 @@ class Blockgrid {
void Set_blockgrid_coordinates(Blockgrid_coordinates* bg_coord_ ) { bg_coord = bg_coord_; }; void Set_blockgrid_coordinates(Blockgrid_coordinates* bg_coord_ ) { bg_coord = bg_coord_; };
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_quadrangle(int id, int i, int j, bool invert = false) const;
D3vector Give_coord_edge(int id, int i) const; D3vector Give_coord_edge(int id, int i) const;
D3vector Give_coord_point(int id) const; D3vector Give_coord_point(int id) const;
......
...@@ -450,6 +450,7 @@ D3vector transform_left_lens_diag_NE_quad ( double t1, double t2, double* pointe ...@@ -450,6 +450,7 @@ D3vector transform_left_lens_diag_NE_quad ( double t1, double t2, double* pointe
} }
D3vector transform_diag_inner_faces_NE_quad( double t1, double t2, double* pointer_global_data) { D3vector transform_diag_inner_faces_NE_quad( double t1, double t2, double* pointer_global_data) {
return D3vector{0,0,0};
double R_global_data = pointer_global_data[0]; double R_global_data = pointer_global_data[0];
double r_global_data = pointer_global_data[1]; double r_global_data = pointer_global_data[1];
double curvatureLeft_global_data = pointer_global_data[2]; double curvatureLeft_global_data = pointer_global_data[2];
...@@ -507,6 +508,7 @@ D3vector transform_diag_inner_faces_NE_quad( double t1, double t2, double* point ...@@ -507,6 +508,7 @@ D3vector transform_diag_inner_faces_NE_quad( double t1, double t2, double* point
} }
D3vector transform_diag_inner_faces_NE_quad_cut( double t1, double t2, double* pointer_global_data) { D3vector transform_diag_inner_faces_NE_quad_cut( double t1, double t2, double* pointer_global_data) {
return D3vector{0,0,0};
double R_global_data = pointer_global_data[0]; double R_global_data = pointer_global_data[0];
double r_global_data = pointer_global_data[1]; double r_global_data = pointer_global_data[1];
double curvatureLeft_global_data = pointer_global_data[2]; double curvatureLeft_global_data = pointer_global_data[2];
...@@ -562,6 +564,7 @@ D3vector transform_diag_inner_faces_NE_quad_cut( double t1, double t2, double* p ...@@ -562,6 +564,7 @@ D3vector transform_diag_inner_faces_NE_quad_cut( double t1, double t2, double* p
zLeft = 0; zLeft = 0;
} }
double z = zRight * (t1) + zLeft * (1-t1) ; double z = zRight * (t1) + zLeft * (1-t1) ;
if (( t1 == 0 && t2 == 0)||( t1 == 1 && t2 == 0)||( t1 == 0 && t2 == 1)||( t1 == 1 && t2 == 1)) if (( t1 == 0 && t2 == 0)||( t1 == 1 && t2 == 0)||( t1 == 0 && t2 == 1)||( t1 == 1 && t2 == 1))
{ {
...@@ -658,8 +661,12 @@ D3vector transform_right_lens_diag_NE_quad ( double t1, double t2, double* point ...@@ -658,8 +661,12 @@ D3vector transform_right_lens_diag_NE_quad ( double t1, double t2, double* point
radiusSquared) radiusSquared)
+ sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data))))); + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data)))));
//std::cout << "z inner " << z << std::endl;
if (std::isnan(x) || std::isnan(y) || std::isnan(z)) if (std::isnan(x) || std::isnan(y) || std::isnan(z))
{ z = 0;} { z = 0;}
// std::cout << "transform_right_lens_diag_NE_quad ";
// D3vector ( x,y,z).Print();
// std::cout << std::endl;
return D3vector ( x,y,z); return D3vector ( x,y,z);
} }
...@@ -1060,13 +1067,15 @@ D3vector transform_left_lens_diag_NE_quad_cut ( double t1, double t2, double* po ...@@ -1060,13 +1067,15 @@ D3vector transform_left_lens_diag_NE_quad_cut ( double t1, double t2, double* po
x = x*t2; x = x*t2;
double y = R_global_data * (cos ( t1*0.5*M_PI )- ( 1-t1 )) ; double y = R_global_data * (cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
y = y*t2; y = y*t2;
double xT = x;
double yT = y;
//added due to bended inner part : //added due to bended inner part :
double xAdd = r_global_data * ( sin ( t1*0.5*M_PI )-t1); double xAdd = r_global_data * ( sin ( t1*0.5*M_PI )-t1);
double yAdd = r_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ; double yAdd = r_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
x += (1-t2) * xAdd; x += (1-t2) * xAdd;
y += (1-t2) * yAdd; y += (1-t2) * yAdd;
double actualX = 0 + r_global_data * t1 + 0 * t2 + (R_global_data-r_global_data) * t1 * t2; double actualX = 0 + r_global_data * t1 + 0 * t2 + (R_global_data-r_global_data) * t1 * t2;
double actualY = r_global_data - r_global_data * t1 + (R_global_data-r_global_data) * t2 - (R_global_data-r_global_data) * t1 * t2; double actualY = r_global_data - r_global_data * t1 + (R_global_data-r_global_data) * t2 - (R_global_data-r_global_data) * t1 * t2;
...@@ -1079,12 +1088,14 @@ D3vector transform_left_lens_diag_NE_quad_cut ( double t1, double t2, double* po ...@@ -1079,12 +1088,14 @@ D3vector transform_left_lens_diag_NE_quad_cut ( double t1, double t2, double* po
if (z_left_inner_global_data == z_left_outer_global_data) if (z_left_inner_global_data == z_left_outer_global_data)
{ {
//std::cout << "set to zero again!\n";
z = 0; z = 0;
} }
if (std::isnan(z)) if (std::isnan(z))
{z = 0;} {z = 0;}
//return D3vector ( xT,yT,z);
return D3vector ( x,y,z); return D3vector ( x,y,z);
} }
...@@ -1121,6 +1132,8 @@ D3vector transform_right_lens_diag_NE_quad_cut ( double t1, double t2, double* p ...@@ -1121,6 +1132,8 @@ D3vector transform_right_lens_diag_NE_quad_cut ( double t1, double t2, double* p
double y = R_global_data * (cos ( t1*0.5*M_PI )- ( 1-t1 )) ; double y = R_global_data * (cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
y = y*t2; y = y*t2;
double xT = x;
double yT = y;
//added due to bended inner part : //added due to bended inner part :
double xAdd = r_global_data * ( sin ( t1*0.5*M_PI )-t1); double xAdd = r_global_data * ( sin ( t1*0.5*M_PI )-t1);
double yAdd = r_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ; double yAdd = r_global_data *(cos ( t1*0.5*M_PI )- ( 1-t1 )) ;
...@@ -1133,11 +1146,18 @@ D3vector transform_right_lens_diag_NE_quad_cut ( double t1, double t2, double* p ...@@ -1133,11 +1146,18 @@ D3vector transform_right_lens_diag_NE_quad_cut ( double t1, double t2, double* p
double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY)); double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));
// std::cout << "(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data))) " << (-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data))) << "\n";
// std::cout << "t1,2 :: " << t1 << " " << t2 << "\n";
//std::cout << "radius " << sqrt(radiusSquared) << std::endl;
double z = offsetZ_global_data+ double z = offsetZ_global_data+
sign*( ( sqrt(pow(curvatureRight_global_data,2)- sign*( ( sqrt(pow(curvatureRight_global_data,2)-
radiusSquared) radiusSquared)
+ sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data))))); + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data)))));
// z =radiusSquared;
//z *= 10;
//z = 1;
//std::cout << "z outer" << z << std::endl;
// double ADAD = offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-pow(r_global_data * (1-t) + R_global_data * t,2)) +sign *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * z_left_outer_global_data))));
if (z_right_inner_global_data == z_right_outer_global_data) if (z_right_inner_global_data == z_right_outer_global_data)
{ {
...@@ -1146,7 +1166,19 @@ D3vector transform_right_lens_diag_NE_quad_cut ( double t1, double t2, double* p ...@@ -1146,7 +1166,19 @@ D3vector transform_right_lens_diag_NE_quad_cut ( double t1, double t2, double* p
if (std::isnan(x) || std::isnan(y) || std::isnan(z)) if (std::isnan(x) || std::isnan(y) || std::isnan(z))
{ z = 0;} { z = 0;}