From 7c3d193bdc37cb3f755938953436be1bf115ec7c Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Wed, 10 Nov 2010 13:25:51 +0000
Subject: [PATCH] Tensor views : eigen values (min/max) and eigen vectors,
 computation of extrema are not (yet) implemented => use 'custom' instead,
 ellipses are coming ...

---
 Common/Options.cpp         |  2 +-
 Fltk/optionWindow.cpp      |  3 +++
 Post/PViewOptions.h        |  5 +++-
 Post/PViewVertexArrays.cpp | 50 +++++++++++++++++++++++++++++++++++---
 4 files changed, 55 insertions(+), 5 deletions(-)

diff --git a/Common/Options.cpp b/Common/Options.cpp
index b413d9370e..de59be40af 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -8011,7 +8011,7 @@ double opt_view_tensor_type(OPT_ARGS_NUM)
   GET_VIEW(0.);
   if(action & GMSH_SET) {
     opt->tensorType = (int)val;
-    if(opt->tensorType != 1)
+    if(opt->tensorType > 4)
       opt->tensorType = 1;
     if(view) view->setChanged(true);
   }
diff --git a/Fltk/optionWindow.cpp b/Fltk/optionWindow.cpp
index 3a079b2dea..aabaf00e19 100644
--- a/Fltk/optionWindow.cpp
+++ b/Fltk/optionWindow.cpp
@@ -3128,6 +3128,9 @@ optionWindow::optionWindow(int deltaFontSize)
       
       static Fl_Menu_Item menu_tensor[] = {
         {"Von-Mises", 0, 0, 0},
+        {"Maximum eigen value", 0, 0, 0},
+        {"Minimum eigen value", 0, 0, 0},
+        {"Eigen vectors", 0, 0, 0},
         {0}
       };
       view.choice[4] = new Fl_Choice
diff --git a/Post/PViewOptions.h b/Post/PViewOptions.h
index febec84dc1..795b300251 100644
--- a/Post/PViewOptions.h
+++ b/Post/PViewOptions.h
@@ -34,7 +34,10 @@ class PViewOptions {
     Displacement = 5
   };
   enum TensorType {
-    VonMises = 1
+    VonMises = 1,
+    MaxEigenValue = 2,
+    MinEigenValue = 3,
+    EigenVectors = 4
   };
   enum GlyphLocation {
     COG = 1,
diff --git a/Post/PViewVertexArrays.cpp b/Post/PViewVertexArrays.cpp
index b5eee4f1b7..139cc9736d 100644
--- a/Post/PViewVertexArrays.cpp
+++ b/Post/PViewVertexArrays.cpp
@@ -944,16 +944,60 @@ static void addVectorElement(PView *p, int ient, int iele, int numNodes,
   }
 }
 
-static void addTensorElement(PView *p, int numNodes, int type,
+static void addTensorElement(PView *p, int iEnt, int iEle, int numNodes, int type,
                              double xyz[PVIEW_NMAX][3], double val[PVIEW_NMAX][9], 
                              bool pre)
 {
   PViewOptions *opt = p->getOptions();
+  fullMatrix <double> tensor(3,3);
+  fullVector<double> S(3), imS (3);
+  fullMatrix<double> V(3,3);
+  fullMatrix <double> rightV(3,3);
 
   if(opt->tensorType == PViewOptions::VonMises){
     for(int i = 0; i < numNodes; i++)
       val[i][0] = ComputeVonMises(val[i]);
     addScalarElement(p, type, xyz, val, pre, numNodes);
+  } else if (opt->tensorType == PViewOptions::MinEigenValue) {
+    for(int i = 0; i < numNodes; i++) {
+      for (int j = 0; j < 3; j++) {
+        tensor(j,0) = val [i][0+j*3];
+        tensor(j,1) = val [i][1+j*3];
+        tensor(j,2) = val [i][2+j*3];
+      }
+      tensor.eig(S, imS, V, rightV, true);
+      val[i][0] = S(0);
+    }
+    addScalarElement(p, type, xyz, val, pre, numNodes);
+  } else if (opt->tensorType == PViewOptions::MaxEigenValue) {
+    for(int i = 0; i < numNodes; i++) {
+      for (int j = 0; j < 3; j++) {
+        tensor(j,0) = val [i][0+j*3];
+        tensor(j,1) = val [i][1+j*3];
+        tensor(j,2) = val [i][2+j*3];
+      }
+      tensor.eig(S, imS, V, rightV, true);
+      val[i][0] = S(2);
+    }
+    addScalarElement(p, type, xyz, val, pre, numNodes);
+  } else if (opt->tensorType == PViewOptions::EigenVectors) {
+    double vval[3][PVIEW_NMAX][9];
+    for(int i = 0; i < numNodes; i++) {
+      for (int j = 0; j < 3; j++) {
+        tensor(j,0) = val [i][0+j*3];
+        tensor(j,1) = val [i][1+j*3];
+        tensor(j,2) = val [i][2+j*3];
+      }
+      tensor.eig(S, imS, V, rightV, false);
+      for (int j = 0; j < 3; j++) {
+        vval[j][i][0] = V(j,0)*S(j);
+        vval[j][i][1] = V(j,1)*S(j);
+        vval[j][i][2] = V(j,2)*S(j);
+      }
+    }
+    addVectorElement(p, iEnt, iEle, numNodes, type, xyz, vval[0], pre);
+    addVectorElement(p, iEnt, iEle, numNodes, type, xyz, vval[1], pre);
+    addVectorElement(p, iEnt, iEle, numNodes, type, xyz, vval[2], pre);
   }
 }
 
@@ -1034,7 +1078,7 @@ static void addElementsInArrays(PView *p, bool preprocessNormalsOnly)
             else if(numComp == 3 && opt->drawVectors)
               addVectorElement(p, ent, i, 1, TYPE_PNT, xyz2, val2, preprocessNormalsOnly);
             else if(numComp == 9 && opt->drawTensors)
-              addTensorElement(p, 1, TYPE_PNT, xyz2, val2, preprocessNormalsOnly);
+              addTensorElement(p, ent, i, 1, TYPE_PNT, xyz2, val2, preprocessNormalsOnly);
           }
         }
         else if(numComp == 1 && opt->drawScalars)
@@ -1042,7 +1086,7 @@ static void addElementsInArrays(PView *p, bool preprocessNormalsOnly)
         else if(numComp == 3 && opt->drawVectors)
           addVectorElement(p, ent, i, numNodes, type, xyz, val, preprocessNormalsOnly);
         else if(numComp == 9 && opt->drawTensors)
-          addTensorElement(p, numNodes, type, xyz, val, preprocessNormalsOnly);
+          addTensorElement(p, ent, i, numNodes, type, xyz, val, preprocessNormalsOnly);
       }
     }
   }
-- 
GitLab