From a02a07ab883df64f4d7454e7568a62402d9c2df0 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 16 Nov 2001 19:35:26 +0000
Subject: [PATCH] *** empty log message ***

---
 Mesh/3D_Mesh.cpp         |  6 ++--
 Mesh/Create.cpp          | 73 ++++++++++++++++++++++------------------
 Mesh/Interpolation.cpp   |  9 +++--
 benchmarks/3d/Sphere.geo |  2 +-
 benchmarks/3d/Torus.geo  |  2 ++
 benchmarks/3d/p19.geo    |  2 +-
 6 files changed, 54 insertions(+), 40 deletions(-)

diff --git a/Mesh/3D_Mesh.cpp b/Mesh/3D_Mesh.cpp
index d9ee2ccf48..a41360c6dd 100644
--- a/Mesh/3D_Mesh.cpp
+++ b/Mesh/3D_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: 3D_Mesh.cpp,v 1.32 2001-10-31 16:33:46 remacle Exp $
+// $Id: 3D_Mesh.cpp,v 1.33 2001-11-16 19:35:26 remacle Exp $
 
 /*
  
@@ -614,7 +614,7 @@ void Convex_Hull_Mesh (List_T * Points, Mesh * m){
   SimplexInteriorCheck ick;
   SimplexInBox         ibx (LC3D/30.);
 
-  fastSearch = new AOMD::mGeomSearch < Simplex, SimplexInBox, SimplexInteriorCheck >
+    fastSearch = new AOMD::mGeomSearch < Simplex, SimplexInBox, SimplexInteriorCheck >
     (  m->Grid.min.X,m->Grid.max.X,
        m->Grid.min.Y,m->Grid.max.Y,
        m->Grid.min.Z,m->Grid.max.Z,
@@ -640,7 +640,7 @@ void Convex_Hull_Mesh (List_T * Points, Mesh * m){
 		break;
 	      }
 	  }
-      }
+      }    
     else
       {
 	nbFast++;
diff --git a/Mesh/Create.cpp b/Mesh/Create.cpp
index 4b294ae9f6..fa15641a2e 100644
--- a/Mesh/Create.cpp
+++ b/Mesh/Create.cpp
@@ -1,4 +1,4 @@
-// $Id: Create.cpp,v 1.26 2001-10-29 08:52:20 geuzaine Exp $
+// $Id: Create.cpp,v 1.27 2001-11-16 19:35:26 remacle Exp $
 
 #include "Gmsh.h"
 #include "Numeric.h"
@@ -273,8 +273,7 @@ void End_Curve (Curve * c){
     mat[1][2] = Curve->Circle.invmat[2][1] = m[2];
     mat[0][0] = Curve->Circle.invmat[0][0] = dir12[0];
     mat[0][1] = Curve->Circle.invmat[1][0] = dir12[1];
-    mat[0][2] = Curve->Circle.invmat[2][0] = dir12[2];
-    
+    mat[0][2] = Curve->Circle.invmat[2][0] = dir12[2];    
     if(CTX.geom.old_circle){
       if(n[0] == 0.0 && n[1] == 0.0){
         mat[2][0] = Curve->Circle.invmat[0][2] = 0;
@@ -311,37 +310,44 @@ void End_Curve (Curve * c){
       A4 = angle_02pi (A4);
     if (A1 >= A3)
       A3 += DP;
-    if (A4 > A1)
-      A4 -= DP;
+    //    if (A4 > A1)
+    //      A4 -= DP;
     
     if (v[4]){
-      AX = (A1 - A4);
-      sys[0][0] = cos (AX) * cos (A4);
-      sys[0][1] = -sin (AX) * sin (A4);
-      sys[1][0] = cos (AX) * sin (A4);
-      sys[1][1] = sin (AX) * cos (A4);
-      rhs[0] = v1.Pos.X;
-      rhs[1] = v1.Pos.Y;
-      det = sys[0][0] * sys[1][1] - sys[1][0] * sys[0][1];
-      if (det < 1.e-12){
-        AX = (A3 - A4);
-        sys[0][0] = cos (AX) * cos (A4);
-        sys[0][1] = -sin (AX) * sin (A4);
-        sys[1][0] = cos (AX) * sin (A4);
-        sys[1][1] = sin (AX) * cos (A4);
-        rhs[0] = v3.Pos.X;
-        rhs[1] = v3.Pos.Y;
-        det = sys[0][0] * sys[1][1] - sys[1][0] * sys[0][1];
-      }
-      if (det < 1.e-12){
-        f1 = DMAX (R, R2);
-        f2 = DMIN (R, R2);
-      }
-      else{
-        sys2x2 (sys, rhs, sol);
-        f1 = sol[0];
-        f2 = sol[1];
-      }
+      AX = (A4);
+      double x1 = v1.Pos.X * cos (AX) + v1.Pos.Y * sin(AX);
+      double y1 = -v1.Pos.X * sin (AX) + v1.Pos.Y * cos(AX); 
+      double x3 = v3.Pos.X * cos (AX) + v3.Pos.Y * sin(AX);
+      double y3 = -v3.Pos.X * sin (AX) + v3.Pos.Y * cos(AX); 
+      sys[0][0] = x1 * x1;
+      sys[0][1] = y1 * y1;
+      sys[1][0] = x3 * x3;
+      sys[1][1] = y3 * y3;
+
+      rhs[0] = 1;
+      rhs[1] = 1;
+
+      //      printf("AX = %lf\n",AX*180./M_PI);
+      //      printf("%lf %lf %lf %lf\n", v1.Pos.X , v1.Pos.Y,x1,y1);
+      //      printf("%lf %lf %lf %lf\n", v3.Pos.X , v3.Pos.Y,x3,y3);
+
+
+      sys2x2 (sys, rhs, sol);
+      //      printf("%lf %lf   %lf = %lf \n",sys[0][0],sys[0][1],rhs[0],sol[0]);
+      //      printf("%lf %lf   %lf = %lf\n",sys[1][0],sys[1][1],rhs[1],sol[1]);
+      if(sol[0] < 0)Msg(FATAL, "Ellipsis %d invalid", Curve->Num);	
+      if(sol[1] < 0)Msg(FATAL, "Ellipsis %d invalid", Curve->Num);	
+      f1 = sqrt(1./sol[0]);
+      f2 = sqrt(1./sol[1]);
+      if(x1 > 0)
+	A1 = asin(y1/f2) + A4; 
+      else 
+	A1 = -asin(y1/f2) + A4; 
+
+      if(x3 > 0)
+	A3 = asin(y3/f2) + A4; 
+      else 
+	A3 = -asin(y3/f2) + A4; 
     }
     else{
       f1 = f2 = R;
@@ -356,6 +362,9 @@ void End_Curve (Curve * c){
     for (i = 0; i < 4; i++)
       Curve->Circle.v[i] = v[i];
 
+    //    double xxx = 180./M_PI;
+    //    printf("%d %lf %lf %lf %lf %lf\n",Curve->Num,f1,f2,A1*xxx,A3*xxx,A4*xxx);
+
     /*
     if (!c->Circle.done){
       float proj[4][4];
diff --git a/Mesh/Interpolation.cpp b/Mesh/Interpolation.cpp
index 7281f88d4d..a2e9c5c114 100644
--- a/Mesh/Interpolation.cpp
+++ b/Mesh/Interpolation.cpp
@@ -1,4 +1,4 @@
-// $Id: Interpolation.cpp,v 1.14 2001-08-17 10:40:23 geuzaine Exp $
+// $Id: Interpolation.cpp,v 1.15 2001-11-16 19:35:26 remacle Exp $
 
 #include "Gmsh.h"
 #include "Numeric.h"
@@ -77,9 +77,12 @@ Vertex InterpolateCurve (Curve * Curve, double u, int derivee){
     /* pour les ellipses */
     teta -= Curve->Circle.incl;
     
