From 7427a72994ad8f71b3d4c042a97cb09552831bed Mon Sep 17 00:00:00 2001 From: Christophe Geuzaine <cgeuzaine@ulg.ac.be> Date: Thu, 25 Oct 2012 20:09:37 +0000 Subject: [PATCH] fix tolerances in isInside --- Geo/MElementOctree.cpp | 3 ++- Geo/MPoint.h | 3 +++ Geo/MQuadrangle.h | 3 ++- Geo/MTriangle.h | 2 +- Plugin/Integrate.cpp | 9 ++++----- Post/shapeFunctions.cpp | 3 +-- Post/shapeFunctions.h | 28 +++++++++++++++------------- demos/indheat.geo | 5 ++++- 8 files changed, 32 insertions(+), 24 deletions(-) diff --git a/Geo/MElementOctree.cpp b/Geo/MElementOctree.cpp index 2c7f841ff7..beca98ee89 100644 --- a/Geo/MElementOctree.cpp +++ b/Geo/MElementOctree.cpp @@ -123,7 +123,8 @@ MElementOctree::~MElementOctree() } -std::vector<MElement *> MElementOctree::findAll(double x, double y, double z, int dim, bool strict) +std::vector<MElement *> MElementOctree::findAll(double x, double y, double z, + int dim, bool strict) { double maxTol = 1.; double tolIncr = 10.; diff --git a/Geo/MPoint.h b/Geo/MPoint.h index 1b09436051..17eb4e6562 100644 --- a/Geo/MPoint.h +++ b/Geo/MPoint.h @@ -68,6 +68,9 @@ class MPoint : public MElement { } virtual bool isInside(double u, double v, double w) { + double tol = _isInsideTolerance; + if(fabs(u) > tol || fabs(v) > tol || fabs(w) > tol) + return false; return true; } virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h index 65d1fde7cd..0f8df873e5 100644 --- a/Geo/MQuadrangle.h +++ b/Geo/MQuadrangle.h @@ -140,7 +140,8 @@ class MQuadrangle : public MElement { virtual bool isInside(double u, double v, double w) { double tol = _isInsideTolerance; - if(u < -(1. + tol) || v < -(1. + tol) || u > (1. + tol) || v > (1. + tol) || fabs(w) > tol) + if(u < -(1. + tol) || v < -(1. + tol) || u > (1. + tol) || v > (1. + tol) || + fabs(w) > tol) return false; return true; } diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h index 54a298904e..963a054aa9 100644 --- a/Geo/MTriangle.h +++ b/Geo/MTriangle.h @@ -143,7 +143,7 @@ class MTriangle : public MElement { virtual bool isInside(double u, double v, double w) { double tol = _isInsideTolerance; - if(u < (-tol) || v < (-tol) || u > ((1. + tol) - v)) + if(u < (-tol) || v < (-tol) || u > ((1. + tol) - v) || fabs(w) > tol) return false; return true; } diff --git a/Plugin/Integrate.cpp b/Plugin/Integrate.cpp index e4a2f1d222..eafc71792f 100644 --- a/Plugin/Integrate.cpp +++ b/Plugin/Integrate.cpp @@ -28,7 +28,8 @@ std::string GMSH_IntegratePlugin::getHelp() const "line/surface elements.\n\n" "If `View' < 0, the plugin is run on the current view.\n\n" "Plugin(Integrate) creates one new view." - "If `OverTime' = 1 , the plugin integrates the scalar view over time instead of over space.\n\n" + "If `OverTime' = 1 , the plugin integrates the scalar view " + "over time instead of over space.\n\n" "Plugin(Integrate) creates one new view."; } @@ -118,9 +119,7 @@ PView *GMSH_IntegratePlugin::execute(PView * v) } else{ - printf("EMI INTEGRATE OVER TIME \n"); - - int timeBeg = data1->getFirstNonEmptyTimeStep(); + int timeBeg = data1->getFirstNonEmptyTimeStep(); int timeEnd = data1->getNumTimeSteps(); for(int ent = 0; ent < data1->getNumEntities(timeBeg); ent++){ for(int ele = 0; ele < data1->getNumElements(timeBeg, ent); ele++){ @@ -147,7 +146,7 @@ PView *GMSH_IntegratePlugin::execute(PView * v) for(int nod = 0; nod < numNodes; nod++){ double val; data1->getValue(step, ent, ele, nod, 0, val); - timeIntegral[nod] += val*dt; + timeIntegral[nod] += val * dt; } } for(int nod = 0; nod < numNodes; nod++) diff --git a/Post/shapeFunctions.cpp b/Post/shapeFunctions.cpp index e027825957..06aa7c3d54 100644 --- a/Post/shapeFunctions.cpp +++ b/Post/shapeFunctions.cpp @@ -5,5 +5,4 @@ #include "shapeFunctions.h" -double element::ONE = 1. + 1.e-6; -double element::ZERO = -1.e-6; +double element::TOL = 1.e-6; diff --git a/Post/shapeFunctions.h b/Post/shapeFunctions.h index 1ebb8d0837..2d73c4d9d8 100644 --- a/Post/shapeFunctions.h +++ b/Post/shapeFunctions.h @@ -15,7 +15,7 @@ class element{ protected: bool _ownData; double *_x, *_y, *_z; - static double ONE, ZERO; + static double TOL; public: element(double *x, double *y, double *z, int numNodes=0) { @@ -50,8 +50,8 @@ public: y = _y[num]; z = _z[num]; } - static void setTolerance (const double tol){ ONE = 1. + tol; ZERO = -tol; } - static double getTolerance () { return -ZERO; } + static void setTolerance (const double tol){ TOL = tol; } + static double getTolerance () { return TOL; } virtual int getDimension() = 0; virtual int getNumNodes() = 0; virtual void getNode(int num, double &u, double &v, double &w) = 0; @@ -305,6 +305,8 @@ public: } int isInside(double u, double v, double w) { + if(fabs(u) > TOL || fabs(v) > TOL || fabs(w) > TOL) + return 0; return 1; } }; @@ -364,7 +366,7 @@ public: } int isInside(double u, double v, double w) { - if(u < -ONE || u > ONE) + if(u < -(1. + TOL) || u > (1. + TOL) || fabs(v) > TOL || fabs(w) > TOL) return 0; return 1; } @@ -460,7 +462,7 @@ public: #endif int isInside(double u, double v, double w) { - if(u < ZERO || v < ZERO || u > (ONE - v)) + if(u < -TOL || v < -TOL || u > ((1. + TOL) - v) || fabs(w) > TOL) return 0; return 1; } @@ -541,7 +543,8 @@ public: } int isInside(double u, double v, double w) { - if(u < -ONE || v < -ONE || u > ONE || v > ONE) + if(u < -(1. + TOL) || v < -(1. + TOL) || u > (1. + TOL) || v > (1. + TOL) || + fabs(w) > TOL) return 0; return 1; } @@ -628,7 +631,7 @@ public: } int isInside(double u, double v, double w) { - if(u < ZERO || v < ZERO || w < ZERO || u > (ONE - v - w)) + if(u < (-TOL) || v < (-TOL) || w < (-TOL) || u > ((1. + TOL) - v - w)) return 0; return 1; } @@ -735,7 +738,8 @@ public: } int isInside(double u, double v, double w) { - if(u < -ONE || v < -ONE || w < -ONE || u > ONE || v > ONE || w > ONE) + if(u < -(1. + TOL) || v < -(1. + TOL) || w < -(1. + TOL) || + u > (1. + TOL) || v > (1. + TOL) || w > (1. + TOL)) return 0; return 1; } @@ -829,7 +833,7 @@ public: } int isInside(double u, double v, double w) { - if(w > ONE || w < -ONE || u < ZERO || v < ZERO || u > (ONE - v)) + if(w > (1. + TOL) || w < -(1. + TOL) || u < (-TOL) || v < (-TOL) || u > ((1. + TOL) - v)) return 0; return 1; } @@ -951,8 +955,8 @@ public: } int isInside(double u, double v, double w) { - if(u < (w - ONE) || u > (ONE - w) || v < (w - ONE) || v > (ONE - w) || - w < ZERO || w > ONE) + if(u < (w - (1. + TOL)) || u > ((1. + TOL) - w) || v < (w - (1. + TOL)) || + v > ((1. + TOL) - w) || w < (-TOL) || w > (1. + TOL)) return 0; return 1; } @@ -981,8 +985,6 @@ class elementFactory{ } }; -#undef ONE -#undef ZERO #undef SQU #endif diff --git a/demos/indheat.geo b/demos/indheat.geo index a587b4e5f3..d1163445ce 100644 --- a/demos/indheat.geo +++ b/demos/indheat.geo @@ -3,7 +3,7 @@ nn = 40; // mesh subdivisions per turn DefineConstant [ - turns = {10, Label "Number of coil turns"}, + turns = {5, Label "Number of coil turns"}, r = {0.11, Label "Coil radius"}, rc = {0.01, Label "Coil wire radius"}, hc = {0.25, Label "Coil height"}, @@ -142,3 +142,6 @@ Physical Surface(SKIN_TUBE) = CombinedBoundary{ Volume{vol_tube[]}; }; Physical Surface(IN) = in; Physical Surface(OUT) = out; Physical Surface(INF) = {s:s+5}; + +Homology(2) {{COIL,TUBE},{SKIN_COIL, SKIN_TUBE}}; +Cohomology(1) {{AIR},{}}; -- GitLab