diff --git a/Mesh/2D_Elliptic.cpp b/Mesh/2D_Elliptic.cpp
index 7c8be057f1e374faef2d40abe3b869b535b7885b..f3225c0fa4f360f8b83a7e84ac57da159920d221 100644
--- a/Mesh/2D_Elliptic.cpp
+++ b/Mesh/2D_Elliptic.cpp
@@ -1,4 +1,4 @@
-// $Id: 2D_Elliptic.cpp,v 1.19 2004-05-27 20:49:03 geuzaine Exp $
+// $Id: 2D_Elliptic.cpp,v 1.20 2004-12-17 04:12:13 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -27,128 +27,54 @@
extern Mesh *THEM;
-int MeshEllipticSurface(Surface * sur)
+static void getVertices(Vertex *start, Vertex *end, Curve *c[4], List_T *list)
{
- int i, j, k, nb, N1, N2, N3, N4;
- Curve *GG[444];
- Vertex **list;
- double alpha, beta, gamma, u, v, lc, x, y, z;
- Vertex *S[4], *v11, *v12, *v13, *v21, *v22, *v23, *v31, *v32, *v33;
- List_T *l1, *l2, *l3, *l4;
+ for(int i = 0; i < 4; i++) {
+ if(!compareVertex(&c[i]->beg, &start) &&
+ !compareVertex(&c[i]->end, &end)){
+ for(int j = 0; j < List_Nbr(c[i]->Vertices); j++)
+ List_Add(list, List_Pointer(c[i]->Vertices, j));
+ return;
+ }
+ if(!compareVertex(&c[i]->beg, &end) &&
+ !compareVertex(&c[i]->end, &start)){
+ for(int j = 0; j < List_Nbr(c[i]->Vertices); j++)
+ List_Add(list, List_Pointer(c[i]->Vertices, List_Nbr(c[i]->Vertices) - 1 - j));
+ return;
+ }
+ }
+}
+int MeshEllipticSurface(Surface * sur)
+{
if(sur->Method != ELLIPTIC)
return 0;
- nb = List_Nbr(sur->Generatrices);
-
- if(nb < 4)
- return 0;
-
- if(List_Nbr(sur->TrsfPoints) != 4)
+ if(List_Nbr(sur->Generatrices) != 4 ||
+ List_Nbr(sur->TrsfPoints) != 4)
return 0;
- for(i = 0; i < nb; i++) {
+ Curve *GG[4];
+ Vertex *S[4];
+ for(int i = 0; i < 4; i++) {
List_Read(sur->Generatrices, i, &GG[i]);
List_Read(sur->TrsfPoints, i, &S[i]);
}
- N1 = N2 = N3 = N4 = 0;
-
- // step 1
- l1 = List_Create(20, 10, sizeof(Vertex *));
- for(i = 0; i < nb; i++) {
- if(!compareVertex(&GG[i]->beg, &S[0])) {
- j = i;
- do {
- if(!List_Nbr(l1))
- List_Add(l1, List_Pointer(GG[j]->Vertices, 0));
- for(k = 1; k < List_Nbr(GG[j]->Vertices); k++)
- List_Add(l1, List_Pointer(GG[j]->Vertices, k));
- j = (j + 1 < nb) ? j + 1 : 0;
- if(j == i){
- Msg(WARNING, "Wrong definition of Elliptic Surface %d", sur->Num);
- List_Delete(l1);
- return 0;
- }
- }
- while(compareVertex(&GG[j]->beg, &S[1]));
- }
- }
-
- // step 2
- l2 = List_Create(20, 10, sizeof(Vertex *));
- for(i = 0; i < nb; i++) {
- if(!compareVertex(&GG[i]->beg, &S[1])) {
- j = i;
- do {
- if(!List_Nbr(l2))
- List_Add(l2, List_Pointer(GG[j]->Vertices, 0));
- for(k = 1; k < List_Nbr(GG[j]->Vertices); k++)
- List_Add(l2, List_Pointer(GG[j]->Vertices, k));
- j = (j + 1 < nb) ? j + 1 : 0;
- if(j == i){
- Msg(WARNING, "Wrong definition of Elliptic Surface %d", sur->Num);
- List_Delete(l1);
- List_Delete(l2);
- return 0;
- }
- }
- while(compareVertex(&GG[j]->beg, &S[2]));
- }
- }
-
- // step 3
- l3 = List_Create(20, 10, sizeof(Vertex *));
- for(i = 0; i < nb; i++) {
- if(!compareVertex(&GG[i]->beg, &S[2])) {
- j = i;
- do {
- if(!List_Nbr(l3))
- List_Add(l3, List_Pointer(GG[j]->Vertices, 0));
- for(k = 1; k < List_Nbr(GG[j]->Vertices); k++)
- List_Add(l3, List_Pointer(GG[j]->Vertices, k));
- j = (j + 1 < nb) ? j + 1 : 0;
- if(j == i){
- Msg(WARNING, "Wrong definition of Elliptic Surface %d", sur->Num);
- List_Delete(l1);
- List_Delete(l2);
- List_Delete(l3);
- return 0;
- }
- }
- while(compareVertex(&GG[j]->beg, &S[3]));
- }
- }
-
- // step 4
- l4 = List_Create(20, 10, sizeof(Vertex *));
- for(i = 0; i < nb; i++) {
- if(!compareVertex(&GG[i]->beg, &S[3])) {
- j = i;
- do {
- if(!List_Nbr(l4))
- List_Add(l4, List_Pointer(GG[j]->Vertices, 0));
- for(k = 1; k < List_Nbr(GG[j]->Vertices); k++)
- List_Add(l4, List_Pointer(GG[j]->Vertices, k));
- j = (j + 1 < nb) ? j + 1 : 0;
- if(j == i){
- Msg(WARNING, "Wrong definition of Elliptic Surface %d", sur->Num);
- List_Delete(l1);
- List_Delete(l2);
- List_Delete(l3);
- List_Delete(l4);
- return 0;
- }
- }
- while(compareVertex(&GG[j]->beg, &S[0]));
- }
- }
-
- N1 = List_Nbr(l1);
- N2 = List_Nbr(l2);
- N3 = List_Nbr(l3);
- N4 = List_Nbr(l4);
- if(N1 != N3 || N2 != N4) {
+ List_T *l1 = List_Create(20, 10, sizeof(Vertex *));
+ getVertices(S[0], S[1], GG, l1);
+ List_T *l2 = List_Create(20, 10, sizeof(Vertex *));
+ getVertices(S[1], S[2], GG, l2);
+ List_T *l3 = List_Create(20, 10, sizeof(Vertex *));
+ getVertices(S[2], S[3], GG, l3);
+ List_T *l4 = List_Create(20, 10, sizeof(Vertex *));
+ getVertices(S[3], S[0], GG, l4);
+
+ int N1 = List_Nbr(l1);
+ int N2 = List_Nbr(l2);
+ int N3 = List_Nbr(l3);
+ int N4 = List_Nbr(l4);
+ if(!N1 || !N2 || !N3 || !N4 || (N1 != N3) || (N2 != N4)) {
Msg(WARNING, "Wrong definition of Elliptic Surface %d", sur->Num);
List_Delete(l1);
List_Delete(l2);
@@ -159,10 +85,10 @@ int MeshEllipticSurface(Surface * sur)
sur->Nu = N1;
sur->Nv = N2;
- list = (Vertex **) Malloc(N1 * N2 * sizeof(Vertex *));
+ Vertex **list = (Vertex **) Malloc(N1 * N2 * sizeof(Vertex *));
- for(i = 0; i < N1; i++) {
- for(j = 0; j < N2; j++) {
+ for(int i = 0; i < N1; i++) {
+ for(int j = 0; j < N2; j++) {
if(i == 0) {
List_Read(l4, N2 - j - 1, &list[i + N1 * j]);
}
@@ -176,24 +102,24 @@ int MeshEllipticSurface(Surface * sur)
List_Read(l3, N1 - i - 1, &list[i + N1 * j]);
}
else {
- u = 1. - 2. * (double)i / double (N1 - 1);
- v = 1. - 2. * (double)j / double (N2 - 1);
- x = 0.25 * ((S[0]->Pos.X * (1 + u) * (1. + v)) +
- (S[1]->Pos.X * (1 - u) * (1. + v)) +
- (S[2]->Pos.X * (1 - u) * (1. - v)) +
- (S[3]->Pos.X * (1 + u) * (1. - v)));
- y = 0.25 * ((S[0]->Pos.Y * (1 + u) * (1. + v)) +
- (S[1]->Pos.Y * (1 - u) * (1. + v)) +
- (S[2]->Pos.Y * (1 - u) * (1. - v)) +
- (S[3]->Pos.Y * (1 + u) * (1. - v)));
- z = 0.25 * ((S[0]->Pos.Z * (1 + u) * (1. + v)) +
- (S[1]->Pos.Z * (1 - u) * (1. + v)) +
- (S[2]->Pos.Z * (1 - u) * (1. - v)) +
- (S[3]->Pos.Z * (1 + u) * (1. - v)));
- lc = 0.25 * ((S[0]->lc * (1 + u) * (1. + v)) +
- (S[1]->lc * (1 - u) * (1. + v)) +
- (S[2]->lc * (1 - u) * (1. - v)) +
- (S[3]->lc * (1 + u) * (1. - v)));
+ double u = 1. - 2. * (double)i / double (N1 - 1);
+ double v = 1. - 2. * (double)j / double (N2 - 1);
+ double x = 0.25 * ((S[0]->Pos.X * (1 + u) * (1. + v)) +
+ (S[1]->Pos.X * (1 - u) * (1. + v)) +
+ (S[2]->Pos.X * (1 - u) * (1. - v)) +
+ (S[3]->Pos.X * (1 + u) * (1. - v)));
+ double y = 0.25 * ((S[0]->Pos.Y * (1 + u) * (1. + v)) +
+ (S[1]->Pos.Y * (1 - u) * (1. + v)) +
+ (S[2]->Pos.Y * (1 - u) * (1. - v)) +
+ (S[3]->Pos.Y * (1 + u) * (1. - v)));
+ double z = 0.25 * ((S[0]->Pos.Z * (1 + u) * (1. + v)) +
+ (S[1]->Pos.Z * (1 - u) * (1. + v)) +
+ (S[2]->Pos.Z * (1 - u) * (1. - v)) +
+ (S[3]->Pos.Z * (1 + u) * (1. - v)));
+ double lc = 0.25 * ((S[0]->lc * (1 + u) * (1. + v)) +
+ (S[1]->lc * (1 - u) * (1. + v)) +
+ (S[2]->lc * (1 - u) * (1. - v)) +
+ (S[3]->lc * (1 + u) * (1. - v)));
list[i + N1 * j] = Create_Vertex(++THEM->MaxPointNum, x, y, z, lc, 0.0);
Tree_Insert(sur->Vertices, &list[i + N1 * j]);
List_Add(sur->TrsfVertices, &list[i + N1 * j]);
@@ -206,31 +132,31 @@ int MeshEllipticSurface(Surface * sur)
List_Delete(l3);
List_Delete(l4);
- k = 0;
+ int k = 0;
while(1) {
k++;
if(k > 1000)
break;
- for(i = 1; i < N1 - 1; i++) {
- for(j = 1; j < N2 - 1; j++) {
- v11 = list[i - 1 + N1 * (j - 1)];
- v12 = list[i + N1 * (j - 1)];
- v13 = list[i + 1 + N1 * (j - 1)];
- v21 = list[i - 1 + N1 * (j)];
- v22 = list[i + N1 * (j)];
- v23 = list[i + 1 + N1 * (j)];
- v31 = list[i - 1 + N1 * (j + 1)];
- v32 = list[i + N1 * (j + 1)];
- v33 = list[i + 1 + N1 * (j + 1)];
-
- alpha = 0.25 * (DSQR(v23->Pos.X - v21->Pos.X) +
- DSQR(v23->Pos.Y - v21->Pos.Y));
- gamma = 0.25 * (DSQR(v32->Pos.X - v12->Pos.X) +
- DSQR(v32->Pos.Y - v12->Pos.Y));
- beta = 0.0625 * ((v32->Pos.X - v12->Pos.X) *
- (v23->Pos.X - v21->Pos.X) +
- (v32->Pos.Y - v12->Pos.Y) *
- (v23->Pos.Y - v21->Pos.Y));
+ for(int i = 1; i < N1 - 1; i++) {
+ for(int j = 1; j < N2 - 1; j++) {
+ Vertex *v11 = list[i - 1 + N1 * (j - 1)];
+ Vertex *v12 = list[i + N1 * (j - 1)];
+ Vertex *v13 = list[i + 1 + N1 * (j - 1)];
+ Vertex *v21 = list[i - 1 + N1 * (j)];
+ Vertex *v22 = list[i + N1 * (j)];
+ Vertex *v23 = list[i + 1 + N1 * (j)];
+ Vertex *v31 = list[i - 1 + N1 * (j + 1)];
+ Vertex *v32 = list[i + N1 * (j + 1)];
+ Vertex *v33 = list[i + 1 + N1 * (j + 1)];
+
+ double alpha = 0.25 * (DSQR(v23->Pos.X - v21->Pos.X) +
+ DSQR(v23->Pos.Y - v21->Pos.Y));
+ double gamma = 0.25 * (DSQR(v32->Pos.X - v12->Pos.X) +
+ DSQR(v32->Pos.Y - v12->Pos.Y));
+ double beta = 0.0625 * ((v32->Pos.X - v12->Pos.X) *
+ (v23->Pos.X - v21->Pos.X) +
+ (v32->Pos.Y - v12->Pos.Y) *
+ (v23->Pos.Y - v21->Pos.Y));
v22->Pos.X = 0.5 * (alpha * (v32->Pos.X + v12->Pos.X) +
gamma * (v23->Pos.X + v21->Pos.X) -
@@ -250,8 +176,8 @@ int MeshEllipticSurface(Surface * sur)
}
}
}
- for(i = 0; i < N1 - 1; i++) {
- for(j = 0; j < N2 - 1; j++) {
+ for(int i = 0; i < N1 - 1; i++) {
+ for(int j = 0; j < N2 - 1; j++) {
if(sur->Recombine) {
Quadrangle *q = Create_Quadrangle
(list[(i) + N1 * (j)], list[(i + 1) + N1 * (j)],
diff --git a/benchmarks/bugs/bug_elliptic.geo b/benchmarks/bugs/bug_elliptic.geo
deleted file mode 100644
index 3710f5cc07d908a4a4780dad63bc34b58051ebba..0000000000000000000000000000000000000000
--- a/benchmarks/bugs/bug_elliptic.geo
+++ /dev/null
@@ -1,45 +0,0 @@
-size = 0.005;
-Point(1) = {0, 0, 0, size};
-Point(2) = {0.0025, 0, 0, size};
-Point(3) = {0.016, 0, 0, size};
-Point(4) = {0.016, 0.015, 0, size};
-Point(5) = {0, 0.015, 0, size};
-Point(6) = {0, 0.0025, 0, size};
-Point(7) = {0.001767766953, 0.001767766953, 0, size};
-Point(8) = {0.016, 0.025, 0, size};
-Point(9) = {0, 0.025, 0, size};
-Line(1) = {2,3};
-Line(2) = {3,4};
-Line(3) = {4,5};
-Line(4) = {5,6};
-Line(5) = {7,4};
-Line(6) = {4,8};
-Line(7) = {8,9};
-Line(8) = {9,5};
-Circle(9) = {2,1,7};
-Circle(10) = {7,1,6};
-Line Loop(11) = {-2,-1,9,5};
-Plane Surface(12) = {11};
-Line Loop(13) = {4,-10,5,3};
-Plane Surface(14) = {13};
-Line Loop(15) = {-3,6,7,8};
-Plane Surface(16) = {15};
-Transfinite Line {1} = 15 Using Progression 1.1;
-Transfinite Line {4} = 15 Using Progression 0.9;
-Transfinite Line {5} = 15 Using Progression 1.15;
-Transfinite Line {9,2} = 11;
-Transfinite Line {10,3,7} = 10;
-Transfinite Line {6,8} = 6;
-
-//Elliptic Surface {12} = {7,4,3,2};
-Elliptic Surface {12} = {2,3,4,7};// aieieieie
-Elliptic Surface {14} = {6,7,4,5};
-Elliptic Surface {16} = {5,4,8,9};
-
-Recombine Surface {12,14,16};
-Physical Surface(1) = {-12,14,16};
-Physical Line(11) = {1};
-Physical Line(12) = {2,6};
-Physical Line(13) = {7};
-Physical Line(14) = {4,8};
-Physical Line(15) = {9,10};