diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 9fc9230f32565f1d8a0a858d528916d48c39ec7f..2dbc76443f2c08a5b36ca6608b01fc290727e895 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -848,9 +848,9 @@ StringXNumber MeshOptions_Number[] = {
// "Target quality for tetrahedral elements (not fully functional)" },
{ F|O, "RadiusInf" , opt_mesh_radius_inf , 0.0 ,
- "Only display elements whose Radius is greater than RadiusInf" },
+ "Only display elements whose longest edge is greater than RadiusInf" },
{ F|O, "RadiusSup" , opt_mesh_radius_sup , 0.0 ,
- "Only display elements whose Radius is smaller than RadiusSup" },
+ "Only display elements whose longest edge is smaller than RadiusSup" },
{ F|O, "RandomFactor" , opt_mesh_rand_factor , 1.e-4 ,
"Random factor used in 2D and 3D meshing algorithm (test other values when the algorithm fails)" },
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index 06f9c1d5322dc24947028fc370309d65cf8a4bb1..dfe72994fa385d0b622916f825dd1641adefae19 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.382 2004-11-18 23:44:53 geuzaine Exp $
+// $Id: GUI.cpp,v 1.383 2004-11-19 18:26:47 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -2087,7 +2087,7 @@ void GUI::create_option_window()
mesh_value[6] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 8 * BH, IW / 2, BH);
mesh_value[6]->align(FL_ALIGN_RIGHT);
- mesh_value[7] = new Fl_Value_Input(L + 2 * WB + IW / 2, 2 * WB + 8 * BH, IW / 2, BH, "Tetrahedra size range");
+ mesh_value[7] = new Fl_Value_Input(L + 2 * WB + IW / 2, 2 * WB + 8 * BH, IW / 2, BH, "Element size range");
mesh_value[7]->align(FL_ALIGN_RIGHT);
mesh_value[8] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 9 * BH, IW, BH, "Normals");
diff --git a/Graphics/Mesh.cpp b/Graphics/Mesh.cpp
index 757e6f29f5eab50128822da77d4ba19c43977484..b375603d229344510c142e8f173ff44055882b54 100644
--- a/Graphics/Mesh.cpp
+++ b/Graphics/Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Mesh.cpp,v 1.112 2004-10-30 03:07:29 geuzaine Exp $
+// $Id: Mesh.cpp,v 1.113 2004-11-19 18:26:47 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -285,6 +285,7 @@ void Draw_Mesh_Volume(void *a, void *b)
v->QuadVertexArray = new VertexArray(4, 1000);
v->QuadVertexArray->fill = 1;
Tree_Action(v->Simplexes, Draw_Mesh_Tetrahedron);
+ Tree_Action(v->SimplexesBase, Draw_Mesh_Tetrahedron);
Tree_Action(v->Hexahedra, Draw_Mesh_Hexahedron);
Tree_Action(v->Prisms, Draw_Mesh_Prism);
Tree_Action(v->Pyramids, Draw_Mesh_Pyramid);
@@ -312,6 +313,7 @@ void Draw_Mesh_Volume(void *a, void *b)
CTX.mesh.normals){
Msg(DEBUG, "classic volume data path");
Tree_Action(v->Simplexes, Draw_Mesh_Tetrahedron);
+ Tree_Action(v->SimplexesBase, Draw_Mesh_Tetrahedron);
Tree_Action(v->Hexahedra, Draw_Mesh_Hexahedron);
Tree_Action(v->Prisms, Draw_Mesh_Prism);
Tree_Action(v->Pyramids, Draw_Mesh_Pyramid);
@@ -336,6 +338,7 @@ void Draw_Mesh_Surface(void *a, void *b)
s->normals = new smooth_normals;
preproNormals = 1;
Tree_Action(s->Simplexes, Draw_Mesh_Triangle);
+ Tree_Action(s->SimplexesBase, Draw_Mesh_Triangle);
Tree_Action(s->Quadrangles, Draw_Mesh_Quadrangle);
preproNormals = 0;
}
@@ -345,9 +348,11 @@ void Draw_Mesh_Surface(void *a, void *b)
Msg(DEBUG, "regenerate surface mesh vertex arrays");
// triangles
if(s->TriVertexArray) delete s->TriVertexArray;
- s->TriVertexArray = new VertexArray(3, Tree_Nbr(s->Simplexes));
+ s->TriVertexArray = new VertexArray(3, Tree_Nbr(s->Simplexes) +
+ Tree_Nbr(s->SimplexesBase));
s->TriVertexArray->fill = 1;
Tree_Action(s->Simplexes, Draw_Mesh_Triangle);
+ Tree_Action(s->SimplexesBase, Draw_Mesh_Triangle);
if(s->TriVertexArray){
Msg(DEBUG, "%d triangles in surface vertex array", s->TriVertexArray->num);
s->TriVertexArray->fill = 0;
@@ -374,6 +379,7 @@ void Draw_Mesh_Surface(void *a, void *b)
CTX.mesh.points_per_element || CTX.mesh.normals){
Msg(DEBUG, "classic triangle data path");
Tree_Action(s->Simplexes, Draw_Mesh_Triangle);
+ Tree_Action(s->SimplexesBase, Draw_Mesh_Triangle);
}
if(!s->QuadVertexArray || CTX.mesh.dual || CTX.mesh.surfaces_num ||
@@ -405,6 +411,7 @@ void Draw_Mesh_Curve(void *a, void *b)
thePhysical = getFirstPhysical(MSH_PHYSICAL_LINE, c->Num);
Tree_Action(c->Simplexes, Draw_Mesh_Line);
+ Tree_Action(c->SimplexesBase, Draw_Mesh_Line);
}
void Draw_Mesh_Point(int num, double x, double y, double z, int degree, int visible)
@@ -454,7 +461,7 @@ void Draw_Mesh_Line(void *a, void *b)
double Xc = 0.0, Yc = 0.0, Zc = 0.0, m[3];
char Num[100];
- Simplex *s = *(Simplex **) a;
+ SimplexBase *s = *(SimplexBase **) a;
if(!(s->Visible & VIS_MESH))
return;
@@ -463,6 +470,12 @@ void Draw_Mesh_Line(void *a, void *b)
if(part && !(*part)->Visible)
return;
+ if(CTX.mesh.radius_sup) {
+ double r = s->maxEdge();
+ if(r < CTX.mesh.radius_inf || r > CTX.mesh.radius_sup)
+ return;
+ }
+
if(!CTX.mesh.cut_plane_only_volume && intersectCutPlane(2, s->V) < 0)
return;
@@ -721,7 +734,7 @@ void Draw_Mesh_Triangle(void *a, void *b)
double n[3];
char Num[256];
- Simplex *s = *(Simplex **) a;
+ SimplexBase *s = *(SimplexBase **) a;
if(!s->V[2]) // the old extrusion algo puts lines in the simp_surf tree...
return;
@@ -733,6 +746,12 @@ void Draw_Mesh_Triangle(void *a, void *b)
if(part && !(*part)->Visible)
return;
+ if(CTX.mesh.radius_sup) {
+ double r = s->maxEdge();
+ if(r < CTX.mesh.radius_inf || r > CTX.mesh.radius_sup)
+ return;
+ }
+
if(!CTX.mesh.cut_plane_only_volume && intersectCutPlane(3, s->V) < 0)
return;
@@ -886,6 +905,12 @@ void Draw_Mesh_Quadrangle(void *a, void *b)
if(part && !(*part)->Visible)
return;
+ if(CTX.mesh.radius_sup) {
+ double r = q->maxEdge();
+ if(r < CTX.mesh.radius_inf || r > CTX.mesh.radius_sup)
+ return;
+ }
+
if(!CTX.mesh.cut_plane_only_volume && intersectCutPlane(4, q->V) < 0)
return;
@@ -1032,7 +1057,7 @@ void Draw_Mesh_Tetrahedron(void *a, void *b)
char Num[100];
double tmp, X[4], Y[4], Z[4], X2[6], Y2[6], Z2[6];
- Simplex *s = *(Simplex **) a;
+ SimplexBase *s = *(SimplexBase **) a;
if(!s->V[3] || !(s->Visible & VIS_MESH))
return;
@@ -1048,10 +1073,11 @@ void Draw_Mesh_Tetrahedron(void *a, void *b)
}
if(CTX.mesh.radius_sup) {
- if(s->Radius < CTX.mesh.radius_inf || s->Radius > CTX.mesh.radius_sup)
+ double r = s->maxEdge();
+ if(r < CTX.mesh.radius_inf || r > CTX.mesh.radius_sup)
return;
}
-
+
int edges = CTX.mesh.volumes_edges;
int faces = CTX.mesh.volumes_faces;
@@ -1213,6 +1239,12 @@ void Draw_Mesh_Hexahedron(void *a, void *b)
if(part && !(*part)->Visible)
return;
+ if(CTX.mesh.radius_sup) {
+ double r = h->maxEdge();
+ if(r < CTX.mesh.radius_inf || r > CTX.mesh.radius_sup)
+ return;
+ }
+
int edges = CTX.mesh.volumes_edges;
int faces = CTX.mesh.volumes_faces;
@@ -1388,6 +1420,12 @@ void Draw_Mesh_Prism(void *a, void *b)
if(part && !(*part)->Visible)
return;
+ if(CTX.mesh.radius_sup) {
+ double r = p->maxEdge();
+ if(r < CTX.mesh.radius_inf || r > CTX.mesh.radius_sup)
+ return;
+ }
+
int edges = CTX.mesh.volumes_edges;
int faces = CTX.mesh.volumes_faces;
@@ -1578,6 +1616,12 @@ void Draw_Mesh_Pyramid(void *a, void *b)
if(part && !(*part)->Visible)
return;
+ if(CTX.mesh.radius_sup) {
+ double r = p->maxEdge();
+ if(r < CTX.mesh.radius_inf || r > CTX.mesh.radius_sup)
+ return;
+ }
+
int edges = CTX.mesh.volumes_edges;
int faces = CTX.mesh.volumes_faces;
diff --git a/Makefile b/Makefile
index d08750d2c7da871b3f0f97b75064b709d6b677c0..a0b4ba0f326e03fdd166f77cd13e1c113216955c 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.379 2004-11-14 04:38:45 geuzaine Exp $
+# $Id: Makefile,v 1.380 2004-11-19 18:26:46 geuzaine Exp $
#
# Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
#
@@ -22,8 +22,8 @@
include variables
GMSH_MAJOR_VERSION = 1
-GMSH_MINOR_VERSION = 56
-GMSH_PATCH_VERSION = 4
+GMSH_MINOR_VERSION = 57
+GMSH_PATCH_VERSION = 0
GMSH_EXTRA_VERSION = "-cvs"
GMSH_VERSION = ${GMSH_MAJOR_VERSION}.${GMSH_MINOR_VERSION}.${GMSH_PATCH_VERSION}${GMSH_EXTRA_VERSION}
diff --git a/Mesh/Create.cpp b/Mesh/Create.cpp
index 9fdd5cfd77e2d6d87ec9d4655e40130b7f2c5d52..b68adf100b224063c78404a2c3758af3cd103fd6 100644
--- a/Mesh/Create.cpp
+++ b/Mesh/Create.cpp
@@ -1,4 +1,4 @@
-// $Id: Create.cpp,v 1.63 2004-08-12 16:57:32 geuzaine Exp $
+// $Id: Create.cpp,v 1.64 2004-11-19 18:26:47 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -534,6 +534,7 @@ Curve *Create_Curve(int Num, int Typ, int Order, List_T * Liste,
pC->Num = Num;
THEM->MaxLineNum = IMAX(THEM->MaxLineNum, Num);
pC->Simplexes = Tree_Create(sizeof(Simplex *), compareSimplex);
+ pC->SimplexesBase = Tree_Create(sizeof(SimplexBase *), compareSimplexBase);
pC->Method = LIBRE;
pC->degre = Order;
pC->Circle.n[0] = 0.0;
@@ -630,6 +631,8 @@ void Free_Curve(void *a, void *b)
List_Delete(pC->Vertices);
Tree_Action(pC->Simplexes, Free_Simplex);
Tree_Delete(pC->Simplexes);
+ Tree_Action(pC->SimplexesBase, Free_SimplexBase);
+ Tree_Delete(pC->SimplexesBase);
Free(pC->k);
List_Delete(pC->Control_Points);
Free(pC->cp);
@@ -658,6 +661,7 @@ Surface *Create_Surface(int Num, int Typ)
pS->RecombineAngle = 75;
pS->TrsfPoints = List_Create(4, 4, sizeof(Vertex *));
pS->Simplexes = Tree_Create(sizeof(Simplex *), compareQuality);
+ pS->SimplexesBase = Tree_Create(sizeof(SimplexBase *), compareSimplexBase);
pS->Quadrangles = Tree_Create(sizeof(Quadrangle *), compareQuadrangle);
pS->Vertices = Tree_Create(sizeof(Vertex *), compareVertex);
pS->TrsfVertices = List_Create(1, 10, sizeof(Vertex *));
@@ -682,6 +686,8 @@ void Free_Surface(void *a, void *b)
List_Delete(pS->TrsfPoints);
Tree_Action(pS->Simplexes, Free_Simplex);
Tree_Delete(pS->Simplexes);
+ Tree_Action(pS->SimplexesBase, Free_SimplexBase);
+ Tree_Delete(pS->SimplexesBase);
Tree_Action(pS->Quadrangles, Free_Quadrangle);
Tree_Delete(pS->Quadrangles);
Tree_Delete(pS->Vertices);
@@ -724,6 +730,7 @@ Volume *Create_Volume(int Num, int Typ)
pV->Surfaces = List_Create(1, 2, sizeof(Surface *));
pV->SurfacesOrientations = List_Create(1, 2, sizeof(int));
pV->Simplexes = Tree_Create(sizeof(Simplex *), compareQuality);
+ pV->SimplexesBase = Tree_Create(sizeof(Simplex *), compareSimplexBase);
pV->Vertices = Tree_Create(sizeof(Vertex *), compareVertex);
pV->Hexahedra = Tree_Create(sizeof(Hexahedron *), compareHexahedron);
pV->Prisms = Tree_Create(sizeof(Prism *), comparePrism);
@@ -744,6 +751,7 @@ void Free_Volume(void *a, void *b)
Volume *pV = *(Volume **) a;
if(pV) {
Tree_Action(pV->Simplexes, Free_Simplex);
+ Tree_Action(pV->SimplexesBase, Free_SimplexBase);
Tree_Action(pV->Hexahedra, Free_Hexahedron);
Tree_Action(pV->Prisms, Free_Prism);
Tree_Action(pV->Pyramids, Free_Pyramid);
@@ -760,6 +768,7 @@ void Free_Volume_But_Not_Elements(void *a, void *b)
List_Delete(pV->Surfaces); // surfaces freed elsewhere
List_Delete(pV->SurfacesOrientations);
Tree_Delete(pV->Simplexes);
+ Tree_Delete(pV->SimplexesBase);
Tree_Delete(pV->Simp_Surf); // for old extrusion mesh generator
Tree_Delete(pV->Quad_Surf); // for old extrusion mesh generator
Tree_Delete(pV->Vertices); // vertices freed elsewhere
diff --git a/Mesh/Element.cpp b/Mesh/Element.cpp
index 11b7f96eb786f1af4a29e1413c1488d5b75e9ad0..bffc464fe35e30173150813be1b9c25d78a42fad 100644
--- a/Mesh/Element.cpp
+++ b/Mesh/Element.cpp
@@ -1,4 +1,4 @@
-// $Id: Element.cpp,v 1.2 2004-07-21 22:19:56 geuzaine Exp $
+// $Id: Element.cpp,v 1.3 2004-11-19 18:26:47 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -24,48 +24,56 @@
#include "Numeric.h"
#include "Element.h"
+extern int edges_quad[4][2];
+extern int edges_hexa[12][2];
+extern int edges_prism[9][2];
+extern int edges_pyramid[8][2];
+
int Element::TotalNumber = 0;
Element::Element()
- : Num(0), iEnt(-1), iPart(-1), Visible(VIS_MESH), VSUP(NULL)
+ : iEnt(-1), iPart(-1), Visible(VIS_MESH), VSUP(NULL)
{
+ Num = ++TotalNumber;
}
Element::~Element()
{
+ if(VSUP) Free(VSUP);
+}
+
+double Element::lij(Vertex *Vi, Vertex *Vj)
+{
+ return sqrt(DSQR(Vi->Pos.X - Vj->Pos.X) +
+ DSQR(Vi->Pos.Y - Vj->Pos.Y) +
+ DSQR(Vi->Pos.Z - Vj->Pos.Z));
}
// Quads
Quadrangle::Quadrangle()
+ : Element()
{
- iEnt = -1;
- iPart = -1;
- Num = ++TotalNumber;
- Visible = VIS_MESH;
for(int i = 0; i < 4; i++) V[i] = NULL;
- VSUP = NULL;
}
Quadrangle::Quadrangle(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4)
+ : Element()
{
- iEnt = -1;
- iPart = -1;
- Num = ++TotalNumber;
- Visible = VIS_MESH;
V[0] = v1; V[1] = v2; V[2] = v3; V[3] = v4;
- VSUP = NULL;
}
-Quadrangle::~Quadrangle()
+double Quadrangle::maxEdge()
{
- if(VSUP) Free(VSUP);
+ double maxlij = 0.;
+ for(int i = 0; i < 4; i++)
+ maxlij = DMAX(maxlij, lij(V[edges_quad[i][0]], V[edges_quad[i][1]]));
+ return maxlij;
}
Quadrangle *Create_Quadrangle(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4)
{
- Quadrangle *q = new Quadrangle(v1, v2, v3, v4);
- return q;
+ return new Quadrangle(v1, v2, v3, v4);
}
void Free_Quadrangle(void *a, void *b)
@@ -87,28 +95,21 @@ int compareQuadrangle(const void *a, const void *b)
// Hexas
Hexahedron::Hexahedron()
+ : Element()
{
- iEnt = -1;
- iPart = -1;
- Num = ++TotalNumber;
- Visible = VIS_MESH;
for(int i = 0; i < 8; i++) V[i] = NULL;
- VSUP = NULL;
}
Hexahedron::Hexahedron(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4,
Vertex *v5, Vertex *v6, Vertex *v7, Vertex *v8)
+ : Element()
{
- iEnt = -1;
- iPart = -1;
- Num = ++TotalNumber;
- Visible = VIS_MESH;
V[0] = v1; V[1] = v2; V[2] = v3; V[3] = v4;
V[4] = v5; V[5] = v6; V[6] = v7; V[7] = v8;
- VSUP = NULL;
}
-double Hexahedron::Orientation(){
+double Hexahedron::Orientation()
+{
double mat[3][3];
mat[0][0] = V[1]->Pos.X - V[0]->Pos.X;
mat[0][1] = V[3]->Pos.X - V[0]->Pos.X;
@@ -122,17 +123,19 @@ double Hexahedron::Orientation(){
return det3x3(mat);
}
-Hexahedron::~Hexahedron()
+double Hexahedron::maxEdge()
{
- if(VSUP) Free(VSUP);
+ double maxlij = 0.;
+ for(int i = 0; i < 12; i++)
+ maxlij = DMAX(maxlij, lij(V[edges_hexa[i][0]], V[edges_hexa[i][1]]));
+ return maxlij;
}
Hexahedron *Create_Hexahedron(Vertex * v1, Vertex * v2, Vertex * v3,
Vertex * v4, Vertex * v5, Vertex * v6,
Vertex * v7, Vertex * v8)
{
- Hexahedron *h = new Hexahedron(v1, v2, v3, v4, v5, v6, v7, v8);
- return h;
+ return new Hexahedron(v1, v2, v3, v4, v5, v6, v7, v8);
}
void Free_Hexahedron(void *a, void *b)
@@ -154,28 +157,21 @@ int compareHexahedron(const void *a, const void *b)
// Prisms
Prism::Prism()
+ : Element()
{
- iEnt = -1;
- iPart = -1;
- Num = 0;
- Visible = VIS_MESH;
for(int i = 0; i < 6; i++) V[i] = NULL;
- VSUP = NULL;
}
Prism::Prism(Vertex *v1, Vertex *v2, Vertex *v3,
Vertex *v4, Vertex *v5, Vertex *v6)
+ : Element()
{
- iEnt = -1;
- iPart = -1;
- Num = ++TotalNumber;
- Visible = VIS_MESH;
V[0] = v1; V[1] = v2; V[2] = v3;
V[3] = v4; V[4] = v5; V[5] = v6;
- VSUP = NULL;
}
-double Prism::Orientation(){
+double Prism::Orientation()
+{
double mat[3][3];
mat[0][0] = V[1]->Pos.X - V[0]->Pos.X;
mat[0][1] = V[2]->Pos.X - V[0]->Pos.X;
@@ -189,16 +185,18 @@ double Prism::Orientation(){
return det3x3(mat);
}
-Prism::~Prism()
+double Prism::maxEdge()
{
- if(VSUP) Free(VSUP);
+ double maxlij = 0.;
+ for(int i = 0; i < 9; i++)
+ maxlij = DMAX(maxlij, lij(V[edges_prism[i][0]], V[edges_prism[i][1]]));
+ return maxlij;
}
Prism *Create_Prism(Vertex * v1, Vertex * v2, Vertex * v3,
Vertex * v4, Vertex * v5, Vertex * v6)
{
- Prism *p = new Prism(v1, v2, v3, v4, v5, v6);
- return p;
+ return new Prism(v1, v2, v3, v4, v5, v6);
}
void Free_Prism(void *a, void *b)
@@ -220,16 +218,19 @@ int comparePrism(const void *a, const void *b)
// Pyramids
Pyramid::Pyramid()
+ : Element()
{
- iEnt = -1;
- iPart = -1;
- Num = 0;
- Visible = VIS_MESH;
for(int i = 0; i < 5; i++) V[i] = NULL;
- VSUP = NULL;
}
-double Pyramid::Orientation(){
+Pyramid::Pyramid(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4, Vertex *v5)
+ : Element()
+{
+ V[0] = v1; V[1] = v2; V[2] = v3; V[3] = v4; V[4] = v5;
+}
+
+double Pyramid::Orientation()
+{
double mat[3][3];
mat[0][0] = V[1]->Pos.X - V[0]->Pos.X;
mat[0][1] = V[3]->Pos.X - V[0]->Pos.X;
@@ -243,26 +244,18 @@ double Pyramid::Orientation(){
return det3x3(mat);
}
-Pyramid::Pyramid(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4, Vertex *v5)
+double Pyramid::maxEdge()
{
- iEnt = -1;
- iPart = -1;
- Num = ++TotalNumber;
- Visible = VIS_MESH;
- V[0] = v1; V[1] = v2; V[2] = v3; V[3] = v4; V[4] = v5;
- VSUP = NULL;
-}
-
-Pyramid::~Pyramid()
-{
- if(VSUP) Free(VSUP);
+ double maxlij = 0.;
+ for(int i = 0; i < 8; i++)
+ maxlij = DMAX(maxlij, lij(V[edges_pyramid[i][0]], V[edges_pyramid[i][1]]));
+ return maxlij;
}
Pyramid *Create_Pyramid(Vertex * v1, Vertex * v2, Vertex * v3,
Vertex * v4, Vertex * v5)
{
- Pyramid *p = new Pyramid(v1, v2, v3, v4, v5);
- return p;
+ return new Pyramid(v1, v2, v3, v4, v5);
}
void Free_Pyramid(void *a, void *b)
diff --git a/Mesh/Element.h b/Mesh/Element.h
index 8bfe99f042ff007f0dec0f7061f6c861b1c276cc..20dafdaac964c4cac3fd258951db20e1aed31d7e 100644
--- a/Mesh/Element.h
+++ b/Mesh/Element.h
@@ -32,6 +32,7 @@ class Element {
static int TotalNumber;
Element();
virtual ~Element();
+ double lij(Vertex *Vi, Vertex *Vj);
};
class Quadrangle : public Element{
@@ -39,7 +40,8 @@ class Quadrangle : public Element{
Vertex *V[4];
Quadrangle();
Quadrangle(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4);
- ~Quadrangle();
+ ~Quadrangle() {}
+ double maxEdge();
};
class Hexahedron : public Element{
@@ -48,8 +50,9 @@ class Hexahedron : public Element{
Hexahedron();
Hexahedron(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4,
Vertex *v5, Vertex *v6, Vertex *v7, Vertex *v8);
- ~Hexahedron();
+ ~Hexahedron() {}
double Orientation();
+ double maxEdge();
};
class Prism : public Element{
@@ -58,8 +61,9 @@ class Prism : public Element{
Prism();
Prism(Vertex *v1, Vertex *v2, Vertex *v3,
Vertex *v4, Vertex *v5, Vertex *v6);
- ~Prism();
+ ~Prism() {}
double Orientation();
+ double maxEdge();
};
class Pyramid : public Element{
@@ -67,8 +71,9 @@ class Pyramid : public Element{
Vertex *V[5];
Pyramid();
Pyramid(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4, Vertex *v5);
- ~Pyramid();
+ ~Pyramid() {}
double Orientation();
+ double maxEdge();
};
// C interface
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index f92199c5f9e5578075b7a6a888a1a10582370f73..16789b5004d77576e48f8ad912f7d266861eceba 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.60 2004-08-13 20:59:44 geuzaine Exp $
+// $Id: Generator.cpp,v 1.61 2004-11-19 18:26:47 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -69,7 +69,7 @@ void GetStatistics(double stat[50])
Surface *s;
List_Read(surfaces, i, &s);
stat[5] += Tree_Nbr(s->Vertices);
- stat[7] += Tree_Nbr(s->Simplexes);
+ stat[7] += Tree_Nbr(s->Simplexes) + Tree_Nbr(s->SimplexesBase);
stat[8] += Tree_Nbr(s->Quadrangles);
}
List_Delete(surfaces);
@@ -82,7 +82,7 @@ void GetStatistics(double stat[50])
Volume *v;
List_Read(volumes, i, &v);
stat[6] += Tree_Nbr(v->Vertices);
- stat[9] += Tree_Nbr(v->Simplexes);
+ stat[9] += Tree_Nbr(v->Simplexes) + Tree_Nbr(v->SimplexesBase);
stat[10] += Tree_Nbr(v->Hexahedra);
stat[11] += Tree_Nbr(v->Prisms);
stat[12] += Tree_Nbr(v->Pyramids);
diff --git a/Mesh/Mesh.h b/Mesh/Mesh.h
index 35bc8b95a2f71ec42a14a50af2b63c7dcaa49656..c05895ea3446fe106c1b2685161475af049b854c 100644
--- a/Mesh/Mesh.h
+++ b/Mesh/Mesh.h
@@ -242,7 +242,7 @@ struct _Surf{
double a, b, c, d;
List_T *Orientations;
List_T *Contours;
- Tree_T *Simplexes;
+ Tree_T *Simplexes, *SimplexesBase;
Tree_T *Quadrangles;
Tree_T *Vertices;
List_T *TrsfVertices;
@@ -307,7 +307,7 @@ typedef struct {
Tree_T *Vertices;
Tree_T *Edges;
Tree_T *Faces;
- Tree_T *Simplexes;
+ Tree_T *Simplexes, *SimplexesBase;
Tree_T *Simp_Surf; // for old extrusion mesh generator
Tree_T *Quad_Surf; // for old extrusion mesh generator
Tree_T *Hexahedra;
@@ -371,7 +371,7 @@ typedef struct{
double ubeg, uend;
List_T *Control_Points;
List_T *Vertices;
- Tree_T *Simplexes;
+ Tree_T *Simplexes, *SimplexesBase;
ExtrudeParams *Extrude;
float *k, *cp;
int degre;
diff --git a/Mesh/MeshQuality.cpp b/Mesh/MeshQuality.cpp
index 1025d88d75fba1b17257d5597624c84a9ab5039e..822188b3431a2629b91816d0fb4399737b64e5e3 100644
--- a/Mesh/MeshQuality.cpp
+++ b/Mesh/MeshQuality.cpp
@@ -1,4 +1,4 @@
-// $Id: MeshQuality.cpp,v 1.14 2004-05-25 04:10:05 geuzaine Exp $
+// $Id: MeshQuality.cpp,v 1.15 2004-11-19 18:26:47 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -37,7 +37,7 @@ static int NbCalcGamma, *Histogram;
void CalculateGamma(void *a, void *b)
{
- Simplex *s = *(Simplex **) a;
+ SimplexBase *s = *(SimplexBase **) a;
double gamma = s->GammaShapeMeasure();
NbCalcGamma++;
GAMMAMAX = DMAX(GAMMAMAX, gamma);
@@ -51,7 +51,7 @@ void CalculateGamma(void *a, void *b)
void CalculateEta(void *a, void *b)
{
- Simplex *s = *(Simplex **) a;
+ SimplexBase *s = *(SimplexBase **) a;
double gamma = s->EtaShapeMeasure();
NbCalcGamma++;
GAMMAMAX = DMAX(GAMMAMAX, gamma);
@@ -65,7 +65,7 @@ void CalculateEta(void *a, void *b)
void CalculateR(void *a, void *b)
{
- Simplex *s = *(Simplex **) a;
+ SimplexBase *s = *(SimplexBase **) a;
double gamma = s->RhoShapeMeasure();
NbCalcGamma++;
GAMMAMAX = DMAX(GAMMAMAX, gamma);
@@ -77,12 +77,11 @@ void CalculateR(void *a, void *b)
Histogram[i]++;
}
-
-
static void g(void *a, void *b)
{
Volume *v = *(Volume **) a;
Tree_Action(v->Simplexes, CalculateGamma);
+ Tree_Action(v->SimplexesBase, CalculateGamma);
Msg(DEBUG, "Gamma computed in volume %d (%d values)", v->Num, NbCalcGamma);
}
@@ -108,6 +107,7 @@ static void e(void *a, void *b)
{
Volume *v = *(Volume **) a;
Tree_Action(v->Simplexes, CalculateEta);
+ Tree_Action(v->SimplexesBase, CalculateEta);
Msg(DEBUG, "Eta computed in volume %d (%d values)", v->Num, NbCalcGamma);
}
@@ -132,6 +132,7 @@ static void r(void *a, void *b)
{
Volume *v = *(Volume **) a;
Tree_Action(v->Simplexes, CalculateR);
+ Tree_Action(v->SimplexesBase, CalculateR);
Msg(DEBUG, "Rho computed in volume %d (%d values)", v->Num, NbCalcGamma);
}
diff --git a/Mesh/Print_Mesh.cpp b/Mesh/Print_Mesh.cpp
index 2b81cff4963668405cdd37336b304830c1ada264..f0a746b85416395c3dd32bdda214899941b1e09d 100644
--- a/Mesh/Print_Mesh.cpp
+++ b/Mesh/Print_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Print_Mesh.cpp,v 1.55 2004-06-08 02:10:32 geuzaine Exp $
+// $Id: Print_Mesh.cpp,v 1.56 2004-11-19 18:26:47 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -102,7 +102,7 @@ static void _msh_print_simplex(void *a, void *b)
{
int i, type, nbn, nbs = 0;
- Simplex *s = *(Simplex **) a;
+ SimplexBase *s = *(SimplexBase **) a;
if(MSH_VOL_NUM && (MSH_VOL_NUM != s->iEnt))
return;
@@ -428,6 +428,7 @@ static void _msh_print_elements(Mesh *M)
MSH_LIN_NUM = abs(Num);
MSH_PHYSICAL_ORI = sign(Num);
Tree_Action(pc->Simplexes, _msh_print_simplex);
+ Tree_Action(pc->SimplexesBase, _msh_print_simplex);
}
}
}
@@ -454,6 +455,7 @@ static void _msh_print_elements(Mesh *M)
MSH_SUR_NUM = abs(Num);
MSH_PHYSICAL_ORI = sign(Num);
Tree_Action(ps->Simplexes, _msh_print_simplex);
+ Tree_Action(ps->SimplexesBase, _msh_print_simplex);
Tree_Action(ps->Quadrangles, _msh_print_quadrangle);
}
}
@@ -468,6 +470,7 @@ static void _msh_print_elements(Mesh *M)
MSH_VOL_NUM = abs(Num);
MSH_PHYSICAL_ORI = sign(Num);
Tree_Action(pV->Simplexes, _msh_print_simplex);
+ Tree_Action(pV->SimplexesBase, _msh_print_simplex);
Tree_Action(pV->Hexahedra, _msh_print_hexahedron);
Tree_Action(pV->Prisms, _msh_print_prism);
Tree_Action(pV->Pyramids, _msh_print_pyramid);
@@ -491,12 +494,14 @@ static void _msh_print_all_curves(void *a, void *b)
{
Curve *c = *(Curve **) a;
Tree_Action(c->Simplexes, _msh_print_simplex);
+ Tree_Action(c->SimplexesBase, _msh_print_simplex);
}
static void _msh_print_all_surfaces(void *a, void *b)
{
Surface *s = *(Surface **) a;
Tree_Action(s->Simplexes, _msh_print_simplex);
+ Tree_Action(s->SimplexesBase, _msh_print_simplex);
Tree_Action(s->Quadrangles, _msh_print_quadrangle);
}
@@ -511,6 +516,7 @@ static void _msh_print_all_volumes(void *a, void *b)
{
Volume *v = *(Volume **) a;
Tree_Action(v->Simplexes, _msh_print_simplex);
+ Tree_Action(v->SimplexesBase, _msh_print_simplex);
Tree_Action(v->Hexahedra, _msh_print_hexahedron);
Tree_Action(v->Prisms, _msh_print_prism);
Tree_Action(v->Pyramids, _msh_print_pyramid);
@@ -685,10 +691,10 @@ static void _unv_process_1D_elements(Mesh *m)
List_Read(ListSurfaces, i, &surf);
for(int j = 0; j < List_Nbr(surf->Generatrices); j++) {
List_Read(surf->Generatrices, j, &c);
- if(Tree_Nbr(c->Simplexes))
+ if(Tree_Nbr(c->Simplexes) || Tree_Nbr(c->SimplexesBase))
List_Add(AllCurves, &c);
c = FindCurve(-c->Num, m);
- if(Tree_Nbr(c->Simplexes))
+ if(Tree_Nbr(c->Simplexes) || Tree_Nbr(c->SimplexesBase))
List_Add(AllCurves, &c);
}
}
@@ -696,15 +702,17 @@ static void _unv_process_1D_elements(Mesh *m)
for(int i = 0; i < List_Nbr(ListCurves); i++) {
List_Read(ListCurves, i, &c);
if(!List_Search(AllCurves, &c, compareCurve)) {
- Elements = Tree2List(c->Simplexes);
- for(int j = 0; j < List_Nbr(Elements); j++) {
- List_Read(Elements, j, &sx);
- if(sx->VSUP)
- _unv_print_record(sx->Num, BEAM2, c->Num, 2, 2, &sx->V[0], sx->VSUP);
- else
- _unv_print_record(sx->Num, BEAM, c->Num, 2, 0, &sx->V[0], NULL);
+ for(int simtype = 0; simtype < 2; simtype ++){
+ Elements = (!simtype) ? Tree2List(c->Simplexes) : Tree2List(c->SimplexesBase);
+ for(int j = 0; j < List_Nbr(Elements); j++) {
+ List_Read(Elements, j, &sx);
+ if(sx->VSUP)
+ _unv_print_record(sx->Num, BEAM2, c->Num, 2, 2, &sx->V[0], sx->VSUP);
+ else
+ _unv_print_record(sx->Num, BEAM, c->Num, 2, 0, &sx->V[0], NULL);
+ }
+ List_Delete(Elements);
}
- List_Delete(Elements);
}
}
@@ -728,7 +736,8 @@ static void _unv_process_2D_elements(Mesh *m)
List_Read(ListVolumes, i, &vol);
for(int j = 0; j < List_Nbr(vol->Surfaces); j++) {
List_Read(vol->Surfaces, j, &s);
- if(Tree_Nbr(s->Simplexes) || Tree_Nbr(s->Quadrangles))
+ if(Tree_Nbr(s->Simplexes) || Tree_Nbr(s->SimplexesBase) ||
+ Tree_Nbr(s->Quadrangles))
List_Add(AllSurfaces, &s);
}
}
@@ -738,16 +747,18 @@ static void _unv_process_2D_elements(Mesh *m)
if(!List_Search(AllSurfaces, &s, compareSurface)) {
// triangles
- Elements = Tree2List(s->Simplexes);
- for(int j = 0; j < List_Nbr(Elements); j++) {
- List_Read(Elements, j, &sx);
- if(sx->VSUP)
- _unv_print_record(abs(sx->Num), THINSHLL, s->Num, 3, 3, &sx->V[0], sx->VSUP);
- else
- _unv_print_record(abs(sx->Num), THINSHLL, s->Num, 3, 0, &sx->V[0], NULL);
+ for(int simtype = 0; simtype < 2; simtype++){
+ Elements = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
+ for(int j = 0; j < List_Nbr(Elements); j++) {
+ List_Read(Elements, j, &sx);
+ if(sx->VSUP)
+ _unv_print_record(abs(sx->Num), THINSHLL, s->Num, 3, 3, &sx->V[0], sx->VSUP);
+ else
+ _unv_print_record(abs(sx->Num), THINSHLL, s->Num, 3, 0, &sx->V[0], NULL);
+ }
+ List_Delete(Elements);
}
- List_Delete(Elements);
-
+
// quadrangles
Elements = Tree2List(s->Quadrangles);
for(int j = 0; j < List_Nbr(Elements); j++) {
@@ -778,29 +789,31 @@ static void _unv_process_3D_elements(Mesh *m)
List_Read(ListVolumes, i, &v);
// Tets
- Elements = Tree2List(v->Simplexes);
- for(int j = 0; j < List_Nbr(Elements); j++) {
- List_Read(Elements, j, &sx);
- if(sx->Volume_Simplexe() < 0) {
- Vertex *temp;
- temp = sx->V[0];
- sx->V[0] = sx->V[1];
- sx->V[1] = temp;
- if(sx->VSUP){
- temp = sx->VSUP[1];
- sx->VSUP[1] = sx->VSUP[2];
- sx->VSUP[2] = temp;
- temp = sx->VSUP[5];
- sx->VSUP[5] = sx->VSUP[3];
- sx->VSUP[3] = temp;
+ for(int simtype = 0; simtype < 2; simtype++){
+ Elements = (!simtype) ? Tree2List(v->Simplexes) : Tree2List(v->SimplexesBase);
+ for(int j = 0; j < List_Nbr(Elements); j++) {
+ List_Read(Elements, j, &sx);
+ if(sx->Volume_Simplexe() < 0) {
+ Vertex *temp;
+ temp = sx->V[0];
+ sx->V[0] = sx->V[1];
+ sx->V[1] = temp;
+ if(sx->VSUP){
+ temp = sx->VSUP[1];
+ sx->VSUP[1] = sx->VSUP[2];
+ sx->VSUP[2] = temp;
+ temp = sx->VSUP[5];
+ sx->VSUP[5] = sx->VSUP[3];
+ sx->VSUP[3] = temp;
+ }
}
+ if(sx->VSUP)
+ _unv_print_record(ELEMENT_ID++, SOLIDFEM, v->Num, 4, 6, &sx->V[0], sx->VSUP);
+ else
+ _unv_print_record(ELEMENT_ID++, SOLIDFEM, v->Num, 4, 0, &sx->V[0], NULL);
}
- if(sx->VSUP)
- _unv_print_record(ELEMENT_ID++, SOLIDFEM, v->Num, 4, 6, &sx->V[0], sx->VSUP);
- else
- _unv_print_record(ELEMENT_ID++, SOLIDFEM, v->Num, 4, 0, &sx->V[0], NULL);
+ List_Delete(Elements);
}
- List_Delete(Elements);
// Prisms
Elements = Tree2List(v->Prisms);
@@ -903,6 +916,7 @@ static void _unv_process_groups(Mesh *m)
for(int j = 0; j < List_Nbr(p->Entities); j++) {
List_Read(p->Entities, j, &UNV_VOL_NUM);
Tree_Action(pV->Simplexes, _unv_add_simplex_vertices);
+ Tree_Action(pV->SimplexesBase, _unv_add_simplex_vertices);
Tree_Action(pV->Hexahedra, _unv_add_hexahedron_vertices);
Tree_Action(pV->Prisms, _unv_add_prism_vertices);
Tree_Action(pV->Pyramids, _unv_add_pyramid_vertices);
@@ -974,17 +988,21 @@ static void _gref_consecutive_nodes(Mesh *M, Tree_T *ConsecutiveNTree,
List_T *ListSurfaces = Tree2List(M->Surfaces);
for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
List_Read(ListSurfaces, i, &s);
- List_T *Triangles = Tree2List(s->Simplexes);
- for(int j = 0; j < List_Nbr(Triangles); j++) {
- List_Read(Triangles, j, &sx);
- for(int k = 0; k < 3; k++) {
- if(sx->V[k]->Frozen >= 0) {
- sx->V[k]->Frozen = --newnum;
- Tree_Insert(ConsecutiveNTree, &(sx->V[k]));
- }
+ // triangles
+ for(int simtype = 0; simtype < 2; simtype++){
+ List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
+ for(int j = 0; j < List_Nbr(Triangles); j++) {
+ List_Read(Triangles, j, &sx);
+ for(int k = 0; k < 3; k++) {
+ if(sx->V[k]->Frozen >= 0) {
+ sx->V[k]->Frozen = --newnum;
+ Tree_Insert(ConsecutiveNTree, &(sx->V[k]));
+ }
+ }
}
+ List_Delete(Triangles);
}
- List_Delete(Triangles);
+ // quads
List_T *Quadrangles = Tree2List(s->Quadrangles);
for(int j = 0; j < List_Nbr(Quadrangles); j++) {
List_Read(Quadrangles, j, &qx);
@@ -1002,17 +1020,21 @@ static void _gref_consecutive_nodes(Mesh *M, Tree_T *ConsecutiveNTree,
for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
List_Read(ListSurfaces, i, &s);
- List_T *Triangles = Tree2List(s->Simplexes);
- for(int j = 0; j < List_Nbr(Triangles); j++) {
- List_Read(Triangles, j, &sx);
- for(int k = 0; k < 3; k++) {
- if(sx->VSUP[k]->Frozen >= 0) {
- sx->VSUP[k]->Frozen = --newnum;
- Tree_Insert(ConsecutiveETree, &(sx->VSUP[k]));
+ // triangles
+ for(int simtype = 0; simtype < 2; simtype++){
+ List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
+ for(int j = 0; j < List_Nbr(Triangles); j++) {
+ List_Read(Triangles, j, &sx);
+ for(int k = 0; k < 3; k++) {
+ if(sx->VSUP[k]->Frozen >= 0) {
+ sx->VSUP[k]->Frozen = --newnum;
+ Tree_Insert(ConsecutiveETree, &(sx->VSUP[k]));
+ }
}
}
+ List_Delete(Triangles);
}
- List_Delete(Triangles);
+ // quads
List_T *Quadrangles = Tree2List(s->Quadrangles);
for(int j = 0; j < List_Nbr(Quadrangles); j++) {
List_Read(Quadrangles, j, &qx);
@@ -1042,15 +1064,19 @@ static void _gref_end_consecutive_nodes(Mesh *M)
List_T *ListSurfaces = Tree2List(M->Surfaces);
for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
List_Read(ListSurfaces, i, &s);
- List_T *Triangles = Tree2List(s->Simplexes);
- for(int j = 0; j < List_Nbr(Triangles); j++) {
- List_Read(Triangles, j, &sx);
- for(int k = 0; k < 3; k++)
- sx->V[k]->Frozen = 0;
- for(int k = 0; k < 3; k++)
- sx->VSUP[k]->Frozen = 0;
+ // triangles
+ for(int simtype = 0; simtype < 2; simtype++){
+ List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
+ for(int j = 0; j < List_Nbr(Triangles); j++) {
+ List_Read(Triangles, j, &sx);
+ for(int k = 0; k < 3; k++)
+ sx->V[k]->Frozen = 0;
+ for(int k = 0; k < 3; k++)
+ sx->VSUP[k]->Frozen = 0;
+ }
+ List_Delete(Triangles);
}
- List_Delete(Triangles);
+ // quads
List_T *Quadrangles = Tree2List(s->Quadrangles);
for(int j = 0; j < List_Nbr(Quadrangles); j++) {
List_Read(Quadrangles, j, &qx);
@@ -1087,7 +1113,7 @@ static int _gref_process_nodes(Mesh *M, Tree_T *ConsecutiveNTree,
nbtriqua = 0;
for(i = 0; i < List_Nbr(ListSurfaces); i++) {
List_Read(ListSurfaces, i, &s);
- nbtriqua += Tree_Nbr(s->Simplexes);
+ nbtriqua += Tree_Nbr(s->Simplexes) + Tree_Nbr(s->SimplexesBase);
nbtriqua += Tree_Nbr(s->Quadrangles);
}
List_Delete(ListSurfaces);
@@ -1189,13 +1215,17 @@ static void _gref_process_elements(Mesh *M, int nn)
for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
List_Read(ListSurfaces, i, &s);
- List_T *Triangles = Tree2List(s->Simplexes);
- for(int j = 0; j < List_Nbr(Triangles); j++) {
- List_Read(Triangles, j, &sx);
- fprintf(GREFFILE, "%d %d %d\n", -sx->V[0]->Frozen,
- -sx->V[1]->Frozen, -sx->V[2]->Frozen);
+ // triangles
+ for(int simtype = 0; simtype < 2; simtype++){
+ List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
+ for(int j = 0; j < List_Nbr(Triangles); j++) {
+ List_Read(Triangles, j, &sx);
+ fprintf(GREFFILE, "%d %d %d\n", -sx->V[0]->Frozen,
+ -sx->V[1]->Frozen, -sx->V[2]->Frozen);
+ }
+ List_Delete(Triangles);
}
- List_Delete(Triangles);
+ // quads
List_T *Quadrangles = Tree2List(s->Quadrangles);
for(int j = 0; j < List_Nbr(Quadrangles); j++) {
List_Read(Quadrangles, j, &qx);
@@ -1208,23 +1238,27 @@ static void _gref_process_elements(Mesh *M, int nn)
// Degres de Liberte
for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
List_Read(ListSurfaces, i, &s);
- List_T *Triangles = Tree2List(s->Simplexes);
- for(int j = 0; j < List_Nbr(Triangles); j++) {
- List_Read(Triangles, j, &sx);
- fprintf(GREFFILE, "%d %d %d %d %d %d %d %d %d %d %d %d\n",
- -2 * sx->V[0]->Frozen - 1,
- -2 * sx->V[0]->Frozen,
- -2 * sx->VSUP[0]->Frozen - 1,
- -2 * sx->VSUP[0]->Frozen,
- -2 * sx->V[1]->Frozen - 1,
- -2 * sx->V[1]->Frozen,
- -2 * sx->VSUP[1]->Frozen - 1,
- -2 * sx->VSUP[1]->Frozen,
- -2 * sx->V[2]->Frozen - 1,
- -2 * sx->V[2]->Frozen,
- -2 * sx->VSUP[2]->Frozen - 1, -2 * sx->VSUP[2]->Frozen);
+ // triangles
+ for(int simtype = 0; simtype < 2; simtype++){
+ List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
+ for(int j = 0; j < List_Nbr(Triangles); j++) {
+ List_Read(Triangles, j, &sx);
+ fprintf(GREFFILE, "%d %d %d %d %d %d %d %d %d %d %d %d\n",
+ -2 * sx->V[0]->Frozen - 1,
+ -2 * sx->V[0]->Frozen,
+ -2 * sx->VSUP[0]->Frozen - 1,
+ -2 * sx->VSUP[0]->Frozen,
+ -2 * sx->V[1]->Frozen - 1,
+ -2 * sx->V[1]->Frozen,
+ -2 * sx->VSUP[1]->Frozen - 1,
+ -2 * sx->VSUP[1]->Frozen,
+ -2 * sx->V[2]->Frozen - 1,
+ -2 * sx->V[2]->Frozen,
+ -2 * sx->VSUP[2]->Frozen - 1, -2 * sx->VSUP[2]->Frozen);
+ }
+ List_Delete(Triangles);
}
- List_Delete(Triangles);
+ // quads
List_T *Quadrangles = Tree2List(s->Quadrangles);
for(int j = 0; j < List_Nbr(Quadrangles); j++) {
List_Read(Quadrangles, j, &qx);
@@ -1346,6 +1380,7 @@ static void _wrl_print_all_curves(void *a, void *b)
fprintf(WRLFILE, "DEF Curve%d IndexedLineSet {\n", c->Num);
fprintf(WRLFILE, " coordIndex [\n");
Tree_Action(c->Simplexes, _wrl_print_line);
+ Tree_Action(c->SimplexesBase, _wrl_print_line);
fprintf(WRLFILE, " ]\n");
fprintf(WRLFILE, "}\n");
}
@@ -1356,6 +1391,7 @@ static void _wrl_print_all_surfaces(void *a, void *b)
fprintf(WRLFILE, "DEF Surface%d IndexedFaceSet {\n", s->Num);
fprintf(WRLFILE, " coordIndex [\n");
Tree_Action(s->Simplexes, _wrl_print_triangle);
+ Tree_Action(s->SimplexesBase, _wrl_print_triangle);
Tree_Action(s->Quadrangles, _wrl_print_quadrangle);
fprintf(WRLFILE, " ]\n");
fprintf(WRLFILE, "}\n");
diff --git a/Mesh/Read_Mesh.cpp b/Mesh/Read_Mesh.cpp
index 4f83b8617bb27571910f97746e9a5b37ec1f3691..7f3295066c6ecf837cc015b52764ae18f9396c84 100644
--- a/Mesh/Read_Mesh.cpp
+++ b/Mesh/Read_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Read_Mesh.cpp,v 1.79 2004-11-18 23:42:19 geuzaine Exp $
+// $Id: Read_Mesh.cpp,v 1.80 2004-11-19 18:26:47 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -134,7 +134,7 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
int NbTags, Tag;
double x, y, z, lc1, lc2;
Vertex *vert, verts[NB_NOD_MAX_ELM], *vertsp[NB_NOD_MAX_ELM], **vertspp;
- Simplex *simp;
+ SimplexBase *simp;
Quadrangle *quad;
Hexahedron *hex;
Prism *pri;
@@ -306,7 +306,7 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
case LGN2:
c = addElementaryCurve(M, abs(Elementary));
addPhysicalGroup(M, MSH_PHYSICAL_LINE, Physical, abs(Elementary));
- simp = Create_Simplex_Fast(vertsp[0], vertsp[1], NULL, NULL);
+ simp = Create_SimplexBase(vertsp[0], vertsp[1], NULL, NULL);
simp->Num = Num;
simp->iEnt = Elementary;
simp->iPart = Add_MeshPartition(Partition, M);
@@ -315,16 +315,16 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
simp->VSUP[0] = vertsp[2];
simp->VSUP[0]->Degree = 2;
}
- if(!Tree_Insert(c->Simplexes, &simp)){
+ if(!Tree_Insert(c->SimplexesBase, &simp)){
Msg(GERROR, "Line element %d already exists", simp->Num);
- Free_Simplex(&simp, 0);
+ Free_SimplexBase(&simp, 0);
}
break;
case TRI1:
case TRI2:
s = addElementarySurface(M, Elementary);
addPhysicalGroup(M, MSH_PHYSICAL_SURFACE, Physical, Elementary);
- simp = Create_Simplex_Fast(vertsp[0], vertsp[1], vertsp[2], NULL);
+ simp = Create_SimplexBase(vertsp[0], vertsp[1], vertsp[2], NULL);
simp->Num = Num;
simp->iEnt = Elementary;
simp->iPart = Add_MeshPartition(Partition, M);
@@ -335,9 +335,9 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
simp->VSUP[i]->Degree = 2;
}
}
- if(!Tree_Insert(s->Simplexes, &simp)){
+ if(!Tree_Insert(s->SimplexesBase, &simp)){
Msg(GERROR, "Triangle %d already exists", simp->Num);
- Free_Simplex(&simp, 0);
+ Free_SimplexBase(&simp, 0);
}
break;
case QUA1:
@@ -357,14 +357,14 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
}
if(!Tree_Insert(s->Quadrangles, &quad)){
Msg(GERROR, "Quadrangle %d already exists", simp->Num);
- Free_Simplex(&simp, 0);
+ Free_SimplexBase(&simp, 0);
}
break;
case TET1:
case TET2:
v = addElementaryVolume(M, Elementary);
addPhysicalGroup(M, MSH_PHYSICAL_VOLUME, Physical, Elementary);
- simp = Create_Simplex_Fast(vertsp[0], vertsp[1], vertsp[2], vertsp[3]);
+ simp = Create_SimplexBase(vertsp[0], vertsp[1], vertsp[2], vertsp[3]);
simp->Num = Num;
simp->iEnt = Elementary;
simp->iPart = Add_MeshPartition(Partition, M);
@@ -375,9 +375,9 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
simp->VSUP[i]->Degree = 2;
}
}
- if(!Tree_Insert(v->Simplexes, &simp)){
+ if(!Tree_Insert(v->SimplexesBase, &simp)){
Msg(GERROR, "Tetrahedron %d already exists", simp->Num);
- Free_Simplex(&simp, 0);
+ Free_SimplexBase(&simp, 0);
}
break;
case HEX1:
diff --git a/Mesh/Simplex.cpp b/Mesh/Simplex.cpp
index 8d8c5783421e590e59c689fcf0112afbc5ef58d3..d26d9cdd63ec56529885f4f9ab84323654ac6a62 100644
--- a/Mesh/Simplex.cpp
+++ b/Mesh/Simplex.cpp
@@ -1,4 +1,4 @@
-// $Id: Simplex.cpp,v 1.35 2004-11-18 23:42:19 geuzaine Exp $
+// $Id: Simplex.cpp,v 1.36 2004-11-19 18:26:47 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -30,127 +30,35 @@ extern Context_T CTX;
extern Mesh *THEM, *LOCAL;
extern Simplex MyNewBoundary;
+extern int edges_tetra[6][2];
int FACE_DIMENSION = 2;
-Simplex::Simplex()
-{
- VSUP = NULL;
- V[0] = V[1] = V[2] = V[3] = NULL;
- S[0] = S[1] = S[2] = S[3] = NULL;
- iEnt = -1;
- iPart = -1;
- Quality = 0.;
- Num = ++TotalNumber;
- Visible = VIS_MESH;
-}
+// Basic simplex (contains only pointers to the nodes)
-Simplex::Simplex(Vertex * v1, Vertex * v2, Vertex * v3, Vertex * v4)
+SimplexBase::SimplexBase()
+ : Element()
{
- VSUP = NULL;
- S[0] = S[1] = S[2] = S[3] = NULL;
- Quality = 0.;
- Fourre_Simplexe(v1, v2, v3, v4);
- Num = ++TotalNumber;
- iEnt = -1;
- iPart = -1;
- Visible = VIS_MESH;
-}
-
-Simplex::~Simplex()
-{
- if(VSUP) Free(VSUP);
+ V[0] = V[1] = V[2] = V[3] = NULL;
}
-int Simplex::CircumCircle(double x1, double y1,
- double x2, double y2,
- double x3, double y3, double *xc, double *yc)
+SimplexBase::SimplexBase(Vertex * v1, Vertex * v2, Vertex * v3, Vertex * v4)
+ : Element()
{
- double d, a1, a2, a3;
-
- d = 2. * (double)(y1 * (x2 - x3) + y2 * (x3 - x1) + y3 * (x1 - x2));
- if(d == 0.0) {
- *xc = *yc = -99999.;
- Msg(WARNING, "Degenerated simplex");
- return 0;
- }
-
- a1 = x1 * x1 + y1 * y1;
- a2 = x2 * x2 + y2 * y2;
- a3 = x3 * x3 + y3 * y3;
- *xc = (double)((a1 * (y3 - y2) + a2 * (y1 - y3) + a3 * (y2 - y1)) / d);
- *yc = (double)((a1 * (x2 - x3) + a2 * (x3 - x1) + a3 * (x1 - x2)) / d);
-
- return 1;
+ V[0] = v1; V[1] = v2; V[2] = v3; V[3] = v4;
}
-void Simplex::Center_Circum()
+double SimplexBase::Volume_Simplexe()
{
- /* Calcul du centre de la boule circonscrite */
- int i, N;
- double X[4], Y[4], Z[4];
- double res[3];
+ double mat[3][3];
- if(!V[3])
- N = 3;
+ if(V[3])
+ return (matsimpl(mat) / 6.);
else
- N = 4;
-
- for(i = 0; i < N; i++) {
- X[i] = V[i]->Pos.X;
- Y[i] = V[i]->Pos.Y;
- Z[i] = V[i]->Pos.Z;
- }
-
- if(N == 3) {
- CircumCircle(V[0]->Pos.X, V[0]->Pos.Y,
- V[1]->Pos.X, V[1]->Pos.Y,
- V[2]->Pos.X, V[2]->Pos.Y, &Center.X, &Center.Y);
- Center.Z = 0.0;
- if(fabs(Center.X) > 1.e10)
- Center.X = 1.e10;
- if(fabs(Center.Y) > 1.e10)
- Center.Y = 1.e10;
- Radius = sqrt((X[0] - Center.X) * (X[0] - Center.X) +
- (Y[0] - Center.Y) * (Y[0] - Center.Y));
- }
- else {
- center_tet(X, Y, Z, res);
-
- Center.X = res[0];
- Center.Y = res[1];
- Center.Z = res[2];
- Radius = sqrt((X[0] - Center.X) * (X[0] - Center.X) +
- (Y[0] - Center.Y) * (Y[0] - Center.Y) +
- (Z[0] - Center.Z) * (Z[0] - Center.Z));
- }
-}
-
-int Simplex::Pt_In_Ellipse(Vertex * v, double Metric[3][3])
-{
- double eps, d1, d2, x[2];
-
- Center_Ellipsum_2D(Metric);
-
- x[0] = Center.X - v->Pos.X;
- x[1] = Center.Y - v->Pos.Y;
-
- d1 = Radius;
- d2 = sqrt(x[0] * x[0] * Metric[0][0]
- + x[1] * x[1] * Metric[1][1]
- + 2. * x[0] * x[1] * Metric[0][1]);
-
- eps = fabs(d1 - d2) / (d1 + d2);
- if(eps < 1.e-12) {
- return (1); // Ou Zero ???
- }
- if(d2 < d1)
- return 1;
- return 0;
-
+ return (surfsimpl());
}
-double Simplex::Volume_Simplexe2D()
+double SimplexBase::Volume_Simplexe2D()
{
return ((V[1]->Pos.X - V[0]->Pos.X) *
(V[2]->Pos.Y - V[1]->Pos.Y) -
@@ -158,7 +66,7 @@ double Simplex::Volume_Simplexe2D()
(V[1]->Pos.Y - V[0]->Pos.Y));
}
-void Simplex::center_tet(double X[4], double Y[4], double Z[4], double res[3])
+void SimplexBase::center_tet(double X[4], double Y[4], double Z[4], double res[3])
{
double mat[3][3], b[3], dum;
int i;
@@ -192,7 +100,7 @@ void Simplex::center_tet(double X[4], double Y[4], double Z[4], double res[3])
}
-double Simplex::matsimpl(double mat[3][3])
+double SimplexBase::matsimpl(double mat[3][3])
{
mat[0][0] = V[1]->Pos.X - V[0]->Pos.X;
mat[0][1] = V[2]->Pos.X - V[0]->Pos.X;
@@ -206,14 +114,43 @@ double Simplex::matsimpl(double mat[3][3])
return det3x3(mat);
}
-double Simplex::rhoin()
+double SimplexBase::AireFace(Vertex * V[3])
+{
+ double a[3], b[3], c[3];
+
+ a[0] = V[2]->Pos.X - V[1]->Pos.X;
+ a[1] = V[2]->Pos.Y - V[1]->Pos.Y;
+ a[2] = V[2]->Pos.Z - V[1]->Pos.Z;
+
+ b[0] = V[0]->Pos.X - V[1]->Pos.X;
+ b[1] = V[0]->Pos.Y - V[1]->Pos.Y;
+ b[2] = V[0]->Pos.Z - V[1]->Pos.Z;
+
+ prodve(a, b, c);
+ return (0.5 * sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]));
+}
+
+double SimplexBase::surfsimpl()
+{
+ return AireFace(V);
+}
+
+double SimplexBase::rhoin()
{
double s1, s2, s3, s4;
+
if(V[3]) {
- s1 = fabs(AireFace(F[0].V));
- s2 = fabs(AireFace(F[1].V));
- s3 = fabs(AireFace(F[2].V));
- s4 = fabs(AireFace(F[3].V));
+ Vertex *F[3];
+ F[0] = V[0]; F[1] = V[1]; F[2] = V[2]; s1 = fabs(AireFace(F));
+ F[0] = V[0]; F[1] = V[2]; F[2] = V[3]; s2 = fabs(AireFace(F));
+ if(FACE_DIMENSION == 1) {
+ V[0] = V[1]; V[1] = V[2]; V[2] = V[3]; s3 = fabs(AireFace(F));
+ V[0] = V[0]; V[1] = V[1]; V[2] = V[3]; s4 = fabs(AireFace(F));
+ }
+ else {
+ F[0] = V[0]; F[1] = V[1]; F[2] = V[3]; s3 = fabs(AireFace(F));
+ F[0] = V[1]; F[1] = V[2]; F[2] = V[3]; s4 = fabs(AireFace(F));
+ }
return 3. * fabs(Volume_Simplexe()) / (s1 + s2 + s3 + s4);
}
else {
@@ -221,51 +158,48 @@ double Simplex::rhoin()
}
}
-double Simplex::lij(int i, int j)
+double SimplexBase::maxEdge()
{
- return sqrt(DSQR(V[i]->Pos.X - V[j]->Pos.X) +
- DSQR(V[i]->Pos.Y - V[j]->Pos.Y) +
- DSQR(V[i]->Pos.Z - V[j]->Pos.Z));
-}
+ double maxlij = 0.;
+ int N = V[3] ? 6 : 3;
-double Simplex::Volume_Simplexe()
-{
- double mat[3][3];
+ if(V[3] || V[2])
+ for(int i = 0; i < N; i++)
+ maxlij = DMAX(maxlij, lij(V[edges_tetra[i][0]], V[edges_tetra[i][1]]));
+ else if(V[1])
+ maxlij = lij(V[0], V[1]);
- if(V[3])
- return (matsimpl(mat) / 6.);
- else
- return (surfsimpl());
+ return maxlij;
}
-double Simplex::EtaShapeMeasure()
+double SimplexBase::EtaShapeMeasure()
{
int i, j;
double lij2 = 0.0;
for(i = 0; i <= 3; i++) {
for(j = i + 1; j <= 3; j++) {
- lij2 += DSQR(lij(i, j));
+ lij2 += DSQR(lij(V[i], V[j]));
}
}
return 12. * pow(9. / 10. * DSQR(fabs(Volume_Simplexe())), 1./3.) / (lij2);
}
-double Simplex::RhoShapeMeasure()
+double SimplexBase::RhoShapeMeasure()
{
int i, j;
double minlij = 1.e25, maxlij = 0.0;
for(i = 0; i <= 3; i++) {
for(j = i + 1; j <= 3; j++) {
if(i != j) {
- minlij = DMIN(minlij, fabs(lij(i, j)));
- maxlij = DMAX(maxlij, fabs(lij(i, j)));
+ minlij = DMIN(minlij, fabs(lij(V[i], V[j])));
+ maxlij = DMAX(maxlij, fabs(lij(V[i], V[j])));
}
}
}
return minlij / maxlij;
}
-double Simplex::GammaShapeMeasure()
+double SimplexBase::GammaShapeMeasure()
{
int i, j, N;
double maxlij = 0.0;
@@ -278,12 +212,154 @@ double Simplex::GammaShapeMeasure()
for(i = 0; i <= N - 1; i++) {
for(j = i + 1; j <= N - 1; j++) {
if(i != j)
- maxlij = DMAX(maxlij, lij(i, j));
+ maxlij = DMAX(maxlij, lij(V[i], V[j]));
}
}
return 12. * rhoin() / (sqrt(6.) * maxlij);
}
+void SimplexBase::ExportLcField(FILE * f)
+{
+ if(!V[2])
+ fprintf(f, "SL(%f,%f,%f,%f,%f,%f){%12.5E,%12.5E};\n",
+ V[0]->Pos.X, V[0]->Pos.Y, V[0]->Pos.Z, V[1]->Pos.X, V[1]->Pos.Y,
+ V[1]->Pos.Z, V[0]->lc, V[1]->lc);
+ else if(!V[3])
+ fprintf(f, "ST(%f,%f,%f,%f,%f,%f,%f,%f,%f){%12.5E,%12.5E,%12.5E};\n",
+ V[0]->Pos.X, V[0]->Pos.Y, V[0]->Pos.Z, V[1]->Pos.X, V[1]->Pos.Y,
+ V[1]->Pos.Z, V[2]->Pos.X, V[2]->Pos.Y, V[2]->Pos.Z, V[0]->lc,
+ V[1]->lc, V[2]->lc);
+ else
+ fprintf(f, "SS(%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f){%12.5E,%12.5E,%12.5E,%12.5E};\n",
+ V[0]->Pos.X, V[0]->Pos.Y, V[0]->Pos.Z, V[1]->Pos.X, V[1]->Pos.Y,
+ V[1]->Pos.Z, V[2]->Pos.X, V[2]->Pos.Y, V[2]->Pos.Z, V[3]->Pos.X,
+ V[3]->Pos.Y, V[3]->Pos.Z, V[0]->lc, V[1]->lc, V[2]->lc, V[3]->lc);
+}
+
+SimplexBase *Create_SimplexBase(Vertex * v1, Vertex * v2, Vertex * v3, Vertex * v4)
+{
+ return new SimplexBase(v1, v2, v3, v4);
+}
+
+void Free_SimplexBase(void *a, void *b)
+{
+ SimplexBase *s = *(SimplexBase **) a;
+ if(s) {
+ delete s;
+ s = NULL;
+ }
+}
+
+int compareSimplexBase(const void *a, const void *b)
+{
+ SimplexBase *q = *(SimplexBase **) a;
+ SimplexBase *w = *(SimplexBase **) b;
+ return (q->Num - w->Num);
+}
+
+// Mesh simplex (contains all the data necessary for mesh generation)
+
+Simplex::Simplex()
+ : SimplexBase()
+{
+ S[0] = S[1] = S[2] = S[3] = NULL;
+ Quality = 0.;
+}
+
+Simplex::Simplex(Vertex * v1, Vertex * v2, Vertex * v3, Vertex * v4)
+ : SimplexBase(v1, v2, v3, v4)
+{
+ S[0] = S[1] = S[2] = S[3] = NULL;
+ Quality = 0.;
+ Fourre_Simplexe(v1, v2, v3, v4);
+}
+
+int Simplex::CircumCircle(double x1, double y1, double x2, double y2,
+ double x3, double y3, double *xc, double *yc)
+{
+ double d, a1, a2, a3;
+
+ d = 2. * (double)(y1 * (x2 - x3) + y2 * (x3 - x1) + y3 * (x1 - x2));
+ if(d == 0.0) {
+ *xc = *yc = -99999.;
+ Msg(WARNING, "Degenerated simplex");
+ return 0;
+ }
+
+ a1 = x1 * x1 + y1 * y1;
+ a2 = x2 * x2 + y2 * y2;
+ a3 = x3 * x3 + y3 * y3;
+ *xc = (double)((a1 * (y3 - y2) + a2 * (y1 - y3) + a3 * (y2 - y1)) / d);
+ *yc = (double)((a1 * (x2 - x3) + a2 * (x3 - x1) + a3 * (x1 - x2)) / d);
+
+ return 1;
+}
+
+void Simplex::Center_Circum()
+{
+ /* Calcul du centre de la boule circonscrite */
+ int i, N;
+ double X[4], Y[4], Z[4];
+ double res[3];
+
+ if(!V[3])
+ N = 3;
+ else
+ N = 4;
+
+ for(i = 0; i < N; i++) {
+ X[i] = V[i]->Pos.X;
+ Y[i] = V[i]->Pos.Y;
+ Z[i] = V[i]->Pos.Z;
+ }
+
+ if(N == 3) {
+ CircumCircle(V[0]->Pos.X, V[0]->Pos.Y,
+ V[1]->Pos.X, V[1]->Pos.Y,
+ V[2]->Pos.X, V[2]->Pos.Y, &Center.X, &Center.Y);
+ Center.Z = 0.0;
+ if(fabs(Center.X) > 1.e10)
+ Center.X = 1.e10;
+ if(fabs(Center.Y) > 1.e10)
+ Center.Y = 1.e10;
+ Radius = sqrt((X[0] - Center.X) * (X[0] - Center.X) +
+ (Y[0] - Center.Y) * (Y[0] - Center.Y));
+ }
+ else {
+ center_tet(X, Y, Z, res);
+
+ Center.X = res[0];
+ Center.Y = res[1];
+ Center.Z = res[2];
+ Radius = sqrt((X[0] - Center.X) * (X[0] - Center.X) +
+ (Y[0] - Center.Y) * (Y[0] - Center.Y) +
+ (Z[0] - Center.Z) * (Z[0] - Center.Z));
+ }
+}
+
+int Simplex::Pt_In_Ellipse(Vertex * v, double Metric[3][3])
+{
+ double eps, d1, d2, x[2];
+
+ Center_Ellipsum_2D(Metric);
+
+ x[0] = Center.X - v->Pos.X;
+ x[1] = Center.Y - v->Pos.Y;
+
+ d1 = Radius;
+ d2 = sqrt(x[0] * x[0] * Metric[0][0]
+ + x[1] * x[1] * Metric[1][1]
+ + 2. * x[0] * x[1] * Metric[0][1]);
+
+ eps = fabs(d1 - d2) / (d1 + d2);
+ if(eps < 1.e-12) {
+ return (1); // Ou Zero ???
+ }
+ if(d2 < d1)
+ return 1;
+ return 0;
+
+}
void Simplex::Fourre_Simplexe(Vertex * v1, Vertex * v2, Vertex * v3,
Vertex * v4)
@@ -360,18 +436,6 @@ Simplex *Create_Simplex(Vertex * v1, Vertex * v2, Vertex * v3, Vertex * v4)
return new Simplex(v1, v2, v3, v4);
}
-Simplex *Create_Simplex_Fast(Vertex * v1, Vertex * v2, Vertex * v3, Vertex * v4)
-{
- // bypasses Fourre_Simplex (use for visualization only!)
- Simplex *s = new Simplex();
- s->V[0] = v1;
- s->V[1] = v2;
- s->V[2] = v3;
- s->V[3] = v4;
- s->VSUP = NULL;
- return s;
-}
-
void Free_Simplex(void *a, void *b)
{
Simplex *s = *(Simplex **) a;
@@ -525,11 +589,9 @@ void Simplex::Center_Ellipsum_3D(double m[3][3])
+ (x[2] - z1) * (x[2] - z1) * m[2][2]
+ 2. * (x[0] - x1) * (x[1] - y1) * m[0][1]
+ 2. * (x[0] - x1) * (x[2] - z1) * m[0][2]
- + 2. * (x[1] - y1) * (x[2] - z1) * m[1][2]
- );
+ + 2. * (x[1] - y1) * (x[2] - z1) * m[1][2]);
}
-
int Simplex::Pt_In_Simplex_2D(Vertex * v)
{
double Xmin, Xmax, Ymin, Ymax, Xtr[4], Ytr[4], A[2], B[2], X, Y, Signus[3];
@@ -570,45 +632,6 @@ int Simplex::Pt_In_Simplex_2D(Vertex * v)
return 1;
}
-void Simplex::ExportLcField(FILE * f)
-{
- if(!V[2])
- fprintf(f, "SL(%f,%f,%f,%f,%f,%f){%12.5E,%12.5E};\n",
- V[0]->Pos.X, V[0]->Pos.Y, V[0]->Pos.Z, V[1]->Pos.X, V[1]->Pos.Y,
- V[1]->Pos.Z, V[0]->lc, V[1]->lc);
- else if(!V[3])
- fprintf(f, "ST(%f,%f,%f,%f,%f,%f,%f,%f,%f){%12.5E,%12.5E,%12.5E};\n",
- V[0]->Pos.X, V[0]->Pos.Y, V[0]->Pos.Z, V[1]->Pos.X, V[1]->Pos.Y,
- V[1]->Pos.Z, V[2]->Pos.X, V[2]->Pos.Y, V[2]->Pos.Z, V[0]->lc,
- V[1]->lc, V[2]->lc);
- else
- fprintf(f, "SS(%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f){%12.5E,%12.5E,%12.5E,%12.5E};\n",
- V[0]->Pos.X, V[0]->Pos.Y, V[0]->Pos.Z, V[1]->Pos.X, V[1]->Pos.Y,
- V[1]->Pos.Z, V[2]->Pos.X, V[2]->Pos.Y, V[2]->Pos.Z, V[3]->Pos.X,
- V[3]->Pos.Y, V[3]->Pos.Z, V[0]->lc, V[1]->lc, V[2]->lc, V[3]->lc);
-}
-
-double Simplex::AireFace(Vertex * V[3])
-{
- double a[3], b[3], c[3];
-
- a[0] = V[2]->Pos.X - V[1]->Pos.X;
- a[1] = V[2]->Pos.Y - V[1]->Pos.Y;
- a[2] = V[2]->Pos.Z - V[1]->Pos.Z;
-
- b[0] = V[0]->Pos.X - V[1]->Pos.X;
- b[1] = V[0]->Pos.Y - V[1]->Pos.Y;
- b[2] = V[0]->Pos.Z - V[1]->Pos.Z;
-
- prodve(a, b, c);
- return (0.5 * sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]));
-}
-
-double Simplex::surfsimpl()
-{
- return AireFace(V);
-}
-
bool Simplex::VertexIn(Vertex * v)
{
if(!this || this == &MyNewBoundary)
@@ -746,7 +769,6 @@ bool Simplex::SwapEdge(int iFac)
return true;
}
-
bool Simplex::SwapFace(int iFac, List_T * newsimp, List_T * delsimp)
{
Simplex *s = S[iFac], *s1, *s2, *s3;
diff --git a/Mesh/Simplex.h b/Mesh/Simplex.h
index bc7baa00dcb744147094de2bd515c6ac4e6c2836..9885731be7989b951c79e5a6f689528595041254 100644
--- a/Mesh/Simplex.h
+++ b/Mesh/Simplex.h
@@ -25,51 +25,60 @@
#include "Element.h"
#include "Face.h"
-class Simplex : public Element {
+class SimplexBase : public Element {
public:
- Face F[4]; // 4 faces
Vertex *V[4]; // 4 nodes
+ SimplexBase();
+ ~SimplexBase() {}
+ SimplexBase(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4);
+ double Volume_Simplexe();
+ double Volume_Simplexe2D();
+ double matsimpl(double mat[3][3]);
+ void center_tet(double X[4],double Y[4], double Z[4], double res[3]);
+ double AireFace(Vertex *V[3]);
+ double surfsimpl();
+ double GammaShapeMeasure();
+ double RhoShapeMeasure();
+ double EtaShapeMeasure();
+ double rhoin();
+ double maxEdge();
+ void ExportLcField(FILE *f);
+};
+
+SimplexBase *Create_SimplexBase(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4);
+void Free_SimplexBase(void *a, void *b);
+int compareSimplexBase(const void *a, const void *b);
+
+class Simplex : public SimplexBase {
+ public:
+ Face F[4]; // 4 faces
double Quality; // simplex quality
Coord Center; // CC center
double Radius; // CC radius
Simplex *S[4]; // 4 neighbours
Simplex();
- ~Simplex();
+ ~Simplex() {}
Simplex(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4);
void Fourre_Simplexe(Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4);
- int Pt_In_Simplexe (Vertex *v, double uvw[3], double tol);
+ int Pt_In_Simplexe(Vertex *v, double uvw[3], double tol);
int Pt_In_Simplex_2D(Vertex *v);
void Center_Circum();
- double Volume_Simplexe ();
- double matsimpl(double mat[3][3]);
- void center_tet(double X[4],double Y[4], double Z[4], double res[3]);
- double AireFace (Vertex *V[3]);
- double surfsimpl();
int CircumCircle(double x1,double y1,double x2,double y2,double x3,double y3,
double *xc,double *yc);
- double Volume_Simplexe2D();
- void Center_Ellipsum_2D (double m[3][3]);
- int Pt_In_Ellipse (Vertex *v,double m[3][3]);
+ void Center_Ellipsum_2D(double m[3][3]);
+ int Pt_In_Ellipse(Vertex *v,double m[3][3]);
bool VertexIn(Vertex *v);
bool EdgeIn(Vertex *v1, Vertex *v2, Vertex *v[2]);
- bool SwapEdge (int iFac);
- bool SwapFace (int iFac, List_T *newsimp, List_T *delsimp);
- bool ExtractOppositeEdges ( int iFac, Vertex *p[2], Vertex *q[2]);
- void ExportLcField (FILE *f);
- void Center_Ellipsum_3D (double m[3][3]);
- double GammaShapeMeasure ();
- double RhoShapeMeasure ();
- double EtaShapeMeasure ();
- double lij (int i, int j);
- double rhoin ();
+ bool SwapEdge(int iFac);
+ bool SwapFace(int iFac, List_T *newsimp, List_T *delsimp);
+ bool ExtractOppositeEdges( int iFac, Vertex *p[2], Vertex *q[2]);
+ void Center_Ellipsum_3D(double m[3][3]);
};
int compareSimplex(const void *a, const void *b);
-int compareFace (const void *a, const void *b);
+int compareFace(const void *a, const void *b);
Simplex *Create_Simplex (Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4);
-Simplex *Create_Simplex_Fast (Vertex *v1, Vertex *v2, Vertex *v3, Vertex *v4);
-void Free_Simplex (void *a, void *b);
-
+void Free_Simplex(void *a, void *b);
#endif
diff --git a/doc/VERSIONS b/doc/VERSIONS
index 26547d926e78f4ab6c75e66adad302dc48ababfd..37a481281a135be52090b6008a413a5ae3b9b5ac 100644
--- a/doc/VERSIONS
+++ b/doc/VERSIONS
@@ -1,6 +1,6 @@
-$Id: VERSIONS,v 1.267 2004-11-18 23:51:53 geuzaine Exp $
+$Id: VERSIONS,v 1.268 2004-11-19 18:26:47 geuzaine Exp $
-New since 1.56: generalized displacement maps to display arbitrary
+New in 1.57: generalized displacement maps to display arbitrary
view types; the arrows representing a vector field can now also be
colored by the values from other scalar, vector or tensor fields; new
adaptive high order visualization mode; new options
diff --git a/doc/texinfo/gmsh.texi b/doc/texinfo/gmsh.texi
index 162039c5f457abc6255751561ccfae24dabd5749..569a4ebd62aaac5d04e1bb7840bdf68727147ac5 100644
--- a/doc/texinfo/gmsh.texi
+++ b/doc/texinfo/gmsh.texi
@@ -1,5 +1,5 @@
\input texinfo.tex @c -*-texinfo-*-
-@c $Id: gmsh.texi,v 1.147 2004-11-08 17:05:24 geuzaine Exp $
+@c $Id: gmsh.texi,v 1.148 2004-11-19 18:26:48 geuzaine Exp $
@c
@c Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
@c
@@ -37,8 +37,8 @@
@c =========================================================================
@c %**start of header
@setfilename gmsh.info
-@set EDITION 1.27
-@set GMSH-VERSION 1.56
+@set EDITION 1.28
+@set GMSH-VERSION 1.57
@set GMSH-WEB @uref{http://www.geuz.org/gmsh/}
@set COPYRIGHT @copyright{} 1997-2004 Christophe Geuzaine, Jean-Fran@,{c}ois Remacle
@c
diff --git a/doc/texinfo/opt_mesh.texi b/doc/texinfo/opt_mesh.texi
index f5e7731252ee168a13626c2b1eee64834052d13c..687149fdfe4ea08f9edb624af76ba109645de685 100644
--- a/doc/texinfo/opt_mesh.texi
+++ b/doc/texinfo/opt_mesh.texi
@@ -230,12 +230,12 @@ Default value: @code{0}@*
Saved in: @code{General.OptionsFileName}
@item Mesh.RadiusInf
-Only display elements whose Radius is greater than RadiusInf@*
+Only display elements whose longest edge is greater than RadiusInf@*
Default value: @code{0}@*
Saved in: @code{General.OptionsFileName}
@item Mesh.RadiusSup
-Only display elements whose Radius is smaller than RadiusSup@*
+Only display elements whose longest edge is smaller than RadiusSup@*
Default value: @code{0}@*
Saved in: @code{General.OptionsFileName}