From af0b031a4d497b6ece949f70083c563f49df7447 Mon Sep 17 00:00:00 2001
From: Bastien Gorissen <bastien.gorissen@cenaero.be>
Date: Mon, 27 Jun 2011 09:03:55 +0000
Subject: [PATCH] Added plugin for visualizing discretization error (co-written
 by Thomas Bolemann from the IAG at Stuttgart university and I).

---
 Plugin/CMakeLists.txt          |  14 ++--
 Plugin/DiscretizationError.cpp | 149 +++++++++++++++++++++++++++++++++
 Plugin/DiscretizationError.h   |  32 +++++++
 Plugin/PluginManager.cpp       |  22 ++---
 4 files changed, 195 insertions(+), 22 deletions(-)
 create mode 100644 Plugin/DiscretizationError.cpp
 create mode 100644 Plugin/DiscretizationError.h

diff --git a/Plugin/CMakeLists.txt b/Plugin/CMakeLists.txt
index 38315d9c01..0ffd46c933 100644
--- a/Plugin/CMakeLists.txt
+++ b/Plugin/CMakeLists.txt
@@ -1,4 +1,4 @@
-# Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+# Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
 #
 # See the LICENSE.txt file for license information. Please report all
 # bugs and problems to <gmsh@geuz.org>.
@@ -6,16 +6,16 @@
 set(SRC
   Plugin.cpp PluginManager.cpp
   Levelset.cpp
-  CutPlane.cpp CutSphere.cpp Isosurface.cpp 
+    CutPlane.cpp CutSphere.cpp Isosurface.cpp
   Smooth.cpp CutParametric.cpp
   Lambda2.cpp
   Eigenvectors.cpp Eigenvalues.cpp
-  StreamLines.cpp Particles.cpp CutGrid.cpp CutBox.cpp
+  StreamLines.cpp Particles.cpp CutGrid.cpp
   Transform.cpp
     LongitudeLatitude.cpp
   Triangulate.cpp Tetrahedralize.cpp
   Warp.cpp SphericalRaise.cpp
-  Skin.cpp GSHHS.cpp 
+  Skin.cpp GSHHS.cpp
   MathEval.cpp ModifyComponent.cpp ExtractElements.cpp
   MakeSimplex.cpp
   Integrate.cpp Gradient.cpp Curl.cpp Divergence.cpp
@@ -24,9 +24,9 @@ set(SRC
   HarmonicToTime.cpp ModulusPhase.cpp
   HomologyComputation.cpp
   Distance.cpp ExtractEdges.cpp NearestNeighbor.cpp
-  AnalyseCurvedMesh.cpp FieldFromAmplitudePhase.cpp
-  Bubbles.cpp NearToFarField.cpp
+  AnalyseCurvedMesh.cpp
+  DiscretizationError.cpp
 )
 
-file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h) 
+file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
 append_gmsh_src(Plugin "${SRC};${HDR}")
