Added plugin for visualizing discretization error (co-written by Thomas...

Added plugin for visualizing discretization error (co-written by Thomas Bolemann from the IAG at Stuttgart university and I).

Plugin/CMakeLists.txt

-# 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>.
@@ -10,7 +10,7 @@ set(SRC
   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
@@ -24,8 +24,8 @@ 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)

Plugin/DiscretizationError.cpp

+// 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;
+}

Plugin/DiscretizationError.h

+// 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

Plugin/PluginManager.cpp

-// 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
@@ -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()));