1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5
6 #include "Eigenvectors.h"
7 #include "Numeric.h"
8 #include "fullMatrix.h"
9 #include "GmshDefines.h"
10
11 StringXNumber EigenvectorsOptions_Number[] = {
12 {GMSH_FULLRC, "ScaleByEigenvalues", nullptr, 1.},
13 {GMSH_FULLRC, "View", nullptr, -1.}};
14
15 extern "C" {
GMSH_RegisterEigenvectorsPlugin()16 GMSH_Plugin *GMSH_RegisterEigenvectorsPlugin()
17 {
18 return new GMSH_EigenvectorsPlugin();
19 }
20 }
21
getHelp() const22 std::string GMSH_EigenvectorsPlugin::getHelp() const
23 {
24 return "Plugin(Eigenvectors) computes the three (right) "
25 "eigenvectors of each tensor in the view `View' "
26 "and sorts them according to the value of the "
27 "associated eigenvalues.\n\n"
28 "If `ScaleByEigenvalues' is set, each eigenvector is "
29 "scaled by its associated eigenvalue. The plugin "
30 "gives an error if the eigenvectors are complex.\n\n"
31 "If `View' < 0, the plugin is run on the current view.\n\n"
32 "Plugin(Eigenvectors) creates three new list-based vector view.";
33 }
34
getNbOptions() const35 int GMSH_EigenvectorsPlugin::getNbOptions() const
36 {
37 return sizeof(EigenvectorsOptions_Number) / sizeof(StringXNumber);
38 }
39
getOption(int iopt)40 StringXNumber *GMSH_EigenvectorsPlugin::getOption(int iopt)
41 {
42 return &EigenvectorsOptions_Number[iopt];
43 }
44
execute(PView * v)45 PView *GMSH_EigenvectorsPlugin::execute(PView *v)
46 {
47 int scale = (int)EigenvectorsOptions_Number[0].def;
48 int iView = (int)EigenvectorsOptions_Number[1].def;
49
50 PView *v1 = getView(iView, v);
51 if(!v1) return v;
52
53 PViewData *data1 = getPossiblyAdaptiveData(v1);
54 if(data1->hasMultipleMeshes()) {
55 Msg::Error("Eigenvectors plugin cannot be run on multi-mesh views");
56 return v;
57 }
58
59 PView *min = new PView();
60 PView *mid = new PView();
61 PView *max = new PView();
62
63 PViewDataList *dmin = getDataList(min);
64 PViewDataList *dmid = getDataList(mid);
65 PViewDataList *dmax = getDataList(max);
66
67 int nbcomplex = 0;
68 fullMatrix<double> mat(3, 3), vl(3, 3), vr(3, 3);
69 fullVector<double> dr(3), di(3);
70 for(int ent = 0; ent < data1->getNumEntities(0); ent++) {
71 for(int ele = 0; ele < data1->getNumElements(0, ent); ele++) {
72 if(data1->skipElement(0, ent, ele)) continue;
73 int numComp = data1->getNumComponents(0, ent, ele);
74 if(numComp != 9) continue;
75 int type = data1->getType(0, ent, ele);
76 int numNodes = data1->getNumNodes(0, ent, ele);
77 std::vector<double> *outmin = dmin->incrementList(3, type, numNodes);
78 std::vector<double> *outmid = dmid->incrementList(3, type, numNodes);
79 std::vector<double> *outmax = dmax->incrementList(3, type, numNodes);
80 if(!outmin || !outmid || !outmax) continue;
81 double xyz[3][8];
82 for(int nod = 0; nod < numNodes; nod++)
83 data1->getNode(0, ent, ele, nod, xyz[0][nod], xyz[1][nod], xyz[2][nod]);
84 for(int i = 0; i < 3; i++) {
85 for(int nod = 0; nod < numNodes; nod++) {
86 outmin->push_back(xyz[i][nod]);
87 outmid->push_back(xyz[i][nod]);
88 outmax->push_back(xyz[i][nod]);
89 }
90 }
91 for(int step = 0; step < data1->getNumTimeSteps(); step++) {
92 for(int nod = 0; nod < numNodes; nod++) {
93 for(int i = 0; i < 3; i++)
94 for(int j = 0; j < 3; j++)
95 data1->getValue(step, ent, ele, nod, 3 * i + j, mat(i, j));
96 if(mat.eig(dr, di, vl, vr, true)) {
97 if(!scale) dr(0) = dr(1) = dr(2) = 1.;
98 for(int i = 0; i < 3; i++) {
99 double res;
100 res = dr(0) * vr(i, 0);
101 outmin->push_back(res);
102 res = dr(1) * vr(i, 1);
103 outmid->push_back(res);
104 res = dr(2) * vr(i, 2);
105 outmax->push_back(res);
106 }
107 if(di(0) || di(1) || di(2)) nbcomplex++;
108 }
109 else {
110 Msg::Error("Could not compute eigenvalues/vectors");
111 }
112 }
113 }
114 }
115 }
116
117 if(nbcomplex) Msg::Error("%d tensors have complex eigenvalues", nbcomplex);
118
119 for(int i = 0; i < data1->getNumTimeSteps(); i++) {
120 double time = data1->getTime(i);
121 dmin->Time.push_back(time);
122 dmid->Time.push_back(time);
123 dmax->Time.push_back(time);
124 }
125 dmin->setName(data1->getName() + "_MinEigenvectors");
126 dmin->setFileName(data1->getName() + "_MinEigenvectors.pos");
127 dmin->finalize();
128 dmid->setName(data1->getName() + "_MidEigenvectors");
129 dmid->setFileName(data1->getName() + "_MidEigenvectors.pos");
130 dmid->finalize();
131 dmax->setName(data1->getName() + "_MaxEigenvectors");
132 dmax->setFileName(data1->getName() + "_MaxEigenvectors.pos");
133 dmax->finalize();
134
135 return nullptr;
136 }
137