diff --git a/Plugin/CMakeLists.txt b/Plugin/CMakeLists.txt
index a32eda554bf63cccbea7533acf3c428f85f41544..7e4269a4a7fc763e23169f98d33ca2db4d45fe9e 100644
--- a/Plugin/CMakeLists.txt
+++ b/Plugin/CMakeLists.txt
@@ -10,7 +10,7 @@ set(SRC
Smooth.cpp CutParametric.cpp
Lambda2.cpp
Eigenvectors.cpp Eigenvalues.cpp
- StreamLines.cpp CutGrid.cpp
+ StreamLines.cpp Particles.cpp CutGrid.cpp
Transform.cpp
LongitudeLatitude.cpp
Triangulate.cpp Tetrahedralize.cpp
diff --git a/Plugin/Particles.cpp b/Plugin/Particles.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..57623cc7cf8bc6be553dcd273a21371dd11a1e3c
--- /dev/null
+++ b/Plugin/Particles.cpp
@@ -0,0 +1,280 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include <math.h>
+#include "GmshConfig.h"
+#include "Particles.h"
+#include "OctreePost.h"
+#include "Context.h"
+#include "PViewOptions.h"
+
+#if defined(HAVE_OPENGL)
+#include "drawContext.h"
+#endif
+
+StringXNumber ParticlesOptions_Number[] = {
+ {GMSH_FULLRC, "X0", GMSH_ParticlesPlugin::callbackX0, 0.},
+ {GMSH_FULLRC, "Y0", GMSH_ParticlesPlugin::callbackY0, 0.},
+ {GMSH_FULLRC, "Z0", GMSH_ParticlesPlugin::callbackZ0, 0.},
+ {GMSH_FULLRC, "X1", GMSH_ParticlesPlugin::callbackX1, 1.},
+ {GMSH_FULLRC, "Y1", GMSH_ParticlesPlugin::callbackY1, 0.},
+ {GMSH_FULLRC, "Z1", GMSH_ParticlesPlugin::callbackZ1, 0.},
+ {GMSH_FULLRC, "X2", GMSH_ParticlesPlugin::callbackX2, 0.},
+ {GMSH_FULLRC, "Y2", GMSH_ParticlesPlugin::callbackY2, 1.},
+ {GMSH_FULLRC, "Z2", GMSH_ParticlesPlugin::callbackZ2, 0.},
+ {GMSH_FULLRC, "nPointsU", GMSH_ParticlesPlugin::callbackU, 10},
+ {GMSH_FULLRC, "nPointsV", GMSH_ParticlesPlugin::callbackV, 1},
+ {GMSH_FULLRC, "MaxIter", NULL, 100},
+ {GMSH_FULLRC, "DT", NULL, .1},
+ {GMSH_FULLRC, "TimeStep", NULL, 0},
+ {GMSH_FULLRC, "A2", NULL, 1. },
+ {GMSH_FULLRC, "A1", NULL, 0. },
+ {GMSH_FULLRC, "A0", NULL, 0. },
+ {GMSH_FULLRC, "iView", NULL, -1.}
+};
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterParticlesPlugin()
+ {
+ return new GMSH_ParticlesPlugin();
+ }
+}
+
+void GMSH_ParticlesPlugin::draw(void *context)
+{
+#if defined(HAVE_OPENGL)
+ glColor4ubv((GLubyte *) & CTX::instance()->color.fg);
+ drawContext *ctx = (drawContext*)context;
+ double p[3];
+ for(int i = 0; i < getNbU(); ++i){
+ for(int j = 0; j < getNbV(); ++j){
+ getPoint(i, j, p);
+ ctx->drawSphere(CTX::instance()->pointSize, p[0], p[1], p[2], 1);
+ }
+ }
+#endif
+}
+
+double GMSH_ParticlesPlugin::callback(int num, int action, double value, double *opt,
+ double step, double min, double max)
+{
+ switch(action){ // configure the input field
+ case 1: return step;
+ case 2: return min;
+ case 3: return max;
+ default: break;
+ }
+ *opt = value;
+ GMSH_Plugin::setDrawFunction(draw);
+ return 0.;
+}
+
+double GMSH_ParticlesPlugin::callbackX0(int num, int action, double value)
+{
+ return callback(num, action, value, &ParticlesOptions_Number[0].def,
+ CTX::instance()->lc / 100.,
+ -2 * CTX::instance()->lc, 2 * CTX::instance()->lc);
+}
+
+double GMSH_ParticlesPlugin::callbackY0(int num, int action, double value)
+{
+ return callback(num, action, value, &ParticlesOptions_Number[1].def,
+ CTX::instance()->lc / 100.,
+ -2 * CTX::instance()->lc, 2 * CTX::instance()->lc);
+}
+
+double GMSH_ParticlesPlugin::callbackZ0(int num, int action, double value)
+{
+ return callback(num, action, value, &ParticlesOptions_Number[2].def,
+ CTX::instance()->lc / 100.,
+ -2 * CTX::instance()->lc, 2 * CTX::instance()->lc);
+}
+
+double GMSH_ParticlesPlugin::callbackX1(int num, int action, double value)
+{
+ return callback(num, action, value, &ParticlesOptions_Number[3].def,
+ CTX::instance()->lc / 100.,
+ -2 * CTX::instance()->lc, 2 * CTX::instance()->lc);
+}
+
+double GMSH_ParticlesPlugin::callbackY1(int num, int action, double value)
+{
+ return callback(num, action, value, &ParticlesOptions_Number[4].def,
+ CTX::instance()->lc / 100.,
+ -2 * CTX::instance()->lc, 2 * CTX::instance()->lc);
+}
+
+double GMSH_ParticlesPlugin::callbackZ1(int num, int action, double value)
+{
+ return callback(num, action, value, &ParticlesOptions_Number[5].def,
+ CTX::instance()->lc / 100., -2*CTX::instance()->lc, 2*CTX::instance()->lc);
+}
+
+double GMSH_ParticlesPlugin::callbackX2(int num, int action, double value)
+{
+ return callback(num, action, value, &ParticlesOptions_Number[6].def,
+ CTX::instance()->lc / 100., -2*CTX::instance()->lc, 2*CTX::instance()->lc);
+}
+
+double GMSH_ParticlesPlugin::callbackY2(int num, int action, double value)
+{
+ return callback(num, action, value, &ParticlesOptions_Number[7].def,
+ CTX::instance()->lc / 100.,
+ -2 * CTX::instance()->lc, 2 * CTX::instance()->lc);
+}
+
+double GMSH_ParticlesPlugin::callbackZ2(int num, int action, double value)
+{
+ return callback(num, action, value, &ParticlesOptions_Number[8].def,
+ CTX::instance()->lc / 100.,
+ -2 * CTX::instance()->lc, 2 * CTX::instance()->lc);
+}
+
+double GMSH_ParticlesPlugin::callbackU(int num, int action, double value)
+{
+ return callback(num, action, value, &ParticlesOptions_Number[9].def,
+ 1, 1, 100);
+}
+
+double GMSH_ParticlesPlugin::callbackV(int num, int action, double value)
+{
+ return callback(num, action, value, &ParticlesOptions_Number[10].def,
+ 1, 1, 100);
+}
+
+std::string GMSH_ParticlesPlugin::getHelp() const
+{
+ return
+ "Plugin(Particles) computes the trajectory\n"
+ "of a particule in the force field given by the\n"
+ "`TimeStep'-th time step of a vector view\n"
+ "`iView'. The plugin takes as input a grid defined\n"
+ "by the 3 points (`X0',`Y0',`Z0') (origin),\n"
+ "(`X1',`Y1',`Z1') (axis of U) and (`X2',`Y2',`Z2')\n"
+ "(axis of V). The number of points along U and V\n"
+ "that are to be transported is set with the\n"
+ "options `nPointsU' and `nPointsV'. The equation\n"
+ "A2*d^2X(t)/dt^2+A1*dX(t)/dt+A0*X(t)=F(x,y,z) is then\n"
+ "solved with the initial conditions X(t=0) chosen\n"
+ "as the grid, dX/dt(t=0)=0, and with F(x,y,z)\n"
+ "interpolated from the vector view. Time stepping\n"
+ "is done using a Newmark scheme with step size `DT'\n"
+ "and `MaxIter' maximum number of iterations. If\n"
+ "`iView' < 0 the plugin is run on the current view.\n"
+ "\n"
+ "Plugin(Particles) creates one new view. This\n"
+ "view contains multi-step vector points.\n";
+}
+
+int GMSH_ParticlesPlugin::getNbOptions() const
+{
+ return sizeof(ParticlesOptions_Number) / sizeof(StringXNumber);
+}
+
+StringXNumber *GMSH_ParticlesPlugin::getOption(int iopt)
+{
+ return &ParticlesOptions_Number[iopt];
+}
+
+int GMSH_ParticlesPlugin::getNbU()
+{
+ return (int)ParticlesOptions_Number[9].def;
+}
+
+int GMSH_ParticlesPlugin::getNbV()
+{
+ return (int)ParticlesOptions_Number[10].def;
+}
+
+void GMSH_ParticlesPlugin::getPoint(int iU, int iV, double *X)
+{
+ double u = getNbU() > 1 ? (double)iU / (double)(getNbU() - 1.) : 0.;
+ double v = getNbV() > 1 ? (double)iV / (double)(getNbV() - 1.) : 0.;
+ X[0] = ParticlesOptions_Number[0].def +
+ u * (ParticlesOptions_Number[3].def - ParticlesOptions_Number[0].def) +
+ v * (ParticlesOptions_Number[6].def - ParticlesOptions_Number[0].def);
+ X[1] = ParticlesOptions_Number[1].def +
+ u * (ParticlesOptions_Number[4].def - ParticlesOptions_Number[1].def) +
+ v * (ParticlesOptions_Number[7].def - ParticlesOptions_Number[1].def);
+ X[2] = ParticlesOptions_Number[2].def +
+ u * (ParticlesOptions_Number[5].def - ParticlesOptions_Number[2].def) +
+ v * (ParticlesOptions_Number[8].def - ParticlesOptions_Number[2].def);
+}
+
+PView *GMSH_ParticlesPlugin::execute(PView *v)
+{
+ int maxIter = (int)ParticlesOptions_Number[11].def;
+ double DT = ParticlesOptions_Number[12].def;
+ int timeStep = (int)ParticlesOptions_Number[13].def;
+ double A2 = ParticlesOptions_Number[14].def;
+ double A1 = ParticlesOptions_Number[15].def;
+ double A0 = ParticlesOptions_Number[16].def;
+ int iView = (int)ParticlesOptions_Number[17].def;
+
+ PView *v1 = getView(iView, v);
+ if(!v1) return v;
+ PViewData *data1 = v1->getData();
+
+ // sanity checks
+ if(timeStep > data1->getNumTimeSteps() - 1){
+ Msg::Error("Invalid time step (%d) in view[%d]: using 0", v1->getIndex());
+ timeStep = 0;
+ }
+
+ OctreePost o1(v1);
+ double *val2 = 0;
+
+ PView *v3 = new PView();
+ PViewDataList *data3 = getDataList(v3);
+
+ // Solve 'A2 d^2x/dt^2 + A1 dx/dt + A0 x = F' using a Newmark scheme:
+ //
+ // (A2 + gamma DT A1 + beta DT^2 A0) x^{n+1} =
+ // (2 A2 - (1 - 2 gamma) DT A1 - (0.5 + gamma - 2 beta) DT^2 A0) x^n +
+ // (-A2 - (gamma - 1) DT A1 - (0.5 - gamma + beta) DT^2 A0) x^{n-1} +
+ // DT^2 (beta b^{n+1} + (0.5 + gamma - 2 beta) b^n + (0.5 - gamma + beta) b^{n-1})
+ //
+ // coefs for constant acceleration (unconditinonally stable)
+ const double gamma = 0.5, beta = 0.25;
+ double c1 = (A2 + gamma * DT * A1 + beta * DT * DT * A0);
+ double c2 = (2 * A2 - (1 - 2 * gamma) * DT * A1 - (0.5 + gamma - 2 * beta) * DT * DT * A0);
+ double c3 = (-A2 - (gamma - 1) * DT * A1 - (0.5 - gamma + beta) * DT * DT * A0);
+ double c4 = DT * DT * (beta + (0.5 + gamma - 2 * beta) + (0.5 - gamma + beta));
+
+ for(int i = 0; i < getNbU(); ++i){
+ for(int j = 0; j < getNbV(); ++j){
+ double XINIT[3], X0[3], X1[3];
+ getPoint(i, j, XINIT);
+ getPoint(i, j, X0);
+ getPoint(i, j, X1);
+ data3->NbVP++;
+ data3->VP.push_back(XINIT[0]);
+ data3->VP.push_back(XINIT[1]);
+ data3->VP.push_back(XINIT[2]);
+ for(int iter = 0; iter < maxIter; iter++){
+ double F[3], X[3];
+ o1.searchVector(X1[0], X1[1], X1[2], F, timeStep);
+ for(int k = 0; k < 3; k++)
+ X[k] = 1 / c1 * (c2 * X1[k] + c3 * X0[k] + c4 * F[k]);
+ data3->VP.push_back(X[0] - XINIT[0]);
+ data3->VP.push_back(X[1] - XINIT[1]);
+ data3->VP.push_back(X[2] - XINIT[2]);
+ for(int k = 0; k < 3; k++){
+ X0[k] = X1[k];
+ X1[k] = X[k];
+ }
+ }
+ }
+ }
+
+ v3->getOptions()->vectorType = PViewOptions::Displacement;
+
+ data3->setName(data1->getName() + "_Particles");
+ data3->setFileName(data1->getName() + "_Particles.pos");
+ data3->finalize();
+
+ return v3;
+}
diff --git a/Plugin/Particles.h b/Plugin/Particles.h
new file mode 100644
index 0000000000000000000000000000000000000000..b960a17752d9ada78544da2dd397501a38e448ea
--- /dev/null
+++ b/Plugin/Particles.h
@@ -0,0 +1,46 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _PARTICLES_H_
+#define _PARTICLES_H
+
+#include "Plugin.h"
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterParticlesPlugin ();
+}
+
+class GMSH_ParticlesPlugin : public GMSH_PostPlugin
+{
+ static double callback(int num, int action, double value, double *opt,
+ double step, double min, double max);
+ public:
+ GMSH_ParticlesPlugin(){}
+ std::string getName() const { return "Particles"; }
+ std::string getHelp() const;
+ int getNbOptions() const;
+ StringXNumber *getOption(int iopt);
+ PView *execute(PView *);
+
+ static int getNbU();
+ static int getNbV();
+ static void getPoint(int iU, int iV, double *X);
+
+ static double callbackX0(int, int, double);
+ static double callbackY0(int, int, double);
+ static double callbackZ0(int, int, double);
+ static double callbackX1(int, int, double);
+ static double callbackY1(int, int, double);
+ static double callbackZ1(int, int, double);
+ static double callbackX2(int, int, double);
+ static double callbackY2(int, int, double);
+ static double callbackZ2(int, int, double);
+ static double callbackU(int, int, double);
+ static double callbackV(int, int, double);
+ static void draw(void *context);
+};
+
+#endif
diff --git a/Plugin/PluginManager.cpp b/Plugin/PluginManager.cpp
index cb735ab5b7fc246bf7dce2d341da96c72d7806f2..c3c1e1f7cff35582c7e8f53cab1c7dc1a3dc7cfd 100644
--- a/Plugin/PluginManager.cpp
+++ b/Plugin/PluginManager.cpp
@@ -13,6 +13,7 @@
#include "CutMap.h"
#include "CutGrid.h"
#include "StreamLines.h"
+#include "Particles.h"
#include "CutPlane.h"
#include "CutParametric.h"
#include "CutSphere.h"
@@ -147,6 +148,8 @@ void PluginManager::registerDefaultPlugins()
if(CTX::instance()->post.plugins){
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("StreamLines", GMSH_RegisterStreamLinesPlugin()));
+ allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
+ ("Particles", GMSH_RegisterParticlesPlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("CutGrid", GMSH_RegisterCutGridPlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
diff --git a/Plugin/StreamLines.cpp b/Plugin/StreamLines.cpp
index 93244d542c1df5a7261d48e42c13a01df2b12b67..9f556d9da231b5005f2f057c11a66d268caaa54f 100644
--- a/Plugin/StreamLines.cpp
+++ b/Plugin/StreamLines.cpp
@@ -145,10 +145,10 @@ std::string GMSH_StreamLinesPlugin::getHelp() const
"lines. The plugin takes as input a grid defined\n"
"by the 3 points (`X0',`Y0',`Z0') (origin),\n"
"(`X1',`Y1',`Z1') (axis of U) and (`X2',`Y2',`Z2')\n"
- "(axis of V). The number of points that are to\n"
- "be transported along U and V is set with the\n"
+ "(axis of V). The number of points along U and V\n"
+ "that are to be transported is set with the\n"
"options `nPointsU' and `nPointsV'. The equation\n"
- "DX(t)/dt=V(x,y,z) is then solved with the initial\n"
+ "dX(t)/dt=V(x,y,z) is then solved with the initial\n"
"condition X(t=0) chosen as the grid and with V(x,y,z)\n"
"interpolated on the vector view. The time stepping\n"
"scheme is a RK44 with step size `DT' and `MaxIter'\n"
@@ -187,14 +187,14 @@ void GMSH_StreamLinesPlugin::getPoint(int iU, int iV, double *X)
double u = getNbU() > 1 ? (double)iU / (double)(getNbU() - 1.) : 0.;
double v = getNbV() > 1 ? (double)iV / (double)(getNbV() - 1.) : 0.;
X[0] = StreamLinesOptions_Number[0].def +
- u * (StreamLinesOptions_Number[3].def-StreamLinesOptions_Number[0].def) +
- v * (StreamLinesOptions_Number[6].def-StreamLinesOptions_Number[0].def);
+ u * (StreamLinesOptions_Number[3].def - StreamLinesOptions_Number[0].def) +
+ v * (StreamLinesOptions_Number[6].def - StreamLinesOptions_Number[0].def);
X[1] = StreamLinesOptions_Number[1].def +
- u * (StreamLinesOptions_Number[4].def-StreamLinesOptions_Number[1].def) +
- v * (StreamLinesOptions_Number[7].def-StreamLinesOptions_Number[1].def);
+ u * (StreamLinesOptions_Number[4].def - StreamLinesOptions_Number[1].def) +
+ v * (StreamLinesOptions_Number[7].def - StreamLinesOptions_Number[1].def);
X[2] = StreamLinesOptions_Number[2].def +
- u * (StreamLinesOptions_Number[5].def-StreamLinesOptions_Number[2].def) +
- v * (StreamLinesOptions_Number[8].def-StreamLinesOptions_Number[2].def);
+ u * (StreamLinesOptions_Number[5].def - StreamLinesOptions_Number[2].def) +
+ v * (StreamLinesOptions_Number[8].def - StreamLinesOptions_Number[2].def);
}
PView *GMSH_StreamLinesPlugin::execute(PView *v)
@@ -245,7 +245,7 @@ PView *GMSH_StreamLinesPlugin::execute(PView *v)
data3->NbVP++;
data3->VP.push_back(X[0]);
data3->VP.push_back(X[1]);
- data3->VP.push_back(X[2]);
+ data3->VP.push_back(X[2]);
}
int currentTimeStep = 0;
@@ -302,7 +302,7 @@ PView *GMSH_StreamLinesPlugin::execute(PView *v)
else{
data3->VP.push_back(DX[0]);
data3->VP.push_back(DX[1]);
- data3->VP.push_back(DX[2]);
+ data3->VP.push_back(DX[2]);
}
}
}
diff --git a/benchmarks/stl/pelvis.geo b/benchmarks/stl/pelvis.geo
index 7b8aa0c3c1871cb76b205d4820d05d5dba5c19a1..ed4c625bc814c51f6be5ad10d911228acc36c422 100644
--- a/benchmarks/stl/pelvis.geo
+++ b/benchmarks/stl/pelvis.geo
@@ -1,4 +1,4 @@
-Mesh.CharacteristicLengthFactor=0.133;
+//Mesh.CharacteristicLengthFactor=0.133;
Mesh.Algorithm3D = 4; //Frontal (4) Delaunay(1)
Merge "pelvisSMOOTH.stl";