diff --git a/Plugin/DiscretizationError.cpp b/Plugin/DiscretizationError.cpp
new file mode 100644
index 0000000000..a4cf718281
--- /dev/null
+++ b/Plugin/DiscretizationError.cpp
@@ -0,0 +1,149 @@
+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include "DiscretizationError.h"
+#include "Numeric.h"
+#include <GEntity.h>
+#include <GModel.h>
+#include <Context.h>
+
+// only temp for syntax higlighting
+#include <MQuadrangle.h>
+#include <MTriangle.h>
+
+StringXNumber DiscretizationErrorOptions_Number[] = {
+  {GMSH_FULLRC, "SuperSamplingNodes", NULL, 10.}
+};
+
+extern "C"
+{
+  GMSH_Plugin *GMSH_RegisterDiscretizationErrorPlugin() {
+    return new GMSH_DiscretizationErrorPlugin();
+  }
+}
+
+std::string GMSH_DiscretizationErrorPlugin::getHelp() const
+{
+  return "Plugin(DiscretizationError) computes the error between the mesh and the geometry. It does so by supersampling the elements and computing the distance between the supersampled points dans their projection on the geometry.";
+}
+
+int GMSH_DiscretizationErrorPlugin::getNbOptions() const
+{
+  return sizeof(DiscretizationErrorOptions_Number) / sizeof(StringXNumber);
+}
+
+StringXNumber *GMSH_DiscretizationErrorPlugin::getOption(int iopt)
+{
+  return &DiscretizationErrorOptions_Number[iopt];
+}
+
+PView *GMSH_DiscretizationErrorPlugin::execute(PView *v)
+{
+  double tol = CTX::instance()->geom.tolerance;
+  int nEdgeNodes = (int)DiscretizationErrorOptions_Number[0].def;
+  double paramQuandt = 1.0 / (nEdgeNodes -1) - 10*tol;
+  double paramQuandtQuad = 2.0 / (nEdgeNodes - 1) - 10*tol;
+  int i, j, k, counter;
+  double startEstimate[2] = {0.5, 0.5}; // used as a start estimate for u,v when performing an orthogonal projection
+  double dx,dy,dz,distance;
+
+  std::vector< std::pair<SPoint3,double> > quadDist(nEdgeNodes*nEdgeNodes);
+  std::vector< std::pair<SPoint3,double> > triDist((nEdgeNodes + 1)*nEdgeNodes / 2);
+
+  PView *v2 = new PView();
+  PViewDataList *data2 = getDataList(v2);
+
+
+  for (GModel::fiter itFace = GModel::current()->firstFace(); itFace != GModel::current()->lastFace(); ++itFace) {
+
+    // sample quadrangles
+    /* 13 14 15 16
+     * 9  10 11 12
+     * 5  6  7  8
+     * 1  2  3  4
+     */
+    for (std::vector<MQuadrangle*>::iterator itQuad = (*itFace)->quadrangles.begin(); itQuad != (*itFace)->quadrangles.end(); ++itQuad) {
+      for (j = 0; j < nEdgeNodes; j++) { // u
+        for (i = 0; i < nEdgeNodes; i++) { // v
+          (*itQuad)->pnt(-1+5*tol+paramQuandtQuad*i, -1+5*tol+paramQuandtQuad*j, 0.0, quadDist[j*(nEdgeNodes) + i].first);
+          SPoint3 *point = &quadDist[j*(nEdgeNodes) + i].first;
+          GPoint closest = (*itFace)->closestPoint(*point,startEstimate);
+          dx = closest.x() - point->x();
+          dy = closest.y() - point->y();
+          dz = closest.z() - point->z();
+          quadDist[j*(nEdgeNodes) + i].second = sqrt(dx*dx + dy*dy + dz*dz);
+        }
+      }
+    for (j = 0; j < nEdgeNodes - 1; j++) {
+      for (i = 0; i < nEdgeNodes - 1; i++) {
+	std::vector<double> *out = data2->incrementList(1, TYPE_QUA);
+        int indices[4] = {i+j*nEdgeNodes, i+j*nEdgeNodes+1, i+(j+1)*nEdgeNodes+1, i+(j+1)*nEdgeNodes};
+	for (k = 0; k < 4; k++) out->push_back(quadDist[indices[k]].first.x());
+	for (k = 0; k < 4; k++) out->push_back(quadDist[indices[k]].first.y());
+	for (k = 0; k < 4; k++) out->push_back(quadDist[indices[k]].first.z());
+	for (k = 0; k < 4; k++) out->push_back(quadDist[indices[k]].second);
+      }
+    }
+    }
+
+    // sample triangles
+    /* 10
+     * 6  9
+     * 3  5  8
+     * 1  2  4  7
+     */
+    for (std::vector<MTriangle*>::iterator itTri = (*itFace)->triangles.begin(); itTri != (*itFace)->triangles.end(); ++itTri) {
+      counter = 0;
+      for (i = 0; i < nEdgeNodes; i++) {
+        for (j = 0; j < (i + 1); j++) {
+          (*itTri)->pnt(5*tol + paramQuandt*(i - j), 5*tol + paramQuandt*j, 0.0, triDist[counter].first);
+          SPoint3 *point = &triDist[counter].first; // Check : the points are good
+          GPoint closest = (*itFace)->closestPoint(*point,startEstimate);
+          dx = (closest.x() - point->x());
+          dy = (closest.y() - point->y());
+          dz = (closest.z() - point->z());
+          triDist[counter].second = sqrt(dx*dx + dy*dy + dz*dz);
+          counter++;
+        }
+      }
+
+    int indices[3];
+    for (j = 0; j < nEdgeNodes - 1; j++) { // row in the triangle
+      bool odd = false;
+      for (i = 0; i < j*2 + 1; i ++) { // small tri in the row
+	if (!odd) {
+	  indices[0] = i / 2 + (j+1)*j/2;
+	  indices[1] = i / 2 + (j+1)*j/2 + j+1;
+	  indices[2] = i / 2 + (j+1)*j/2 + j+2;
+	} else {
+	  indices[0] = (i-1)/2 + (j+1)*j/2;
+	  indices[1] = (i-1)/2 + (j+1)*j/2 + j+2;
+	  indices[2] = (i-1)/2 + (j+1)*j/2 + 1;
+	}
+	std::vector<double> *out = data2->incrementList(1, TYPE_TRI);
+	for (k = 0; k < 3; k++) out->push_back(triDist[indices[k]].first.x());
+	for (k = 0; k < 3; k++) out->push_back(triDist[indices[k]].first.y());
+	for (k = 0; k < 3; k++) out->push_back(triDist[indices[k]].first.z());
+	for (k = 0; k < 3; k++) out->push_back(triDist[indices[k]].second);
+	odd = !odd;
+      }
+
+    }
+
+
+    }
+
+    //viusalize stuff
+
+
+
+  }
+
+  data2->setName("Discretization Error");
+  data2->setFileName("discretization_err.pos");
+  data2->finalize();
+
+  return v2;
+}
diff --git a/Plugin/DiscretizationError.h b/Plugin/DiscretizationError.h
new file mode 100644
index 0000000000..52336d2b47
--- /dev/null
+++ b/Plugin/DiscretizationError.h
@@ -0,0 +1,32 @@
+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _DISCRETIZATION_ERROR_H_
+#define _DISCRETIZATION_ERROR_H_
+
+#include "Plugin.h"
+
+extern "C"
+{
+  GMSH_Plugin *GMSH_RegisterDiscretizationErrorPlugin();
+}
+
+class GMSH_DiscretizationErrorPlugin : public GMSH_PostPlugin
+{
+public:
+  GMSH_DiscretizationErrorPlugin(){}
+  std::string getName() const { return "DiscretizationError"; }
+  std::string getShortHelp() const
+  {
+    return "Computes spline distance to the real geometry.";
+  }
+  std::string getHelp() const;
+  std::string getAuthor() const { return "B. Gorissen, T. Bolemann"; }
+  int getNbOptions() const;
+  StringXNumber* getOption(int iopt);
+  PView *execute(PView *);
+};
+
+#endif
diff --git a/Plugin/PluginManager.cpp b/Plugin/PluginManager.cpp
index bfb8adb113..4c6f6ae931 100644
--- a/Plugin/PluginManager.cpp
+++ b/Plugin/PluginManager.cpp
@@ -1,4 +1,4 @@
-// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
@@ -17,7 +17,6 @@
 #include "CutPlane.h"
 #include "CutParametric.h"
 #include "CutSphere.h"
