diff --git a/Plugin/Makefile b/Plugin/Makefile
index 16603a2b8f083df303fcdc5f8a82005da0306848..464be49dfdd8a06936e6297c6d017affcbf8c8aa 100644
--- a/Plugin/Makefile
+++ b/Plugin/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.62 2004-12-06 04:59:09 geuzaine Exp $
+# $Id: Makefile,v 1.63 2004-12-08 03:10:06 geuzaine Exp $
#
# Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
#
@@ -29,7 +29,9 @@ CFLAGS = -DMPICH_SKIP_MPICXX ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = Plugin.cpp\
Levelset.cpp\
CutPlane.cpp CutSphere.cpp CutMap.cpp \
- Smooth.cpp CutParametric.cpp Gradient.cpp Lambda2.cpp\
+ Smooth.cpp CutParametric.cpp\
+ Gradient.cpp Lambda2.cpp\
+ PrincipalStresses.cpp\
Octree.cpp OctreeInternals.cpp OctreePost.cpp\
StreamLines.cpp CutGrid.cpp\
Transform.cpp\
@@ -78,7 +80,7 @@ Plugin.o: Plugin.cpp Plugin.h ../Common/Options.h ../Common/Message.h \
../Geo/Geo.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h ../Common/GmshUI.h \
- Evaluate.h ../Common/Context.h
+ PrincipalStresses.h Evaluate.h ../Common/Context.h
Levelset.o: Levelset.cpp Levelset.h Plugin.h ../Common/Options.h \
../Common/Message.h ../Common/Views.h ../Common/ColorTable.h \
../DataStr/List.h ../Common/VertexArray.h ../Common/SmoothNormals.h \
@@ -127,6 +129,11 @@ Lambda2.o: Lambda2.cpp Plugin.h ../Common/Options.h ../Common/Message.h \
../Common/VertexArray.h ../Common/SmoothNormals.h \
../Common/GmshMatrix.h ../Common/AdaptiveViews.h Lambda2.h \
../Numeric/Numeric.h ../Common/Context.h
+PrincipalStresses.o: PrincipalStresses.cpp Plugin.h ../Common/Options.h \
+ ../Common/Message.h ../Common/Views.h ../Common/ColorTable.h \
+ ../DataStr/List.h ../Common/VertexArray.h ../Common/SmoothNormals.h \
+ ../Common/GmshMatrix.h ../Common/AdaptiveViews.h PrincipalStresses.h \
+ ../Common/Context.h ../DataStr/Malloc.h ../Numeric/Numeric.h
Octree.o: Octree.cpp Octree.h OctreeInternals.h
OctreeInternals.o: OctreeInternals.cpp ../Common/Message.h \
OctreeInternals.h
diff --git a/Plugin/Plugin.cpp b/Plugin/Plugin.cpp
index 5f3a0954fa0e022ca31957b50822bfb744fa5a2b..cf1748cbd6906aee2f10fd395171123e5febe220 100644
--- a/Plugin/Plugin.cpp
+++ b/Plugin/Plugin.cpp
@@ -1,4 +1,4 @@
-// $Id: Plugin.cpp,v 1.64 2004-11-13 22:52:46 geuzaine Exp $
+// $Id: Plugin.cpp,v 1.65 2004-12-08 03:10:06 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -49,6 +49,7 @@
#include "SphericalRaise.h"
#include "DisplacementRaise.h"
#include "StructuralSolver.h"
+#include "PrincipalStresses.h"
#include "Evaluate.h"
#include "Context.h"
@@ -194,6 +195,8 @@ void GMSH_PluginManager::registerDefaultPlugins()
("HarmonicToTime", GMSH_RegisterHarmonicToTimePlugin()));
allPlugins.insert(std::pair < char *, GMSH_Plugin * >
("Integrate", GMSH_RegisterIntegratePlugin()));
+ allPlugins.insert(std::pair < char *, GMSH_Plugin * >
+ ("PrincipalStresses", GMSH_RegisterPrincipalStressesPlugin()));
#if defined(HAVE_TRIANGLE)
allPlugins.insert(std::pair < char *, GMSH_Plugin * >
("Triangulate", GMSH_RegisterTriangulatePlugin()));
diff --git a/Plugin/PrincipalStresses.cpp b/Plugin/PrincipalStresses.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e42f81b1bec041e6fc8ac402dc79ee27badf9648
--- /dev/null
+++ b/Plugin/PrincipalStresses.cpp
@@ -0,0 +1,189 @@
+// $Id: PrincipalStresses.cpp,v 1.1 2004-12-08 03:10:06 geuzaine Exp $
+//
+// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "Plugin.h"
+#include "PrincipalStresses.h"
+#include "List.h"
+#include "Views.h"
+#include "Context.h"
+#include "Malloc.h"
+#include "Numeric.h"
+
+extern Context_T CTX;
+
+StringXNumber PrincipalStressesOptions_Number[] = {
+ {GMSH_FULLRC, "iView", NULL, -1.}
+};
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterPrincipalStressesPlugin()
+ {
+ return new GMSH_PrincipalStressesPlugin();
+ }
+}
+
+GMSH_PrincipalStressesPlugin::GMSH_PrincipalStressesPlugin()
+{
+ ;
+}
+
+void GMSH_PrincipalStressesPlugin::getName(char *name) const
+{
+ strcpy(name, "Principal Stresses");
+}
+
+void GMSH_PrincipalStressesPlugin::getInfos(char *author, char *copyright, char *help_text) const
+{
+ strcpy(author, "C. Geuzaine (geuz@geuz.org)");
+ strcpy(copyright, "DGR (www.multiphysics.com)");
+ strcpy(help_text,
+ "Plugin(PrincipalStresses) computes the min,\n"
+ "mid and max principal stresses from the tensor\n"
+ "view `iView'. If `iView' < 0, the plugin is run on\n"
+ "the current view.\n"
+ "\n"
+ "Plugin(PrincipalStresses) creates three new\n"
+ "views.\n");
+}
+
+int GMSH_PrincipalStressesPlugin::getNbOptions() const
+{
+ return sizeof(PrincipalStressesOptions_Number) / sizeof(StringXNumber);
+}
+
+StringXNumber *GMSH_PrincipalStressesPlugin::getOption(int iopt)
+{
+ return &PrincipalStressesOptions_Number[iopt];
+}
+
+void GMSH_PrincipalStressesPlugin::catchErrorMessage(char *errorMessage) const
+{
+ strcpy(errorMessage, "PrincipalStresses failed...");
+}
+
+static int nonzero(double v[3])
+{
+ for(int i = 0; i < 3; i++)
+ if(fabs(v[i]) > 1.e-15) return 1;
+ return 0;
+}
+
+static void principal_stresses(List_T *inList, int inNb, int nbNod, int nbTime,
+ List_T *minList, int *minNb,
+ List_T *midList, int *midNb,
+ List_T *maxList, int *maxNb)
+{
+ int nbcomplex = 0;
+ int nb = List_Nbr(inList) / inNb;
+ for(int i = 0; i < List_Nbr(inList); i += nb) {
+ for(int j = 0; j < 3 * nbNod; j++){
+ List_Add(minList, List_Pointer_Fast(inList, i + j));
+ List_Add(midList, List_Pointer_Fast(inList, i + j));
+ List_Add(maxList, List_Pointer_Fast(inList, i + j));
+ }
+ for(int j = 0; j < nbTime; j++){
+ for(int k = 0; k < nbNod; k++){
+ double *val = (double *)List_Pointer_Fast(inList, i + 3 * nbNod +
+ nbNod * 9 * j + 9 * k);
+ double A[3][3] = { {val[0], val[1], val[2]},
+ {val[3], val[4], val[5]},
+ {val[6], val[7], val[8]} };
+ double wr[3], wi[3], B[3][3];
+ EigenSolve3x3(A, wr, wi, B);
+ nbcomplex += nonzero(wi);
+ //printf("djf=%g %g %g\n", wr[0], wr[1], wr[2]);
+ //printf("vec1=%g %g %g\n", B[0][0], B[1][0], B[2][0]);
+ //printf("vec1=%g %g %g\n", B[0][1], B[1][1], B[2][1]);
+ //printf("vec1=%g %g %g\n", B[0][2], B[1][2], B[2][2]);
+ for(int l = 0; l < 3; l++){
+ double res;
+ // wrong if there are complex eigenvals (B contains both
+ // real and imag parts: cf. explanation in EigenSolve3x3)
+ res = wr[0] * B[l][0]; List_Add(minList, &res);
+ res = wr[1] * B[l][1]; List_Add(midList, &res);
+ res = wr[2] * B[l][2]; List_Add(maxList, &res);
+ }
+ }
+ }
+ (*minNb)++;
+ (*midNb)++;
+ (*maxNb)++;
+ }
+
+ if(nbcomplex)
+ Msg(GERROR, "%d elements have complex eigenvalues/eigenvectors");
+}
+
+Post_View *GMSH_PrincipalStressesPlugin::execute(Post_View * v)
+{
+ int iView = (int)PrincipalStressesOptions_Number[0].def;
+
+ if(iView < 0)
+ iView = v ? v->Index : 0;
+
+ if(!List_Pointer_Test(CTX.post.list, iView)) {
+ Msg(GERROR, "View[%d] does not exist", iView);
+ return v;
+ }
+
+ Post_View *v1 = *(Post_View **)List_Pointer(CTX.post.list, iView);
+ Post_View *min = BeginView(1);
+ Post_View *mid = BeginView(1);
+ Post_View *max = BeginView(1);
+
+ principal_stresses(v1->TP, v1->NbTP, 1, v1->NbTimeStep,
+ min->VP, &min->NbVP, mid->VP, &mid->NbVP, max->VP, &max->NbVP);
+ principal_stresses(v1->TL, v1->NbTL, 2, v1->NbTimeStep,
+ min->VL, &min->NbVL, mid->VL, &mid->NbVL, max->VL, &max->NbVL);
+ principal_stresses(v1->TT, v1->NbTT, 3, v1->NbTimeStep,
+ min->VT, &min->NbVT, mid->VT, &mid->NbVT, max->VT, &max->NbVT);
+ principal_stresses(v1->TQ, v1->NbTQ, 4, v1->NbTimeStep,
+ min->VQ, &min->NbVQ, mid->VQ, &mid->NbVQ, max->VQ, &max->NbVQ);
+ principal_stresses(v1->TS, v1->NbTS, 4, v1->NbTimeStep,
+ min->VS, &min->NbVS, mid->VS, &mid->NbVS, max->VS, &max->NbVS);
+ principal_stresses(v1->TH, v1->NbTH, 8, v1->NbTimeStep,
+ min->VH, &min->NbVH, mid->VH, &mid->NbVH, max->VH, &max->NbVH);
+ principal_stresses(v1->TI, v1->NbTI, 6, v1->NbTimeStep,
+ min->VI, &min->NbVI, mid->VI, &mid->NbVI, max->VI, &max->NbVI);
+ principal_stresses(v1->TY, v1->NbTY, 5, v1->NbTimeStep,
+ min->VY, &min->NbVY, mid->VY, &mid->NbVY, max->VY, &max->NbVY);
+
+ // copy time data
+ for(int i = 0; i < List_Nbr(v1->Time); i++){
+ List_Add(min->Time, List_Pointer(v1->Time, i));
+ List_Add(mid->Time, List_Pointer(v1->Time, i));
+ List_Add(max->Time, List_Pointer(v1->Time, i));
+ }
+ // finalize
+ char name[1024], filename[1024];
+ sprintf(name, "%s_MinPrincipalStress", v1->Name);
+ sprintf(filename, "%s_MinPrincipalStress.pos", v1->Name);
+ EndView(min, 1, filename, name);
+ sprintf(name, "%s_MidPrincipalStress", v1->Name);
+ sprintf(filename, "%s_MidPrincipalStress.pos", v1->Name);
+ EndView(mid, 1, filename, name);
+ sprintf(name, "%s_MaxPrincipalStress", v1->Name);
+ sprintf(filename, "%s_MaxPrincipalStress.pos", v1->Name);
+ EndView(max, 1, filename, name);
+
+ return NULL;
+}
diff --git a/Plugin/PrincipalStresses.h b/Plugin/PrincipalStresses.h
new file mode 100644
index 0000000000000000000000000000000000000000..814cdf5d5c3623e6a3deb4002594f114dc6596aa
--- /dev/null
+++ b/Plugin/PrincipalStresses.h
@@ -0,0 +1,43 @@
+#ifndef _PRINCIPAL_STRESSES_H_
+#define _PRINCIPAL_STRESSES_H_
+
+// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "Plugin.h"
+#include "List.h"
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterPrincipalStressesPlugin();
+}
+
+class GMSH_PrincipalStressesPlugin : public GMSH_Post_Plugin
+{
+ public:
+ GMSH_PrincipalStressesPlugin();
+ void getName(char *name) const;
+ void getInfos(char *author, char *copyright, char *helpText) const;
+ void catchErrorMessage(char *errorMessage) const;
+ int getNbOptions() const;
+ StringXNumber* getOption(int iopt);
+ Post_View *execute(Post_View *);
+};
+
+#endif
diff --git a/TODO b/TODO
index a4460c1f5a4c3ad2ea4aada2d555bfc478f70751..18591925e8bafb3abe56f25175848f79ebb38fd3 100644
--- a/TODO
+++ b/TODO
@@ -1,4 +1,4 @@
-$Id: TODO,v 1.70 2004-12-06 04:59:08 geuzaine Exp $
+$Id: TODO,v 1.71 2004-12-08 03:10:06 geuzaine Exp $
********************************************************************
@@ -171,59 +171,3 @@ bien toutes les vues, mais tu ne sais acceder qu'a la derniere... Pour
corriger ca, il faudrait que l'on garde comme info supplementaire dans
chaque vue quelle etait sa position dans le fichier.
-********************************************************************
-
-osterwischc@asme.org
-
-Using Gmsh to post-process structural analysis results would be more
-convenient if there were a plugin to calculate maximum, mid, and
-minimum principal stress from a tensor view (as well as the von mises
-which currently exists).
-
-I'm not familiar enough with the code to do it myself but I would think
-it would work somewhat like the "extract" plugin and create 3 new
-vector views for the results labeled basename_maxprin,
-basename_midprin, and basename_minprin.
-
-Thank you for maintaining this excellent software.
--Carl
-
-********************************************************************
-
-Tenseurs, Laurent Stainier <l.stainier@ulg.ac.be>
-
-L'approche que je suggérerais est la suivante. Tout d'abord, on peut
-commencer par le cas particulier (mais sans doute le plus courant) des
-tenseurs symétriques (à valeurs réelles pour être précis). Dans ce cas,
-les trois valeurs propres, et les vecteurs propres associés, sont bien
-définies et réelles. On peut également définir les trois invariants du
-tenseur A: I1=tr(A), I2=0.5*[tr(A)^2-tr(A^2)], I3=det(A), ou
-alternativement tr(A), tr(A^2), tr(A^3). On peut alors envisager
-différentes représentations graphiques:
-- la croix de contrainte (généralisée à 3D), composée à partir des
-valeurs propres et des vecteurs propres;
-- les invariants (sous formes de plusieurs champs scalaires). Le choix
-des invariants n'est pas évident: en mécanique, les plus intéressants
-sont certainement I1, J2, I3. Quid en EM ?
-- composante par composante, ce qui serait d'ailleurs bien utile pour
-les vecteurs également (mais peut-être cela peut-il être fait par
-plugin, dont je n'ai pas encore étudié le mécanisme ??).
-
-Là où ça se complique franchement, c'est pour les tenseurs non
-symétriques. Quoique ! Dans ce cas là, plus aucune garantie sur les
-valeurs propres ... Mais il me semble donc que la solution peut être de
-se ramener au cas précédent via une décomposition polaire A=R.U ou
-A=V.R, où U et V sont des tenseurs symétriques et R un tenseur
-anti-symétrique normé (R^t.R=1 => R^-1=R^t ; U^2=A^t.A, V^2=A.A^t), R
-représentant une rotation rigide en mécanique. Ce dernier a trois
-valeurs propres: 1 (avec un vecteur propre correspondant à l'axe de
-rotation) et exp(+-i \theta) où \theta est l'angle de rotation (à
-vérifier !!). Pour la représentation graphique, il y a l'option de
-laisser tomber R. Sinon, il faudrait de fait voir comment représenter U
-(ou V) avec R.
-
-Remarque finale: le calcul des valeurs et vecteurs propres est
-évidemment un peu lourd au niveau CPU, mais bon, les tenseurs peuvent
-être stockés sous cette forme (les invariants sont facilement calculés
-à partir des v.p.), et cela nous rappelle cette belle expression "No
-free meal" !
diff --git a/configure b/configure
index c2a53868e93c8f41353ef3ebcd06a3fead3ca83d..06c3fdbf0d5a75ff4b98603f2cfe6d832862317f 100755
--- a/configure
+++ b/configure
@@ -3163,7 +3163,7 @@ if test "x$enable_gui" != "xno"; then
GMSH_DIRS="${GMSH_DIRS} Graphics Fltk"
GMSH_LIBS="-Llib -lGmshFltk -lGmshParser -lGmshGraphics -lGmshMesh -lGmshGeo"
- GMSH_LIBS="${GMSH_LIBS} -lGmshNumeric -lGmshCommon -lGmshDataStr -lGmshPlugin"
+ GMSH_LIBS="${GMSH_LIBS} -lGmshCommon -lGmshDataStr -lGmshPlugin -lGmshNumeric"
GMSH_LIBS="${GMSH_LIBS} -lGmshParallel"
FLAGS="-DHAVE_FLTK ${FLAGS}"
diff --git a/configure.in b/configure.in
index 251af8d17418448b9d1798d74f29e21a520a5ff3..1bd9fd5a223f22a9b15d6d9c5e2ffa14277ad297 100644
--- a/configure.in
+++ b/configure.in
@@ -1,4 +1,4 @@
-dnl $Id: configure.in,v 1.56 2004-09-28 17:13:48 geuzaine Exp $
+dnl $Id: configure.in,v 1.57 2004-12-08 03:10:06 geuzaine Exp $
dnl
dnl Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
dnl
@@ -126,7 +126,7 @@ if test "x$enable_gui" != "xno"; then
GMSH_DIRS="${GMSH_DIRS} Graphics Fltk"
GMSH_LIBS="-Llib -lGmshFltk -lGmshParser -lGmshGraphics -lGmshMesh -lGmshGeo"
- GMSH_LIBS="${GMSH_LIBS} -lGmshNumeric -lGmshCommon -lGmshDataStr -lGmshPlugin"
+ GMSH_LIBS="${GMSH_LIBS} -lGmshCommon -lGmshDataStr -lGmshPlugin -lGmshNumeric"
GMSH_LIBS="${GMSH_LIBS} -lGmshParallel"
FLAGS="-DHAVE_FLTK ${FLAGS}"
diff --git a/doc/VERSIONS b/doc/VERSIONS
index 42b1752bf4536e187e572b016a1a9f56f63d5594..598c468c0114b0b4301c812c977fcfbc5e1b0b1d 100644
--- a/doc/VERSIONS
+++ b/doc/VERSIONS
@@ -1,4 +1,4 @@
-$Id: VERSIONS,v 1.271 2004-12-07 04:52:27 geuzaine Exp $
+$Id: VERSIONS,v 1.272 2004-12-08 03:10:06 geuzaine Exp $
New in 1.57: generalized displacement maps to display arbitrary view
types; the arrows representing a vector field can now also be colored
@@ -7,12 +7,13 @@ high order visualization mode; new options (Solver.SocketCommand,
Solver.NameCommand, View.ArrowSizeProportional, View.Normals,
View.Tangents and General.ClipFactor); fixed display of undesired
solver plugin popups; enhanced interactive plugin behaviour; new
-plugins (HarmonicToTime and Integrate); tetrahedral mesh file reading
-speedup (50% faster on large meshes); large memory footprint reduction
-(up to 50%) for the visualization of triangular/tetrahadral meshes;
-the solver interface now supports TCP/IP connections; new generalized
-raise mode (allows to use complex expressions to offset
-post-processing maps); various small bug fixes and enhancements;
+plugins (HarmonicToTime, Integrate, PrincipalStresses); tetrahedral
+mesh file reading speedup (50% faster on large meshes); large memory
+footprint reduction (up to 50%) for the visualization of
+triangular/tetrahadral meshes; the solver interface now supports
+TCP/IP connections; new generalized raise mode (allows to use complex
+expressions to offset post-processing maps); various small bug fixes
+and enhancements;
New in 1.56: new post-processing option to draw a scalar view raised
by a displacement view without using Plugin(DisplacementRaise) (makes
diff --git a/doc/texinfo/opt_plugin.texi b/doc/texinfo/opt_plugin.texi
index 26f4bb3893ba533db6af45e15344732100d3b771..83bf090ba65c95790b19006dd75f23a3d17a923c 100644
--- a/doc/texinfo/opt_plugin.texi
+++ b/doc/texinfo/opt_plugin.texi
@@ -297,6 +297,21 @@ Default value: @code{-1}
Default value: @code{0}
@end table
+@item Plugin(PrincipalStresses)
+Plugin(PrincipalStresses) computes the min,
+mid and max principal stresses from the tensor
+view `iView'. If `iView' < 0, the plugin is run on
+the current view.
+
+Plugin(PrincipalStresses) creates three new
+views.
+
+Numeric options:
+@table @code
+@item iView
+Default value: @code{-1}
+@end table
+
@item Plugin(Skin)
Plugin(Skin) extracts the skin (the boundary) of
the view `iView'. If `iView' < 0, the plugin is run