From 655cb7daed5b5d25a160b777acc03d7e6390a660 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Wed, 5 Jun 2013 14:56:34 +0200 Subject: [PATCH] operator: fix in handling curve primitives in cylindrical coordinates Signed-off-by: Thorsten Liebig --- FDTD/operator.cpp | 27 +++++++- FDTD/operator.h | 11 +++ FDTD/operator_cylinder.cpp | 138 +++++++++++++++++++++++++++++++++++++ FDTD/operator_cylinder.h | 4 ++ 4 files changed, 177 insertions(+), 3 deletions(-) diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index f9cdff1..4e21d4f 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -348,14 +348,35 @@ int Operator::SnapBox2Mesh(const double* start, const double* stop, unsigned int return -1; } +int Operator::SnapLine2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh, bool fullMesh) const +{ + bool bStartIn[3]; + bool bStopIn[3]; + SnapToMesh(start, uiStart, dualMesh, fullMesh, bStartIn); + SnapToMesh(stop, uiStop, dualMesh, fullMesh, bStopIn); + + int ret = 0; + if (!(bStartIn[0] && bStartIn[1] && bStartIn[2])) + ret = ret + 1; + if (!(bStopIn[0] && bStopIn[1] && bStopIn[2])) + ret = ret + 2; + if (ret==0) + return ret; + + //fixme, do we need to do something about start or stop being outside the field domain? + //maybe caclulate the intersection point and snap to that? + //it seems to work like this as well... + + return ret; +} + + struct Operator::Grid_Path Operator::FindPath(double start[], double stop[]) { struct Grid_Path path; unsigned int uiStart[3],uiStop[3],currPos[3]; - SnapToMesh(start,uiStart,false,true); - SnapToMesh(stop,uiStop,false,true); - + SnapLine2Mesh(start, stop, uiStart, uiStop, false, true); currPos[0]=uiStart[0]; currPos[1]=uiStart[1]; currPos[2]=uiStart[2]; diff --git a/FDTD/operator.h b/FDTD/operator.h index f690a44..ef7be56 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -140,6 +140,17 @@ public: //! Snap a given box to the FDTD mesh virtual int SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, bool fullMesh=false, int SnapMethod=0, bool* bStartIn=NULL, bool* bStopIn=NULL) const; + //! Snap a given line to the operator mesh + /*! + \param[in] start coorindate of the line + \param[in] stop coorindate of the line + \param[out] uiStart the snapped line-start coorindate index + \param[out] uiStop the snapped line-stop coorindate index + \param[in] dualMesh snap to main or dual mesh (default is main mesh) + \return returns a status, 0 = success, 1 = start outside, 2 = stop outside, 3 = both outside + */ + virtual int SnapLine2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, bool fullMesh=false) const; + virtual void AddExtension(Operator_Extension* op_ext); virtual void DeleteExtension(Operator_Extension* op_ext); virtual size_t GetNumberOfExtentions() const {return m_Op_exts.size();} diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp index 58c437a..0c2bda7 100644 --- a/FDTD/operator_cylinder.cpp +++ b/FDTD/operator_cylinder.cpp @@ -21,6 +21,7 @@ #include "operator_cylinder.h" #include "extensions/operator_extension.h" #include "extensions/operator_ext_cylinder.h" +#include "tools/useful.h" Operator_Cylinder* Operator_Cylinder::New(unsigned int numThreads) { @@ -260,6 +261,143 @@ int Operator_Cylinder::SnapBox2Mesh(const double* start, const double* stop, uns return Operator_Multithread::SnapBox2Mesh(l_start, l_stop, uiStart, uiStop, dualMesh, fullMesh, SnapMethod, bStartIn, bStopIn); } +int Operator_Cylinder::SnapLine2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh, bool fullMesh) const +{ + int ret = Operator_Multithread::SnapLine2Mesh(start, stop, uiStart, uiStop, dualMesh, fullMesh); + + if ((stop[1]>start[1]) && (uiStop[1]uiStart[1]) && (uiStop[1]==GetNumberOfLines(1, fullMesh)-1-(int)CC_closedAlpha)) + uiStop[1] = 0; + + return ret; +} + + +struct Operator::Grid_Path Operator_Cylinder::FindPath(double start[], double stop[]) +{ + double l_start[3]; + double l_stop[3]; + if (stop[1]= a_start) + return Operator_Multithread::FindPath(start, stop); + + if (g_settings.GetVerboseLevel()>2) + { + cerr << __func__ << ": A path was leaving the alpha-direction mesh..." << endl; + } + + // this section comes into play, if the line moves over the angulare mesh-end/start + // we try to have one part of the path on both "ends" of the mesh and stitch them together + + struct Grid_Path path; + struct Grid_Path path1; + struct Grid_Path path2; + + // calculate the intersection of the line with the a-max boundary + double p0[3],p1[3],p2[3]; + for (int n=0;n<3;++n) + { + p0[n] = GetDiscLine(n,0); + p1[n] = p0[n]; + p2[n] = p0[n]; + } + p0[1] = GetDiscLine(1,GetNumberOfLines(1,true)-1-(int)CC_closedAlpha); + p1[1] = p0[1]; + p2[1] = p0[1]; + p1[0] = discLines[0][numLines[0]-1]; + p2[2] = discLines[2][numLines[2]-1]; + + TransformCoordSystem(p0,p0,m_MeshType,CARTESIAN); + TransformCoordSystem(p1,p1,m_MeshType,CARTESIAN); + TransformCoordSystem(p2,p2,m_MeshType,CARTESIAN); + + double c_start[3],c_stop[3]; + TransformCoordSystem(l_start,c_start,m_MeshType,CARTESIAN); + TransformCoordSystem(l_stop,c_stop,m_MeshType,CARTESIAN); + double intersect[3]; + double dist; + int ret = LinePlaneIntersection(p0,p1,p2,c_start,c_stop,intersect,dist); + if (ret<0) + { + cerr << __func__ << ": Error, unable to calculate intersection, this should not happen!" << endl; + return path; // return empty path; + } + + if (ret==0) + { + TransformCoordSystem(intersect,intersect,CARTESIAN,m_MeshType); + path1 = Operator::FindPath(l_start, intersect); + if (g_settings.GetVerboseLevel()>2) + cerr << __func__ << ": Intersection top: " << intersect[0] << "," << intersect[1] << "," << intersect[2] << endl; + } //otherwise the path was not intersecting the upper a-bound... + + if (CC_closedAlpha==false) + { + for (int n=0;n<3;++n) + { + p0[n] = GetDiscLine(n,0); + p1[n] = p0[n]; + p2[n] = p0[n]; + } + p1[0] = discLines[0][numLines[0]-1]; + p2[2] = discLines[2][numLines[2]-1]; + + TransformCoordSystem(p0,p0,m_MeshType,CARTESIAN); + TransformCoordSystem(p1,p1,m_MeshType,CARTESIAN); + TransformCoordSystem(p2,p2,m_MeshType,CARTESIAN); + + TransformCoordSystem(l_start,c_start,m_MeshType,CARTESIAN); + TransformCoordSystem(l_stop,c_stop,m_MeshType,CARTESIAN); + + ret = LinePlaneIntersection(p0,p1,p2,c_start,c_stop,intersect,dist); + TransformCoordSystem(intersect,intersect,CARTESIAN,m_MeshType); + } + + if (ret==0) + { + path2 = Operator::FindPath(intersect, l_stop); + if (g_settings.GetVerboseLevel()>2) + cerr << __func__ << ": Intersection bottom: " << intersect[0] << "," << intersect[1] << "," << intersect[2] << endl; + } + + //combine path + for (size_t t=0; t