-#include "CutBox.h"
 #include "Skin.h"
 #include "AnalyseCurvedMesh.h"
 #include "MathEval.h"
@@ -49,9 +48,7 @@
 #include "GSHHS.h"
 #include "HomologyComputation.h"
 #include "ExtractEdges.h"
-#include "FieldFromAmplitudePhase.h"
-#include "Bubbles.h"
-#include "NearToFarField.h"
+#include "DiscretizationError.h"
 
 // for testing purposes only :-)
 #undef HAVE_DLOPEN
@@ -92,7 +89,7 @@ GMSH_SolverPlugin *PluginManager::findSolverPlugin()
     GMSH_Plugin *p = it->second;
     if(p->getType() == GMSH_Plugin::GMSH_SOLVER_PLUGIN) {
       return (GMSH_SolverPlugin*)(p);
-    }      
+    }
   }
   return 0;
 }
@@ -108,7 +105,7 @@ void PluginManager::action(std::string pluginName, std::string action, void *dat
     throw "Unknown plugin action";
 }
 
-void PluginManager::setPluginOption(std::string pluginName, std::string option, 
+void PluginManager::setPluginOption(std::string pluginName, std::string option,
                                     std::string value)
 {
   GMSH_Plugin *plugin = find(pluginName);
@@ -124,7 +121,7 @@ void PluginManager::setPluginOption(std::string pluginName, std::string option,
   throw "Unknown plugin option name";
 }
 
-void PluginManager::setPluginOption(std::string pluginName, std::string option, 
+void PluginManager::setPluginOption(std::string pluginName, std::string option,
                                     double value)
 {
   GMSH_Plugin *plugin = find(pluginName);
@@ -159,8 +156,6 @@ void PluginManager::registerDefaultPlugins()
                       ("Particles", GMSH_RegisterParticlesPlugin()));
     allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
                       ("CutGrid", GMSH_RegisterCutGridPlugin()));
-    allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
-                      ("CutBox", GMSH_RegisterCutBoxPlugin()));
     allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
                       ("Isosurface", GMSH_RegisterIsosurfacePlugin()));
     allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
@@ -226,11 +221,8 @@ void PluginManager::registerDefaultPlugins()
     allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
                       ("ExtractEdges", GMSH_RegisterExtractEdgesPlugin()));
     allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
-                      ("FieldFromAmplitudePhase", GMSH_RegisterFieldFromAmplitudePhasePlugin()));
-    allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
-                      ("NearToFarField", GMSH_RegisterNearToFarFieldPlugin()));
-    allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
-                      ("Bubbles", GMSH_RegisterBubblesPlugin()));
+                      ("DiscretizationError", GMSH_RegisterDiscretizationErrorPlugin()));
+
 #if defined(HAVE_TETGEN)
     allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
                       ("Tetrahedralize", GMSH_RegisterTetrahedralizePlugin()));
-- 
GitLab