From d62617cb2cacf0265ecc950c62005eb2dfcb3ea9 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Thu, 24 Jan 2002 17:48:28 +0000
Subject: [PATCH] Fix precision problem with asin()
---
Common/Numeric.cpp | 8 +++++++-
Common/Numeric.h | 1 +
Mesh/Create.cpp | 15 +++++++++------
3 files changed, 17 insertions(+), 7 deletions(-)
diff --git a/Common/Numeric.cpp b/Common/Numeric.cpp
index 43a6ac8d5a..176c4162db 100644
--- a/Common/Numeric.cpp
+++ b/Common/Numeric.cpp
@@ -1,4 +1,4 @@
-// $Id: Numeric.cpp,v 1.3 2001-12-03 08:41:43 geuzaine Exp $
+// $Id: Numeric.cpp,v 1.4 2002-01-24 17:48:28 geuzaine Exp $
#include "Gmsh.h"
#include "Numeric.h"
@@ -20,6 +20,12 @@ double myacos (double a){
return acos (a);
}
+double myasin(double a){
+ if(a<=-1.) return -Pi/2.;
+ else if(a>=1.) return Pi/2.;
+ else return asin(a);
+}
+
void prodve (double a[3], double b[3], double c[3]){
c[2] = a[0] * b[1] - a[1] * b[0];
c[1] = -a[0] * b[2] + a[2] * b[0];
diff --git a/Common/Numeric.h b/Common/Numeric.h
index 7d00a54e85..98db2432c1 100644
--- a/Common/Numeric.h
+++ b/Common/Numeric.h
@@ -36,6 +36,7 @@
double myatan2 (double a, double b);
double myacos (double a);
+double myasin (double a);
void prodve (double a[3], double b[3], double c[3]);
void prosca (double a[3], double b[3], double *c);
void norme (double a[3]);
diff --git a/Mesh/Create.cpp b/Mesh/Create.cpp
index ed9a1ab9c4..87f5c7cbec 100644
--- a/Mesh/Create.cpp
+++ b/Mesh/Create.cpp
@@ -1,4 +1,4 @@
-// $Id: Create.cpp,v 1.33 2002-01-23 16:28:00 geuzaine Exp $
+// $Id: Create.cpp,v 1.34 2002-01-24 17:48:28 geuzaine Exp $
#include "Gmsh.h"
#include "Numeric.h"
@@ -318,17 +318,20 @@ void End_Curve (Curve * c){
A1 = A3 = 0.;
f1 = f2 = R ;
}
- else{
+ else{
f1 = sqrt(1./sol[0]);
f2 = sqrt(1./sol[1]);
+ // myasin() permet de contourner les problemes de precision
+ // sur y1/f2 ou y3/f2, qui peuvent legerement etre hors de
+ // [-1,1]
if(x1 < 0)
- A1 = -asin(y1/f2) + A4 + Pi;
+ A1 = -myasin(y1/f2) + A4 + Pi;
else
- A1 = asin(y1/f2) + A4;
+ A1 = myasin(y1/f2) + A4;
if(x3 < 0)
- A3 = -asin(y3/f2) + A4 + Pi;
+ A3 = -myasin(y3/f2) + A4 + Pi;
else
- A3 = asin(y3/f2) + A4;
+ A3 = myasin(y3/f2) + A4;
}
}
else{
--
GitLab