Commit 33c650ce authored by Alexandre Halbach's avatar Alexandre Halbach

Check geo

parent 7e233812
......@@ -37,7 +37,7 @@ double geotools::getdistance(std::vector<double> pt1coords, std::vector<double>
return std::sqrt( std::pow(pt1coords[0]-pt2coords[0],2) + std::pow(pt1coords[1]-pt2coords[1],2) + std::pow(pt1coords[2]-pt2coords[2],2) );
}
std::vector<double> geotools::getplaneangles(std::vector<double> p1, std::vector<double> p2, std::vector<double> p3)
double geotools::getplanerotation(std::string xy, std::vector<double> p1, std::vector<double> p2, std::vector<double> p3)
{
double pi = 3.1415926535897932384;
......@@ -52,7 +52,7 @@ std::vector<double> geotools::getplaneangles(std::vector<double> p1, std::vector
// Shortcut if all points are in a plane parallel to the xy plane.
// This allows colinear points in 2D.
if (std::abs(v12[2])/v12norm < threshold && std::abs(v13[2])/v12norm < threshold)
return std::vector<double>{0.0,0.0};
return 0.0;
// Normal [a b c] to the plane is the cross product:
std::vector<double> vnormal = {v12[1]*v13[2]-v12[2]*v13[1], v12[2]*v13[0]-v12[0]*v13[2], v12[0]*v13[1]-v12[1]*v13[0]};
......@@ -60,14 +60,17 @@ std::vector<double> geotools::getplaneangles(std::vector<double> p1, std::vector
// If in the xz plane:
if (std::abs(vnormal[0])/normalnorm < threshold && std::abs(vnormal[2])/normalnorm < threshold)
return std::vector<double>{0.0,90.0};
return 90;
// If in the yz plane:
if (std::abs(vnormal[1])/normalnorm < threshold && std::abs(vnormal[2])/normalnorm < threshold)
return std::vector<double>{90.0,0.0};
return 90;
// Plane equation is ax + by + cz = d.
// Setting x then y to zero gives the angles we need:
return std::vector<double>{360/2/pi*std::atan(-vnormal[0]/vnormal[2]), 360/2/pi*std::atan(-vnormal[1]/vnormal[2])};
if (xy == "xrot")
return -360/2/pi*std::atan(-vnormal[1]/vnormal[2]);
if (xy == "yrot")
return 360/2/pi*std::atan(-vnormal[0]/vnormal[2]);
}
void geotools::rotate(double alphax, double alphay, double alphaz, std::vector<double>* coords)
......
......@@ -29,8 +29,10 @@ namespace geotools
// Get the distance between two points:
double getdistance(std::vector<double> pt1coords, std::vector<double> pt2coords);
// Get the angle (in degrees) with the x and y axis of the plane defined by the 3 input points
std::vector<double> getplaneangles(std::vector<double> p1, std::vector<double> p2, std::vector<double> p3);
// Get the angle (in degrees) which rotates the plane defined by the 3 input points
// around the x (or y) axis to bring it parallel to the y axis (or x axis).
// Get the angle around the x axis with "xrot" and around the y axis with "yrot"
double getplanerotation(std::string xy, std::vector<double> p1, std::vector<double> p2, std::vector<double> p3);
// Rotate a vector of coordinates by alphax, alphay and alphaz degrees around the x, y and z axis respectively:
void rotate(double alphax, double alphay, double alphaz, std::vector<double>* coords);
......
......@@ -123,13 +123,23 @@ void rawarc::mesh(void)
///// Rotate the arc plane to have it parallel to the xy plane:
std::vector<double> xyrot = geotools::getplaneangles(pccoords, p1coords, p2coords);
double xaxisrot = geotools::getplanerotation("xrot", pccoords, p1coords, p2coords);
if (xaxisrot != 0)
{
sons[0]->rotate(xaxisrot,0,0);
sons[1]->rotate(xaxisrot,0,0);
mycenterpoint->rotate(xaxisrot,0,0);
if (xyrot[0] != 0 || xyrot[1] != 0)
p1coords = *(sons[0]->getcoords());
p2coords = *(sons[1]->getcoords());
pccoords = *(mycenterpoint->getcoords());
}
double yaxisrot = geotools::getplanerotation("yrot", pccoords, p1coords, p2coords);
if (yaxisrot != 0)
{
sons[0]->rotate(-xyrot[1],xyrot[0],0);
sons[1]->rotate(-xyrot[1],xyrot[0],0);
mycenterpoint->rotate(-xyrot[1],xyrot[0],0);
sons[0]->rotate(0,yaxisrot,0);
sons[1]->rotate(0,yaxisrot,0);
mycenterpoint->rotate(0,yaxisrot,0);
p1coords = *(sons[0]->getcoords());
p2coords = *(sons[1]->getcoords());
......@@ -166,7 +176,10 @@ void rawarc::mesh(void)
///// Rotate everything back:
if (xyrot[0] != 0 || xyrot[1] != 0)
rotate(xyrot[1],-xyrot[0],0);
if (xaxisrot != 0 || yaxisrot != 0)
{
rotate(0,-yaxisrot,0);
rotate(-xaxisrot,0,0);
}
}
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