From e60bf1979b1808fd3a370cacb705bb0ddeb630e3 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Mon, 28 Jun 2004 19:00:22 +0000
Subject: [PATCH] Added hooks to optimize our meshes using netgen (in
 particular: the "special" volume 99999 is now cleaned up after use and all
 the elements/vertices are transfered back into the original volumes)

As a bonus, this finally allowed me to remove the ugly hacks in the volume
mesh drawing routines...
---
 Graphics/Mesh.cpp       | 50 +++++------------------------------------
 Mesh/3D_Mesh_Netgen.cpp | 13 ++++++-----
 Mesh/Create.cpp         | 20 ++++++++++++-----
 Mesh/Create.h           |  1 +
 Mesh/Generator.cpp      | 46 ++++++++++++++++++++++++++++++++-----
 Mesh/Makefile           |  4 ++--
 Mesh/Smoothing.cpp      |  6 ++---
 Netgen/Makefile         | 39 ++++++++++++++++++++++++++++++--
 Netgen/nglib_addon.cpp  | 20 +++++++++++++++++
 Netgen/nglib_addon.h    |  6 +++++
 10 files changed, 135 insertions(+), 70 deletions(-)
 create mode 100644 Netgen/nglib_addon.cpp
 create mode 100644 Netgen/nglib_addon.h

