Skip to content
Snippets Groups Projects
Commit dc8db1d0 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files 

fix+cleanup elliptic algorithm (didn't work when line mesh orientation was
opposite of original generatrice)
parent 22ebdcc9
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Mesh/2D_Elliptic.cpp
+ 82
− 156
View file @ dc8db1d0
// $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 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
// //
...@@ -27,128 +27,54 @@ ...@@ -27,128 +27,54 @@
extern Mesh *THEM; 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; for(int i = 0; i < 4; i++) {
Curve *GG[444]; if(!compareVertex(&c[i]->beg, &start) &&
Vertex **list; !compareVertex(&c[i]->end, &end)){
double alpha, beta, gamma, u, v, lc, x, y, z; for(int j = 0; j < List_Nbr(c[i]->Vertices); j++)
Vertex *S[4], *v11, *v12, *v13, *v21, *v22, *v23, *v31, *v32, *v33; List_Add(list, List_Pointer(c[i]->Vertices, j));
List_T *l1, *l2, *l3, *l4; return;
if(sur->Method != ELLIPTIC)
return 0;
nb = List_Nbr(sur->Generatrices);
if(nb < 4)
return 0;
if(List_Nbr(sur->TrsfPoints) != 4)
return 0;
for(i = 0; i < nb; 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;
} }
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;
} }
while(compareVertex(&GG[j]->beg, &S[1]));
} }
} }
// step 2 int MeshEllipticSurface(Surface * sur)
l2 = List_Create(20, 10, sizeof(Vertex *)); {
for(i = 0; i < nb; i++) { if(sur->Method != ELLIPTIC)
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; return 0;
}
}
while(compareVertex(&GG[j]->beg, &S[2]));
}
}
// step 3 if(List_Nbr(sur->Generatrices) != 4 ||
l3 = List_Create(20, 10, sizeof(Vertex *)); List_Nbr(sur->TrsfPoints) != 4)
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; return 0;
}
}
while(compareVertex(&GG[j]->beg, &S[3]));
}
}
// step 4 Curve *GG[4];
l4 = List_Create(20, 10, sizeof(Vertex *)); Vertex *S[4];
for(i = 0; i < nb; i++) { for(int i = 0; i < 4; i++) {
if(!compareVertex(&GG[i]->beg, &S[3])) { List_Read(sur->Generatrices, i, &GG[i]);
j = i; List_Read(sur->TrsfPoints, i, &S[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); List_T *l1 = List_Create(20, 10, sizeof(Vertex *));
N2 = List_Nbr(l2); getVertices(S[0], S[1], GG, l1);
N3 = List_Nbr(l3); List_T *l2 = List_Create(20, 10, sizeof(Vertex *));
N4 = List_Nbr(l4); getVertices(S[1], S[2], GG, l2);
if(N1 != N3 || N2 != N4) { 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); Msg(WARNING, "Wrong definition of Elliptic Surface %d", sur->Num);
List_Delete(l1); List_Delete(l1);
List_Delete(l2); List_Delete(l2);
...@@ -159,10 +85,10 @@ int MeshEllipticSurface(Surface * sur) ...@@ -159,10 +85,10 @@ int MeshEllipticSurface(Surface * sur)
sur->Nu = N1; sur->Nu = N1;
sur->Nv = N2; 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(int i = 0; i < N1; i++) {
for(j = 0; j < N2; j++) { for(int j = 0; j < N2; j++) {
if(i == 0) { if(i == 0) {
List_Read(l4, N2 - j - 1, &list[i + N1 * j]); List_Read(l4, N2 - j - 1, &list[i + N1 * j]);
} }
...@@ -176,21 +102,21 @@ int MeshEllipticSurface(Surface * sur) ...@@ -176,21 +102,21 @@ int MeshEllipticSurface(Surface * sur)
List_Read(l3, N1 - i - 1, &list[i + N1 * j]); List_Read(l3, N1 - i - 1, &list[i + N1 * j]);
} }
else { else {
u = 1. - 2. * (double)i / double (N1 - 1); double u = 1. - 2. * (double)i / double (N1 - 1);
v = 1. - 2. * (double)j / double (N2 - 1); double v = 1. - 2. * (double)j / double (N2 - 1);
x = 0.25 * ((S[0]->Pos.X * (1 + u) * (1. + v)) + double x = 0.25 * ((S[0]->Pos.X * (1 + u) * (1. + v)) +
(S[1]->Pos.X * (1 - u) * (1. + v)) + (S[1]->Pos.X * (1 - u) * (1. + v)) +
(S[2]->Pos.X * (1 - u) * (1. - v)) + (S[2]->Pos.X * (1 - u) * (1. - v)) +
(S[3]->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)) + double y = 0.25 * ((S[0]->Pos.Y * (1 + u) * (1. + v)) +
(S[1]->Pos.Y * (1 - u) * (1. + v)) + (S[1]->Pos.Y * (1 - u) * (1. + v)) +
(S[2]->Pos.Y * (1 - u) * (1. - v)) + (S[2]->Pos.Y * (1 - u) * (1. - v)) +
(S[3]->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)) + double z = 0.25 * ((S[0]->Pos.Z * (1 + u) * (1. + v)) +
(S[1]->Pos.Z * (1 - u) * (1. + v)) + (S[1]->Pos.Z * (1 - u) * (1. + v)) +
(S[2]->Pos.Z * (1 - u) * (1. - v)) + (S[2]->Pos.Z * (1 - u) * (1. - v)) +
(S[3]->Pos.Z * (1 + u) * (1. - v))); (S[3]->Pos.Z * (1 + u) * (1. - v)));
lc = 0.25 * ((S[0]->lc * (1 + u) * (1. + v)) + double lc = 0.25 * ((S[0]->lc * (1 + u) * (1. + v)) +
(S[1]->lc * (1 - u) * (1. + v)) + (S[1]->lc * (1 - u) * (1. + v)) +
(S[2]->lc * (1 - u) * (1. - v)) + (S[2]->lc * (1 - u) * (1. - v)) +
(S[3]->lc * (1 + u) * (1. - v))); (S[3]->lc * (1 + u) * (1. - v)));
...@@ -206,28 +132,28 @@ int MeshEllipticSurface(Surface * sur) ...@@ -206,28 +132,28 @@ int MeshEllipticSurface(Surface * sur)
List_Delete(l3); List_Delete(l3);
List_Delete(l4); List_Delete(l4);
k = 0; int k = 0;
while(1) { while(1) {
k++; k++;
if(k > 1000) if(k > 1000)
break; break;
for(i = 1; i < N1 - 1; i++) { for(int i = 1; i < N1 - 1; i++) {
for(j = 1; j < N2 - 1; j++) { for(int j = 1; j < N2 - 1; j++) {
v11 = list[i - 1 + N1 * (j - 1)]; Vertex *v11 = list[i - 1 + N1 * (j - 1)];
v12 = list[i + N1 * (j - 1)]; Vertex *v12 = list[i + N1 * (j - 1)];
v13 = list[i + 1 + N1 * (j - 1)]; Vertex *v13 = list[i + 1 + N1 * (j - 1)];
v21 = list[i - 1 + N1 * (j)]; Vertex *v21 = list[i - 1 + N1 * (j)];
v22 = list[i + N1 * (j)]; Vertex *v22 = list[i + N1 * (j)];
v23 = list[i + 1 + N1 * (j)]; Vertex *v23 = list[i + 1 + N1 * (j)];
v31 = list[i - 1 + N1 * (j + 1)]; Vertex *v31 = list[i - 1 + N1 * (j + 1)];
v32 = list[i + N1 * (j + 1)]; Vertex *v32 = list[i + N1 * (j + 1)];
v33 = list[i + 1 + N1 * (j + 1)]; Vertex *v33 = list[i + 1 + N1 * (j + 1)];
alpha = 0.25 * (DSQR(v23->Pos.X - v21->Pos.X) + double alpha = 0.25 * (DSQR(v23->Pos.X - v21->Pos.X) +
DSQR(v23->Pos.Y - v21->Pos.Y)); DSQR(v23->Pos.Y - v21->Pos.Y));
gamma = 0.25 * (DSQR(v32->Pos.X - v12->Pos.X) + double gamma = 0.25 * (DSQR(v32->Pos.X - v12->Pos.X) +
DSQR(v32->Pos.Y - v12->Pos.Y)); DSQR(v32->Pos.Y - v12->Pos.Y));
beta = 0.0625 * ((v32->Pos.X - v12->Pos.X) * double beta = 0.0625 * ((v32->Pos.X - v12->Pos.X) *
(v23->Pos.X - v21->Pos.X) + (v23->Pos.X - v21->Pos.X) +
(v32->Pos.Y - v12->Pos.Y) * (v32->Pos.Y - v12->Pos.Y) *
(v23->Pos.Y - v21->Pos.Y)); (v23->Pos.Y - v21->Pos.Y));
...@@ -250,8 +176,8 @@ int MeshEllipticSurface(Surface * sur) ...@@ -250,8 +176,8 @@ int MeshEllipticSurface(Surface * sur)
} }
} }
} }
for(i = 0; i < N1 - 1; i++) { for(int i = 0; i < N1 - 1; i++) {
for(j = 0; j < N2 - 1; j++) { for(int j = 0; j < N2 - 1; j++) {
if(sur->Recombine) { if(sur->Recombine) {
Quadrangle *q = Create_Quadrangle Quadrangle *q = Create_Quadrangle
(list[(i) + N1 * (j)], list[(i + 1) + N1 * (j)], (list[(i) + N1 * (j)], list[(i + 1) + N1 * (j)],
......
benchmarks/bugs/bug_elliptic.geo deleted 100644 → 0
+ 0
− 45
View file @ 22ebdcc9
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};
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment