From df90a675a2afa014a3879cc411174c52441a2fd9 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Wed, 23 Jun 2004 18:52:45 +0000
Subject: [PATCH] Replaced my orientation hack with a rigorous algorithm...
Seems to work nicely (and is general).
---
Graphics/Geom.cpp | 76 +++++++++++++++++++++++++++++++++++----------
Mesh/2D_Mesh.cpp | 67 ++++++++++++++++-----------------------
Mesh/3D_Extrude.cpp | 3 +-
Mesh/Mesh.h | 2 --
doc/VERSIONS | 6 ++--
5 files changed, 91 insertions(+), 63 deletions(-)
diff --git a/Graphics/Geom.cpp b/Graphics/Geom.cpp
index 30c5cffdc8..b439ae5f8c 100644
--- a/Graphics/Geom.cpp
+++ b/Graphics/Geom.cpp
@@ -1,4 +1,4 @@
-// $Id: Geom.cpp,v 1.67 2004-06-23 03:57:43 geuzaine Exp $
+// $Id: Geom.cpp,v 1.68 2004-06-23 18:52:45 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -207,7 +207,7 @@ void Draw_Curve(void *a, void *b)
// Surfaces
-void put_Z(Vertex * v, Surface * s)
+void putZ(Vertex * v, Surface * s)
{
Vertex V;
V.Pos.X = s->a;
@@ -255,6 +255,45 @@ int isPointOnPlanarSurface(Surface * S, double X, double Y, double Z,
return 0;
}
+void getPlaneSurfaceNormal(Surface *s, double x, double y, double z, double n[3])
+{
+ double t1[3], t2[3];
+ Curve *c;
+
+ for(int i = 0; i < 3; i++)
+ n[i] = 0.0;
+
+ if(s->Typ != MSH_SURF_PLAN)
+ return;
+
+ // We cannot use s->plan or InterpolateSurface here (they both rely
+ // on the mean plane, which borks the orientation). So we need to
+ // use this trick:
+
+ if(List_Nbr(s->Generatrices) > 1){
+ List_Read(s->Generatrices, 0, &c);
+ t1[0] = x - c->beg->Pos.X;
+ t1[1] = y - c->beg->Pos.Y;
+ t1[2] = z - c->beg->Pos.Z;
+ for(int i = 1; i < List_Nbr(s->Generatrices); i++){
+ List_Read(s->Generatrices, i, &c);
+ t2[0] = x - c->beg->Pos.X;
+ t2[1] = y - c->beg->Pos.Y;
+ t2[2] = z - c->beg->Pos.Z;
+ prodve(t1, t2, n);
+ if(norme(n))
+ break;
+ }
+ }
+ if(List_Nbr(s->Generatrices) < 2 || !norme(n)){
+ Msg(WARNING, "Reverting to mean plane normal for surface %d", s->Num);
+ n[0] = s->a;
+ n[1] = s->b;
+ n[2] = s->c;
+ norme(n);
+ }
+}
+
void Draw_Triangulated_Surface(Surface * s)
{
int k = 0;
@@ -344,7 +383,7 @@ void Draw_Plane_Surface(Surface * s)
for(i = 0; i < 4; i++) {
V[i].Pos.Z = 0.0;
- put_Z(&V[i], s);
+ putZ(&V[i], s);
}
n[0] = s->plan[2][0];
@@ -441,21 +480,17 @@ void Draw_Plane_Surface(Surface * s)
if(CTX.geom.normals) {
List_Read(s->Orientations, 0, &vv1);
List_Read(s->Orientations, 1, &vv2);
- // don't rely on MeanPlane: teh orientation is arbotrary
- //n[0] = s->plan[2][0];
- //n[1] = s->plan[2][1];
- //n[2] = s->plan[2][2];
- //norme(n);
- Get_SurfaceNormal(s, n);
+ double x = (vv1.Pos.X + vv2.Pos.X) / 2.;
+ double y = (vv1.Pos.Y + vv2.Pos.Y) / 2.;
+ double z = (vv1.Pos.Z + vv2.Pos.Z) / 2.;
+ getPlaneSurfaceNormal(s, x, y, z, n);
n[0] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[0];
n[1] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[1];
n[2] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[2];
glColor4ubv((GLubyte *) & CTX.color.geom.normals);
Draw_Vector(CTX.vector_type, 0, CTX.arrow_rel_head_radius,
CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius,
- (vv2.Pos.X + vv1.Pos.X) / 2., (vv2.Pos.Y + vv1.Pos.Y) / 2.,
- (vv2.Pos.Z + vv1.Pos.Z) / 2., n[0], n[1], n[2], NULL,
- CTX.geom.light);
+ x, y, z, n[0], n[1], n[2], NULL, CTX.geom.light);
}
}
@@ -518,16 +553,25 @@ void Draw_NonPlane_Surface(Surface * s)
}
if(CTX.geom.normals) {
- Vertex v = InterpolateSurface(s, 0.5, 0.5, 0, 0);
- double n[3];
- Get_SurfaceNormal(s, n);
+ Vertex v1 = InterpolateSurface(s, 0.5, 0.5, 0, 0);
+ Vertex v2 = InterpolateSurface(s, 0.6, 0.5, 0, 0);
+ Vertex v3 = InterpolateSurface(s, 0.5, 0.6, 0, 0);
+ double t1[3], t2[3], n[3];
+ t1[0] = v2.Pos.X - v1.Pos.X;
+ t1[1] = v2.Pos.Y - v1.Pos.Y;
+ t1[2] = v2.Pos.Z - v1.Pos.Z;
+ t2[0] = v3.Pos.X - v1.Pos.X;
+ t2[1] = v3.Pos.Y - v1.Pos.Y;
+ t2[2] = v3.Pos.Z - v1.Pos.Z;
+ prodve(t1, t2, n);
+ norme(n);
n[0] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[0];
n[1] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[1];
n[2] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[2];
glColor4ubv((GLubyte *) & CTX.color.geom.normals);
Draw_Vector(CTX.vector_type, 0, CTX.arrow_rel_head_radius,
CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius,
- v.Pos.X, v.Pos.Y, v.Pos.Z, n[0], n[1], n[2], NULL,
+ v1.Pos.X, v1.Pos.Y, v1.Pos.Z, n[0], n[1], n[2], NULL,
CTX.geom.light);
}
}
diff --git a/Mesh/2D_Mesh.cpp b/Mesh/2D_Mesh.cpp
index e42f166d3b..b47cc4483d 100644
--- a/Mesh/2D_Mesh.cpp
+++ b/Mesh/2D_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: 2D_Mesh.cpp,v 1.64 2004-06-23 03:57:43 geuzaine Exp $
+// $Id: 2D_Mesh.cpp,v 1.65 2004-06-23 18:52:45 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -789,6 +789,8 @@ void ActionInvertQua(void *a, void *b)
q->V[3] = tmp;
}
+int isPointOnPlanarSurface(Surface * S, double X, double Y, double Z, double n[3]);
+
int isMiddlePointOnPlanarSurface(Surface *s, Vertex *v1, Vertex *v2)
{
double n[3] = {s->a, s->b, s->c};
@@ -856,50 +858,36 @@ void Get_SurfaceNormal(Surface *s, double n[3])
}
}
-// Horrible (tm) hack to orient the elements correctly. This is
-// *definitely* not the best way to do it, but I don't have time to
-// look into the issue right now.
+// A little routine to re-orient the surface meshe to match the
+// geometrical orientation. It would be better to actually generate
+// the meshes correctly in the first place (!), but this is handy
+// since we use many different algorithms since the mean plane
+// computation doesn't take the original orientation into account.
void ReOrientSurfaceMesh(Surface *s)
{
- double t1[3], t2[3], n1[3], n2[3], res;
- Simplex *simp;
- Quadrangle *quad;
-
- Get_SurfaceNormal(s, n1);
-
- if(Tree_Right(s->Simplexes, &simp)){
- t1[0] = simp->V[1]->Pos.X - simp->V[0]->Pos.X;
- t1[1] = simp->V[1]->Pos.Y - simp->V[0]->Pos.Y;
- t1[2] = simp->V[1]->Pos.Z - simp->V[0]->Pos.Z;
- t2[0] = simp->V[2]->Pos.X - simp->V[0]->Pos.X;
- t2[1] = simp->V[2]->Pos.Y - simp->V[0]->Pos.Y;
- t2[2] = simp->V[2]->Pos.Z - simp->V[0]->Pos.Z;
- prodve(t1, t2, n2);
- norme(n2);
- prosca(n1, n2, &res);
- if(res < 0.0){
- Msg(DEBUG, "Inverting triangles in %s surface %d (res = %g)",
- (s->Typ == MSH_SURF_PLAN) ? "Plane" : "NonPlane", s->Num, res);
+ Curve *c;
+ List_Read(s->Generatrices, 0, &c);
+ Vertex *v1, *v2;
+ List_Read(c->Vertices, 0, &v1);
+ List_Read(c->Vertices, 1, &v2);
+ Edge edge;
+ edge.V[0] = v1;
+ edge.V[1] = v2;
+
+ if(Tree_Nbr(s->Simplexes)){
+ EdgesContainer edges;
+ edges.AddSimplexTree(s->Simplexes);
+ Edge *e = (Edge*)Tree_PQuery(edges.AllEdges, &edge);
+ if(e && e->V[0]->Num == v2->Num && e->V[1]->Num == v1->Num)
Tree_Action(s->Simplexes, ActionInvertTri);
- }
}
-
- if(Tree_Right(s->Quadrangles, &quad)){
- t1[0] = quad->V[1]->Pos.X - quad->V[0]->Pos.X;
- t1[1] = quad->V[1]->Pos.Y - quad->V[0]->Pos.Y;
- t1[2] = quad->V[1]->Pos.Z - quad->V[0]->Pos.Z;
- t2[0] = quad->V[2]->Pos.X - quad->V[0]->Pos.X;
- t2[1] = quad->V[2]->Pos.Y - quad->V[0]->Pos.Y;
- t2[2] = quad->V[2]->Pos.Z - quad->V[0]->Pos.Z;
- prodve(t1, t2, n2);
- norme(n2);
- prosca(n1, n2, &res);
- if(res < 0.0){
- Msg(DEBUG, "Inverting quads in %s surface %d (res = %g)",
- (s->Typ == MSH_SURF_PLAN) ? "Plane" : "NonPlane", s->Num, res);
+ if(Tree_Nbr(s->Quadrangles)){
+ EdgesContainer edges;
+ edges.AddQuadrangleTree(s->Quadrangles);
+ Edge *e = (Edge*)Tree_PQuery(edges.AllEdges, &edge);
+ if(e && e->V[0]->Num == v2->Num && e->V[1]->Num == v1->Num)
Tree_Action(s->Quadrangles, ActionInvertQua);
- }
}
}
@@ -1006,7 +994,6 @@ void Maillage_Surface(void *data, void *dum)
End_Surface(s->Support, 0);
End_Surface(s, 0);
- // FIXME: big hack
ReOrientSurfaceMesh(s);
}
diff --git a/Mesh/3D_Extrude.cpp b/Mesh/3D_Extrude.cpp
index f98b655dd7..c0ab9572f1 100644
--- a/Mesh/3D_Extrude.cpp
+++ b/Mesh/3D_Extrude.cpp
@@ -1,4 +1,4 @@
-// $Id: 3D_Extrude.cpp,v 1.81 2004-06-23 03:57:43 geuzaine Exp $
+// $Id: 3D_Extrude.cpp,v 1.82 2004-06-23 18:52:45 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -1046,7 +1046,6 @@ int Extrude_Mesh(Surface * s)
Tree_Action(s->Simplexes, AddSimVertsInSurf);
Tree_Action(s->Quadrangles, AddQuadVertsInSurf);
- // FIXME: big hack
ReOrientSurfaceMesh(s);
return true;
diff --git a/Mesh/Mesh.h b/Mesh/Mesh.h
index 181cc9c19e..224bc57c86 100644
--- a/Mesh/Mesh.h
+++ b/Mesh/Mesh.h
@@ -476,8 +476,6 @@ void Projette_Plan_Moyen(void *a, void *b);
void Projette_Inverse(void *a, void *b);
void Freeze_Vertex(void *a, void *b);
void deFreeze_Vertex(void *a, void *b);
-int isPointOnPlanarSurface(Surface * S, double X, double Y, double Z, double n[3]);
-void Get_SurfaceNormal(Surface *s, double n[3]);
void ReOrientSurfaceMesh(Surface *s);
double Lc_XYZ(double X, double Y, double Z, Mesh * m);
diff --git a/doc/VERSIONS b/doc/VERSIONS
index 3cf8fe1b3d..8552079b24 100644
--- a/doc/VERSIONS
+++ b/doc/VERSIONS
@@ -1,12 +1,12 @@
-$Id: VERSIONS,v 1.226 2004-06-23 04:02:38 geuzaine Exp $
+$Id: VERSIONS,v 1.227 2004-06-23 18:52:45 geuzaine Exp $
New since 1.53: fixed UNV output; make Layers' region numbering
consistent between lines/surfaces/volumes; fixed home directory
problem on Win98; new Plugin(CutParametric); the default project file
is now created in the home directory if no current directory is
defined (e.g. when double-clicking on the icon on Windows/MacOS);
-(partial) fix of the geometry/mesh surface orientation discrepancies;
-small bug fixes and cleanups.
+fixed the discrepancy between the orientation of geometrical surfaces
+and the associated surface meshes; small bug fixes and cleanups.
New in 1.53: completed support for second order elements in the mesh
module (lines, triangles, quadrangles, tetrahedra, hexahedra, prisms
--
GitLab