diff --git a/Graphics/Mesh.cpp b/Graphics/Mesh.cpp
index eaa4cd0bb3..277b44aae2 100644
--- a/Graphics/Mesh.cpp
+++ b/Graphics/Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Mesh.cpp,v 1.97 2004-06-23 19:53:52 geuzaine Exp $
+// $Id: Mesh.cpp,v 1.98 2004-06-28 19:00:22 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -225,6 +225,10 @@ void Draw_Mesh(Mesh * M)
 void Draw_Mesh_Volume(void *a, void *b)
 {
   Volume *v = *(Volume **) a;
+  if(!(v->Visible & VIS_MESH))
+    return;
+  theColor = v->Color;
+  thePhysical = getFirstPhysical(MSH_PHYSICAL_VOLUME, v->Num);
   Tree_Action(v->Simplexes, Draw_Mesh_Tetrahedron);
   Tree_Action(v->Hexahedra, Draw_Mesh_Hexahedron);
   Tree_Action(v->Prisms, Draw_Mesh_Prism);
@@ -852,17 +856,6 @@ void Draw_Mesh_Tetrahedron(void *a, void *b)
   if(part && !(*part)->Visible)
     return;
 
-  // FIXME: move this in Draw_Mesh_Volume as soon as a coherent
-  // structure exists for volumes
-  Volume *v = FindVolume(s->iEnt, THEM);
-  if(v){
-    if(!(v->Visible & VIS_MESH))
-      return;
-    theColor = v->Color;
-  }
-  if(CTX.mesh.color_carousel == 2)
-    thePhysical = getFirstPhysical(MSH_PHYSICAL_VOLUME, s->iEnt);
-
   if(CTX.mesh.gamma_sup) {
     tmp = s->GammaShapeMeasure();
     if(tmp < CTX.mesh.gamma_inf || tmp > CTX.mesh.gamma_sup)
@@ -1004,17 +997,6 @@ void Draw_Mesh_Hexahedron(void *a, void *b)
   if(part && !(*part)->Visible)
     return;
 
-  // FIXME: move this in Draw_Mesh_Volume as soon as a coherent
-  // structure exists for volumes
-  Volume *v = FindVolume(h->iEnt, THEM);
-  if(v){
-    if(!(v->Visible & VIS_MESH))
-      return;
-    theColor = v->Color;
-  }
-  if(CTX.mesh.color_carousel == 2)
-    thePhysical = getFirstPhysical(MSH_PHYSICAL_VOLUME, h->iEnt);
-
   int edges = CTX.mesh.volumes_edges;
   int faces = CTX.mesh.volumes_faces;
 
@@ -1157,17 +1139,6 @@ void Draw_Mesh_Prism(void *a, void *b)
   if(part && !(*part)->Visible)
     return;
 
-  // FIXME: move this in Draw_Mesh_Volume as soon as a coherent
-  // structure exists for volumes
-  Volume *v = FindVolume(p->iEnt, THEM);
-  if(v){
-    if(!(v->Visible & VIS_MESH))
-      return;
-    theColor = v->Color;
-  }
-  if(CTX.mesh.color_carousel == 2)
-    thePhysical = getFirstPhysical(MSH_PHYSICAL_VOLUME, p->iEnt);
-
   int edges = CTX.mesh.volumes_edges;
   int faces = CTX.mesh.volumes_faces;
 
@@ -1307,17 +1278,6 @@ void Draw_Mesh_Pyramid(void *a, void *b)
   if(part && !(*part)->Visible)
     return;
 
-  // FIXME: move this in Draw_Mesh_Volume as soon as a coherent
-  // structure exists for volumes
-  Volume *v = FindVolume(p->iEnt, THEM);
-  if(v){
-    if(!(v->Visible & VIS_MESH))
-      return;
-    theColor = v->Color;
-  }
-  if(CTX.mesh.color_carousel == 2)
-    thePhysical = getFirstPhysical(MSH_PHYSICAL_VOLUME, p->iEnt);
-
   int edges = CTX.mesh.volumes_edges;
   int faces = CTX.mesh.volumes_faces;
 
diff --git a/Mesh/3D_Mesh_Netgen.cpp b/Mesh/3D_Mesh_Netgen.cpp
index 85bfd9f8a2..15094430e2 100644
--- a/Mesh/3D_Mesh_Netgen.cpp
+++ b/Mesh/3D_Mesh_Netgen.cpp
@@ -1,4 +1,4 @@
-// $Id: 3D_Mesh_Netgen.cpp,v 1.4 2004-06-28 00:56:07 geuzaine Exp $
+// $Id: 3D_Mesh_Netgen.cpp,v 1.5 2004-06-28 19:00:22 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -47,10 +47,8 @@ void Optimize_Netgen(Volume * v)
 
 #else
 
-extern "C"
-{
 #include "nglib.h"
-}
+#include "nglib_addon.h"
 
 class Netgen{
  private:
@@ -200,11 +198,14 @@ void Netgen::MeshVolume()
 
 void Netgen::OptimizeVolume()
 {
+  Ng_Meshing_Parameters mp;
+  mp.maxh = 1;
+  mp.fineness = 1;
+  mp.secondorder = 0;
   // remove the pure volume vertices from THEM->Vertices and v->Vertices
   // remove the tets from THEM->Simplexes
   // reset v->Simplexes
-  //RemoveIllegalElements(*_ngmesh);
-  //OptimizeVolume(mparam, *_ngmesh, NULL);
+  NgAddOn_OptimizeVolumeMesh(_ngmesh, &mp);
   // transfer vertices and tets back into v and THEM
 }
 
diff --git a/Mesh/Create.cpp b/Mesh/Create.cpp
index dfc5f309af..0edf45c090 100644
--- a/Mesh/Create.cpp
+++ b/Mesh/Create.cpp
@@ -1,4 +1,4 @@
-// $Id: Create.cpp,v 1.60 2004-06-26 05:07:47 geuzaine Exp $
+// $Id: Create.cpp,v 1.61 2004-06-28 19:00:22 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -735,24 +735,32 @@ Volume *Create_Volume(int Num, int Typ)
 }
 
 void Free_Volume(void *a, void *b)
+{
+  Volume *pV = *(Volume **) a;
+  if(pV) {
+    Tree_Action(pV->Simplexes, Free_Simplex);
+    Tree_Action(pV->Hexahedra, Free_Hexahedron);
+    Tree_Action(pV->Prisms, Free_Prism);
+    Tree_Action(pV->Pyramids, Free_Pyramid);
+    Tree_Action(pV->Edges, Free_Edge);
+  }
+  Free_Volume_But_Not_Elements(a, b);
+}
+
+void Free_Volume_But_Not_Elements(void *a, void *b)
 {
   Volume *pV = *(Volume **) a;
   if(pV) {
     List_Delete(pV->TrsfPoints);
     List_Delete(pV->Surfaces);  // surfaces freed elsewhere
     List_Delete(pV->SurfacesOrientations);
-    Tree_Action(pV->Simplexes, Free_Simplex);
     Tree_Delete(pV->Simplexes);
     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
-    Tree_Action(pV->Hexahedra, Free_Hexahedron);
     Tree_Delete(pV->Hexahedra);
-    Tree_Action(pV->Prisms, Free_Prism);
     Tree_Delete(pV->Prisms);
-    Tree_Action(pV->Pyramids, Free_Pyramid);
     Tree_Delete(pV->Pyramids);
-    Tree_Action(pV->Edges, Free_Edge);
     Tree_Delete(pV->Edges);
     Tree_Delete(pV->Faces);
     Free(pV);
diff --git a/Mesh/Create.h b/Mesh/Create.h
index b1c4f61468..05ca74f6c7 100644
--- a/Mesh/Create.h
+++ b/Mesh/Create.h
@@ -52,6 +52,7 @@ void Free_PhysicalGroup(void *a, void *b);
 void Free_MeshPartition(void *a, void *b);
 void Free_Surface(void *a, void *b);
 void Free_Volume(void *a, void *b);
+void Free_Volume_But_Not_Elements(void *a, void *b);
 void Free_Curve(void *a, void *b);
 void Free_EdgeLoop(void *a, void *b);
 void Free_SurfaceLoop(void *a, void *b);
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index b6e512944a..8512c6bb46 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.54 2004-06-20 23:25:32 geuzaine Exp $
+// $Id: Generator.cpp,v 1.55 2004-06-28 19:00:22 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -244,6 +244,24 @@ void Maillage_Dimension_2(Mesh * M)
   M->timing[1] = t2 - t1;
 }
 
+static Volume *IVOL;
+
+void TransferData(void *a, void *b)
+{
+  Simplex *s = *(Simplex**)a;
+  if(s->iEnt == IVOL->Num){
+    if(s->Volume_Simplexe() < 0) {
+      Vertex *temp;
+      temp = s->V[0];
+      s->V[0] = s->V[1];
+      s->V[1] = temp;
+    }
+    Tree_Add(IVOL->Simplexes, &s);
+    for(int i = 0; i < 4; i++)
+      Tree_Insert(IVOL->Vertices, &s->V[i]);
+  }
+}
+
 void Maillage_Dimension_3(Mesh * M)
 {
   Volume *v;
@@ -252,8 +270,8 @@ void Maillage_Dimension_3(Mesh * M)
 
   t1 = Cpu();
 
+  // merge all the delaunay parts in a single special volume
   v = Create_Volume(99999, 99999);
-
   List_T *list = Tree2List(M->Volumes);
   for(int i = 0; i < List_Nbr(list); i++) {
     List_Read(list, i, &vol);
@@ -265,16 +283,32 @@ void Maillage_Dimension_3(Mesh * M)
       }
     }
   }
-  List_Delete(list);
   Tree_Insert(M->Volumes, &v);
 
   if(CTX.mesh.oldxtrude) {
-    Extrude_Mesh_Old(M);        // old automatic extrusion algorithm
+    Extrude_Mesh_Old(M); // old extrusion
   }
   else {
-    Extrude_Mesh(M->Volumes);   // new extrusion
-    Tree_Action(M->Volumes, Maillage_Volume);   // delaunay of remaining parts
+    Extrude_Mesh(M->Volumes); // new extrusion
+    Tree_Action(M->Volumes, Maillage_Volume); // delaunay of remaining parts
+  }
+
+  // transfer data back to individual volumes and remove special volume
+  for(int i = 0; i < List_Nbr(list); i++){
+    List_Read(list, i, &IVOL);
+    Tree_Action(v->Simplexes, TransferData);
   }
+  Tree_Suppress(M->Volumes, &v);
+  Free_Volume_But_Not_Elements(&v, NULL);
+
+  // optimize quality
+  //if(CTX.mesh.quality) {
+  //for(int i = 0; i < List_Nbr(list); i++){
+  //  List_Read(list, i, &v);
+  //  Optimize_Netgen(v);
+  // }
+
+  List_Delete(list);
 
   t2 = Cpu();
 
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 19b464d217..adae4675c7 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.67 2004-06-27 00:20:25 geuzaine Exp $
+# $Id: Makefile,v 1.68 2004-06-28 19:00:22 geuzaine Exp $
 #
 # Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 #
@@ -24,7 +24,7 @@ include ../variables
 LIB     = ../lib/libGmshMesh.a
 INCLUDE = -I../Numeric -I../NR -I../Common -I../DataStr -I../Geo -I../Mesh\
           -I../Graphics -I../Parser -I../Fltk -I../Triangle\
-          -I../Netgen/libsrc/include  -I../Netgen/libsrc/interface
+          -I../Netgen -I../Netgen/libsrc/include  -I../Netgen/libsrc/interface
 CFLAGS  = ${OPTIM} ${FLAGS} ${INCLUDE}
 
 SRC = 1D_Mesh.cpp \
diff --git a/Mesh/Smoothing.cpp b/Mesh/Smoothing.cpp
index f0a329d5a1..fc1361cc5b 100644
--- a/Mesh/Smoothing.cpp
+++ b/Mesh/Smoothing.cpp
@@ -1,4 +1,4 @@
-// $Id: Smoothing.cpp,v 1.13 2004-05-25 04:10:05 geuzaine Exp $
+// $Id: Smoothing.cpp,v 1.14 2004-06-28 19:00:22 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -92,7 +92,7 @@ void ActionLiss(void *data, void *dummy)
   for(i = 0; i < List_Nbr(pnxe->Liste); i++) {
     List_Read(pnxe->Liste, i, &s);
     min_quality_old = DMIN(min_quality_old, s->GammaShapeMeasure());
-    volume_before += s->Volume_Simplexe();
+    volume_before += fabs(s->Volume_Simplexe());
     // On Ne Lisse Point Les Points sur les surfaces quand les volumes
     // sont mailles
     if(s->V[3] && pnxe->v->ListSurf)
@@ -127,7 +127,7 @@ void ActionLiss(void *data, void *dummy)
   double volume_after = 0.0;
   for(i = 0; i < List_Nbr(pnxe->Liste); i++) {
     List_Read(pnxe->Liste, i, &s);
-    volume_after += s->Volume_Simplexe();
+    volume_after += fabs(s->Volume_Simplexe());
   }
   if(fabs(volume_after - volume_before) > 1.e-8 * 
      fabs(volume_after + volume_before) || 
diff --git a/Netgen/Makefile b/Netgen/Makefile
index d178bc3492..29e0c73238 100644
--- a/Netgen/Makefile
+++ b/Netgen/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.1 2004-06-26 17:58:14 geuzaine Exp $
+# $Id: Makefile,v 1.2 2004-06-28 19:00:22 geuzaine Exp $
 #
 # Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 #
@@ -121,7 +121,8 @@ SRC = libsrc/opti/linopt.cpp \
 	libsrc/general/bitarray.cpp \
 	libsrc/general/array.cpp \
 	libsrc/general/symbolta.cpp \
-	libsrc/general/mystring.cpp 
+	libsrc/general/mystring.cpp \
+	nglib_addon.cpp \
 
 OBJ = ${SRC:.cpp=.o}
 
@@ -2924,3 +2925,37 @@ mystring.o: libsrc/general/mystring.cpp libsrc/include/mystdlib.h \
   libsrc/gprim/geomfuncs.hpp libsrc/gprim/geom2d.hpp \
   libsrc/gprim/geom3d.hpp libsrc/gprim/geomtest3d.hpp \
   libsrc/gprim/transform3d.hpp libsrc/gprim/adtree.hpp
+nglib_addon.o: nglib_addon.cpp libsrc/include/meshing.hpp \
+  libsrc/meshing/meshing.hpp libsrc/include/myadt.hpp \
+  libsrc/general/myadt.hpp libsrc/include/mystdlib.h \
+  libsrc/include/mydefs.hpp libsrc/general/ngexception.hpp \
+  libsrc/general/parthreads.hpp libsrc/general/moveablemem.hpp \
+  libsrc/general/dynamicmem.hpp libsrc/general/template.hpp \
+  libsrc/general/array.hpp libsrc/general/table.hpp \
+  libsrc/general/hashtabl.hpp libsrc/general/symbolta.hpp \
+  libsrc/general/bitarray.hpp libsrc/general/flags.hpp \
+  libsrc/general/spbita2d.hpp libsrc/general/seti.hpp \
+  libsrc/general/optmem.hpp libsrc/general/autoptr.hpp \
+  libsrc/general/sort.hpp libsrc/general/stack.hpp \
+  libsrc/general/mystring.hpp libsrc/include/gprim.hpp \
+  libsrc/gprim/gprim.hpp libsrc/gprim/geomobjects.hpp \
+  libsrc/gprim/geomops.hpp libsrc/gprim/geomfuncs.hpp \
+  libsrc/gprim/geom2d.hpp libsrc/gprim/geom3d.hpp \
+  libsrc/gprim/geomtest3d.hpp libsrc/gprim/transform3d.hpp \
+  libsrc/gprim/adtree.hpp libsrc/include/linalg.hpp \
+  libsrc/linalg/linalg.hpp libsrc/linalg/vector.hpp \
+  libsrc/linalg/densemat.hpp libsrc/linalg/polynomial.hpp \
+  libsrc/include/opti.hpp libsrc/opti/opti.hpp \
+  libsrc/meshing/msghandler.hpp libsrc/meshing/meshtype.hpp \
+  libsrc/meshing/localh.hpp libsrc/meshing/meshclass.hpp \
+  libsrc/meshing/global.hpp libsrc/meshing/meshtool.hpp \
+  libsrc/meshing/ruler2.hpp libsrc/meshing/adfront2.hpp \
+  libsrc/meshing/meshing2.hpp libsrc/meshing/improve2.hpp \
+  libsrc/meshing/geomsearch.hpp libsrc/meshing/adfront3.hpp \
+  libsrc/meshing/ruler3.hpp libsrc/meshing/meshing3.hpp \
+  libsrc/meshing/improve3.hpp libsrc/meshing/findip.hpp \
+  libsrc/meshing/topology.hpp libsrc/meshing/curvedelems.hpp \
+  libsrc/meshing/bisect.hpp libsrc/meshing/clusters.hpp \
+  libsrc/meshing/meshfunc.hpp libsrc/meshing/hprefinement.hpp \
+  libsrc/meshing/boundarylayer.hpp libsrc/meshing/specials.hpp \
+  libsrc/interface/nglib.h
diff --git a/Netgen/nglib_addon.cpp b/Netgen/nglib_addon.cpp
new file mode 100644
index 0000000000..9b1c88a817
--- /dev/null
+++ b/Netgen/nglib_addon.cpp
@@ -0,0 +1,20 @@
+// small addition(s) to the netgen interface library
+
+#include <meshing.hpp>
+#include "libsrc/interface/nglib.h"
+
+using namespace netgen;
+
+Ng_Result NgAddOn_OptimizeVolumeMesh (Ng_Mesh * mesh, Ng_Meshing_Parameters * mp)
+{
+  Mesh * m = (Mesh*)mesh;
+
+  MeshingParameters mparam;
+  mparam.maxh = mp->maxh;
+  mparam.meshsizefilename = mp->meshsize_filename;
+
+  RemoveIllegalElements (*m);
+  OptimizeVolume (mparam, *m, NULL);
+
+  return NG_OK;
+}
diff --git a/Netgen/nglib_addon.h b/Netgen/nglib_addon.h
new file mode 100644
index 0000000000..9f1f84f51b
--- /dev/null
+++ b/Netgen/nglib_addon.h
@@ -0,0 +1,6 @@
+#ifndef _NGLIB_ADDON_H_
+#define _NGLIB_ADDON_H_
+
+Ng_Result NgAddOn_OptimizeVolumeMesh (Ng_Mesh * mesh, Ng_Meshing_Parameters * mp);
+
+#endif
-- 
GitLab