-    V.Pos.X = Curve->Circle.f1 * cos (teta) * cos (Curve->Circle.incl) -
+    V.Pos.X = 
+      Curve->Circle.f1 * cos (teta) * cos (Curve->Circle.incl) -
       Curve->Circle.f2 * sin (teta) * sin (Curve->Circle.incl);
-    V.Pos.Y = Curve->Circle.f1 * cos (teta) * sin (Curve->Circle.incl) +
+
+    V.Pos.Y = 
+      Curve->Circle.f1 * cos (teta) * sin (Curve->Circle.incl) +
       Curve->Circle.f2 * sin (teta) * cos (Curve->Circle.incl);
     V.Pos.Z = 0.0;
     Projette (&V, Curve->Circle.invmat);
diff --git a/benchmarks/3d/Sphere.geo b/benchmarks/3d/Sphere.geo
index a507718634..cdd6c6cf5c 100644
--- a/benchmarks/3d/Sphere.geo
+++ b/benchmarks/3d/Sphere.geo
@@ -1,4 +1,4 @@
-lc = .1;
+lc = .03;
 Point(1) = {0.0,0.0,0.0,lc};
 Point(2) = {1,0.0,0.0,lc};
 Point(3) = {0,1,0.0,lc};
diff --git a/benchmarks/3d/Torus.geo b/benchmarks/3d/Torus.geo
index 5848a1f2f2..771b08b3ae 100644
--- a/benchmarks/3d/Torus.geo
+++ b/benchmarks/3d/Torus.geo
@@ -15,3 +15,5 @@ Plane Surface(6) = {5};
 Extrude Surface{6, {0.0,1,0}, {0,0.0,0.0}, 1*3.14159/2};         
        
 Coherence;         
+Surface Loop(29) = {6,15,19,23,27,28};
+Volume(30) = {29};
diff --git a/benchmarks/3d/p19.geo b/benchmarks/3d/p19.geo
index 344119dd0d..4e363a9ca1 100644
--- a/benchmarks/3d/p19.geo
+++ b/benchmarks/3d/p19.geo
@@ -1,4 +1,4 @@
-fact     = 0.4 ;
+fact     = 0.1 ;
 rondelle = fact * 0.01;
 iris     = fact * 0.004;
 size     = fact * 0.01;
-- 
GitLab