gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
Lambda2.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 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 "Lambda2.h"
7 #include "Numeric.h"
8 
9 StringXNumber Lambda2Options_Number[] = {
10  {GMSH_FULLRC, "Eigenvalue", nullptr, 2.},
11  {GMSH_FULLRC, "View", nullptr, -1.}};
12 
13 extern "C" {
14 GMSH_Plugin *GMSH_RegisterLambda2Plugin() { return new GMSH_Lambda2Plugin(); }
15 }
16 
17 std::string GMSH_Lambda2Plugin::getHelp() const
18 {
19  return "Plugin(Lambda2) computes the eigenvalues "
20  "Lambda(1,2,3) of the tensor (S_ik S_kj + "
21  "Om_ik Om_kj), where S_ij = 0.5 (ui,j + uj,i) "
22  "and Om_ij = 0.5 (ui,j - uj,i) are respectively "
23  "the symmetric and antisymmetric parts of the "
24  "velocity gradient tensor.\n\n"
25  "Vortices are well represented by regions where "
26  "Lambda(2) is negative.\n\n"
27  "If `View' contains tensor elements, the plugin "
28  "directly uses the tensors as the values of the "
29  "velocity gradient tensor; if `View' contains "
30  "vector elements, the plugin uses them as the "
31  "velocities from which to derive the velocity "
32  "gradient tensor.\n\n"
33  "If `View' < 0, the plugin is run on the current view.\n\n"
34  "Plugin(Lambda2) creates one new list-based view.";
35 }
36 
37 int GMSH_Lambda2Plugin::getNbOptions() const
38 {
39  return sizeof(Lambda2Options_Number) / sizeof(StringXNumber);
40 }
41 
42 StringXNumber *GMSH_Lambda2Plugin::getOption(int iopt)
43 {
44  return &Lambda2Options_Number[iopt];
45 }
46 
47 static int inv3x3tran(double mat[3][3], double inv[3][3], double *det)
48 {
49  double ud;
50 
51  *det = det3x3(mat);
52 
53  if(*det == 0.0) return (0);
54 
55  ud = 1. / (*det);
56 
57  inv[0][0] = ud * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]);
58  inv[0][1] = -ud * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]);
59  inv[0][2] = ud * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
60  inv[1][0] = -ud * (mat[0][1] * mat[2][2] - mat[0][2] * mat[2][1]);
61  inv[1][1] = ud * (mat[0][0] * mat[2][2] - mat[0][2] * mat[2][0]);
62  inv[1][2] = -ud * (mat[0][0] * mat[2][1] - mat[0][1] * mat[2][0]);
63  inv[2][0] = ud * (mat[0][1] * mat[1][2] - mat[0][2] * mat[1][1]);
64  inv[2][1] = -ud * (mat[0][0] * mat[1][2] - mat[0][2] * mat[1][0]);
65  inv[2][2] = ud * (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]);
66  return 1;
67 }
68 
69 static void eigen(std::vector<double> &inList, int inNb,
70  std::vector<double> &outList, int *outNb, int nbTime,
71  int nbNod, int nbComp, int lam)
72 {
73  if(!inNb || (nbComp != 3 && nbComp != 9) || lam < 1 || lam > 3) return;
74 
75  // loop on elements
76  int nb = inList.size() / inNb;
77  for(std::size_t i = 0; i < inList.size(); i += nb) {
78  // copy node coordinates
79  for(int j = 0; j < 3 * nbNod; j++) outList.push_back(inList[i + j]);
80 
81  // loop on time steps
82  for(int j = 0; j < nbTime; j++) {
83  double *x = &inList[i];
84  double *y = &inList[i + nbNod];
85  double *z = &inList[i + 2 * nbNod];
86 
87  double GradVel[3][3];
88 
89  if(nbComp == 9) {
90  // val is the velocity gradient tensor: we assume that it is
91  // constant per element
92  double *v = &inList[i + 3 * nbNod + nbNod * nbComp * j + nbComp * 0];
93  GradVel[0][0] = v[0];
94  GradVel[0][1] = v[1];
95  GradVel[0][2] = v[2];
96  GradVel[1][0] = v[3];
97  GradVel[1][1] = v[4];
98  GradVel[1][2] = v[5];
99  GradVel[2][0] = v[6];
100  GradVel[2][1] = v[7];
101  GradVel[2][2] = v[8];
102  }
103  else if(nbComp == 3) {
104  // FIXME: the following could be greatly simplified and
105  // generalized by using the classes in shapeFunctions.h
106 
107  // val contains the velocities: compute the gradient tensor
108  // from them
109  const int MAX_NOD = 4;
110  double val[3][MAX_NOD];
111  for(int k = 0; k < nbNod; k++) {
112  double *v = &inList[i + 3 * nbNod + nbNod * nbComp * j + nbComp * k];
113  for(int l = 0; l < 3; l++) { val[l][k] = v[l]; }
114  }
115  // compute gradient of shape functions
116  double GradPhi_x[MAX_NOD][3];
117  double GradPhi_ksi[MAX_NOD][3];
118  double dx_dksi[3][3];
119  double dksi_dx[3][3];
120  double det;
121  if(nbNod == 3) { // triangles
122  double a[3], b[3], cross[3];
123  a[0] = x[1] - x[0];
124  a[1] = y[1] - y[0];
125  a[2] = z[1] - z[0];
126  b[0] = x[2] - x[0];
127  b[1] = y[2] - y[0];
128  b[2] = z[2] - z[0];
129  prodve(a, b, cross);
130  dx_dksi[0][0] = x[1] - x[0];
131  dx_dksi[0][1] = x[2] - x[0];
132  dx_dksi[0][2] = cross[0];
133  dx_dksi[1][0] = y[1] - y[0];
134  dx_dksi[1][1] = y[2] - y[0];
135  dx_dksi[1][2] = cross[1];
136  dx_dksi[2][0] = z[1] - z[0];
137  dx_dksi[2][1] = z[2] - z[0];
138  dx_dksi[2][2] = cross[2];
139  inv3x3tran(dx_dksi, dksi_dx, &det);
140  GradPhi_ksi[0][0] = -1;
141  GradPhi_ksi[0][1] = -1;
142  GradPhi_ksi[0][2] = 0;
143  GradPhi_ksi[1][0] = 1;
144  GradPhi_ksi[1][1] = 0;
145  GradPhi_ksi[1][2] = 0;
146  GradPhi_ksi[2][0] = 0;
147  GradPhi_ksi[2][1] = 1;
148  GradPhi_ksi[2][2] = 0;
149  }
150  else if(nbNod == 4) { // tetrahedra
151  dx_dksi[0][0] = x[1] - x[0];
152  dx_dksi[0][1] = x[2] - x[0];
153  dx_dksi[0][2] = x[3] - x[0];
154  dx_dksi[1][0] = y[1] - y[0];
155  dx_dksi[1][1] = y[2] - y[0];
156  dx_dksi[1][2] = y[3] - y[0];
157  dx_dksi[2][0] = z[1] - z[0];
158  dx_dksi[2][1] = z[2] - z[0];
159  dx_dksi[2][2] = z[3] - z[0];
160  inv3x3tran(dx_dksi, dksi_dx, &det);
161  GradPhi_ksi[0][0] = -1;
162  GradPhi_ksi[0][1] = -1;
163  GradPhi_ksi[0][2] = -1;
164  GradPhi_ksi[1][0] = 1;
165  GradPhi_ksi[1][1] = 0;
166  GradPhi_ksi[1][2] = 0;
167  GradPhi_ksi[2][0] = 0;
168  GradPhi_ksi[2][1] = 1;
169  GradPhi_ksi[2][2] = 0;
170  GradPhi_ksi[3][0] = 0;
171  GradPhi_ksi[3][1] = 0;
172  GradPhi_ksi[3][2] = 1;
173  }
174  else {
175  Msg::Error("Lambda2 not ready for this type of element");
176  return;
177  }
178  for(int k = 0; k < nbNod; k++) {
179  for(int l = 0; l < 3; l++) {
180  GradPhi_x[k][l] = 0.0;
181  for(int m = 0; m < 3; m++) {
182  GradPhi_x[k][l] += GradPhi_ksi[k][m] * dksi_dx[l][m];
183  }
184  }
185  }
186  // compute gradient of velocities
187  for(int k = 0; k < 3; k++) {
188  for(int l = 0; l < 3; l++) {
189  GradVel[k][l] = 0.0;
190  for(int m = 0; m < nbNod; m++) {
191  GradVel[k][l] += val[k][m] * GradPhi_x[m][l];
192  }
193  }
194  }
195  }
196  else
197  for(int k = 0; k < 3; k++)
198  for(int l = 0; l < 3; l++) GradVel[k][l] = 0.0;
199 
200  // compute the sym and antisymetric parts
201  double sym[3][3];
202  double asym[3][3];
203  for(int m = 0; m < 3; m++) {
204  for(int n = 0; n < 3; n++) {
205  sym[m][n] = 0.5 * (GradVel[m][n] + GradVel[n][m]);
206  asym[m][n] = 0.5 * (GradVel[m][n] - GradVel[n][m]);
207  }
208  }
209  double a[3][3];
210  for(int m = 0; m < 3; m++) {
211  for(int n = 0; n < 3; n++) {
212  a[m][n] = 0.0;
213  for(int l = 0; l < 3; l++)
214  a[m][n] += sym[m][l] * sym[l][n] + asym[m][l] * asym[l][n];
215  }
216  }
217 
218  // compute the eigenvalues
219  double lambda[3];
220  eigenvalue(a, lambda);
221  for(int k = 0; k < nbNod; k++) outList.push_back(lambda[lam - 1]);
222  }
223 
224  (*outNb)++;
225  }
226 }
227 
228 PView *GMSH_Lambda2Plugin::execute(PView *v)
229 {
230  int ev = (int)Lambda2Options_Number[0].def;
231  int iView = (int)Lambda2Options_Number[1].def;
232 
233  PView *v1 = getView(iView, v);
234  if(!v1) return v;
235 
236  PViewDataList *data1 = getDataList(v1);
237  if(!data1) return v;
238 
239  PView *v2 = new PView();
240 
241  PViewDataList *data2 = getDataList(v2);
242  if(!data2) return v;
243 
244  // assume that the tensors contain the velocity gradient tensor
245  int nts = data1->getNumTimeSteps();
246  eigen(data1->TP, data1->NbTP, data2->SP, &data2->NbSP, nts, 1, 9, ev);
247  eigen(data1->TL, data1->NbTL, data2->SL, &data2->NbSL, nts, 2, 9, ev);
248  eigen(data1->TT, data1->NbTT, data2->ST, &data2->NbST, nts, 3, 9, ev);
249  eigen(data1->TQ, data1->NbTQ, data2->SQ, &data2->NbSQ, nts, 4, 9, ev);
250  eigen(data1->TS, data1->NbTS, data2->SS, &data2->NbSS, nts, 4, 9, ev);
251  eigen(data1->TH, data1->NbTH, data2->SH, &data2->NbSH, nts, 8, 9, ev);
252  eigen(data1->TI, data1->NbTI, data2->SI, &data2->NbSI, nts, 6, 9, ev);
253  eigen(data1->TY, data1->NbTY, data2->SY, &data2->NbSY, nts, 5, 9, ev);
254 
255  // assume that the vectors contain the velocities
256  // FIXME: only implemented for tri/tet at the moment
257  // eigen(data1->VP, data1->NbVP, data2->SP, &data2->NbSP, nts, 1, 3, ev);
258  // eigen(data1->VL, data1->NbVL, data2->SL, &data2->NbSL, nts, 2, 3, ev);
259  eigen(data1->VT, data1->NbVT, data2->ST, &data2->NbST, nts, 3, 3, ev);
260  // eigen(data1->VQ, data1->NbVQ, data2->SQ, &data2->NbSQ, nts, 4, 3, ev);
261  eigen(data1->VS, data1->NbVS, data2->SS, &data2->NbSS, nts, 4, 3, ev);
262  // eigen(data1->VH, data1->NbVH, data2->SH, &data2->NbSH, nts, 8, 3, ev);
263  // eigen(data1->VI, data1->NbVI, data2->SI, &data2->NbSI, nts, 6, 3, ev);
264  // eigen(data1->VY, data1->NbVY, data2->SY, &data2->NbSY, nts, 5, 3, ev);
265 
266  data2->Time = data1->Time;
267  data2->setName(data1->getName() + "_Lambda2");
268  data2->setFileName(data1->getName() + "_Lambda2.pos");
269  data2->finalize();
270 
271  return v2;
272 }
GMSH_Lambda2Plugin::getOption
StringXNumber * getOption(int iopt)
Definition: Lambda2.cpp:42
PViewDataList::TS
std::vector< double > TS
Definition: PViewDataList.h:37
PView
Definition: PView.h:27
PViewDataList::TQ
std::vector< double > TQ
Definition: PViewDataList.h:33
GMSH_RegisterLambda2Plugin
GMSH_Plugin * GMSH_RegisterLambda2Plugin()
Definition: Lambda2.cpp:14
PViewDataList::NbTL
int NbTL
Definition: PViewDataList.h:28
PViewDataList::SS
std::vector< double > SS
Definition: PViewDataList.h:37
PViewDataList::TL
std::vector< double > TL
Definition: PViewDataList.h:29
GMSH_Plugin
Definition: Plugin.h:26
PViewDataList::NbTP
int NbTP
Definition: PViewDataList.h:26
PViewDataList::NbTT
int NbTT
Definition: PViewDataList.h:30
PViewDataList
Definition: PViewDataList.h:17
GMSH_Lambda2Plugin::getNbOptions
int getNbOptions() const
Definition: Lambda2.cpp:37
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
PViewDataList::NbSP
int NbSP
Definition: PViewDataList.h:26
det3x3
double det3x3(double mat[3][3])
Definition: Numeric.cpp:126
StringXNumber
Definition: Options.h:918
PViewDataList::TI
std::vector< double > TI
Definition: PViewDataList.h:41
PViewDataList::NbSQ
int NbSQ
Definition: PViewDataList.h:32
PViewDataList::SY
std::vector< double > SY
Definition: PViewDataList.h:43
PViewDataList::NbSH
int NbSH
Definition: PViewDataList.h:38
inv3x3tran
static int inv3x3tran(double mat[3][3], double inv[3][3], double *det)
Definition: Lambda2.cpp:47
cross
void cross(const double *pt0, const double *pt1, const double *pt2, double *cross)
Definition: gmshLevelset.cpp:210
PViewData::setFileName
virtual void setFileName(const std::string &val)
Definition: PViewData.h:75
Lambda2.h
GMSH_Lambda2Plugin::getHelp
std::string getHelp() const
Definition: Lambda2.cpp:17
PViewDataList::TY
std::vector< double > TY
Definition: PViewDataList.h:43
PViewDataList::TP
std::vector< double > TP
Definition: PViewDataList.h:27
PViewDataList::getNumTimeSteps
int getNumTimeSteps()
Definition: PViewDataList.h:79
PViewDataList::NbTH
int NbTH
Definition: PViewDataList.h:38
PViewDataList::SP
std::vector< double > SP
Definition: PViewDataList.h:27
GMSH_FULLRC
#define GMSH_FULLRC
Definition: Options.h:20
Numeric.h
PViewDataList::SL
std::vector< double > SL
Definition: PViewDataList.h:29
PViewDataList::NbVS
int NbVS
Definition: PViewDataList.h:36
PViewDataList::finalize
bool finalize(bool computeMinMax=true, const std::string &interpolationScheme="")
Definition: PViewDataList.cpp:81
PViewData::setName
virtual void setName(const std::string &val)
Definition: PViewData.h:71
Lambda2Options_Number
StringXNumber Lambda2Options_Number[]
Definition: Lambda2.cpp:9
PViewDataList::ST
std::vector< double > ST
Definition: PViewDataList.h:31
prodve
void prodve(double a[3], double b[3], double c[3])
Definition: Numeric.h:105
PViewDataList::TT
std::vector< double > TT
Definition: PViewDataList.h:31
PViewDataList::NbTS
int NbTS
Definition: PViewDataList.h:36
PViewDataList::NbSL
int NbSL
Definition: PViewDataList.h:28
PViewDataList::TH
std::vector< double > TH
Definition: PViewDataList.h:39
PViewDataList::SH
std::vector< double > SH
Definition: PViewDataList.h:39
PViewDataList::NbTI
int NbTI
Definition: PViewDataList.h:40
PViewDataList::NbSS
int NbSS
Definition: PViewDataList.h:36
PViewDataList::NbTY
int NbTY
Definition: PViewDataList.h:42
eigen
static void eigen(std::vector< double > &inList, int inNb, std::vector< double > &outList, int *outNb, int nbTime, int nbNod, int nbComp, int lam)
Definition: Lambda2.cpp:69
z
const double z
Definition: GaussQuadratureQuad.cpp:56
GMSH_Lambda2Plugin::execute
PView * execute(PView *)
Definition: Lambda2.cpp:228
PViewDataList::NbSI
int NbSI
Definition: PViewDataList.h:40
GMSH_PostPlugin::getView
virtual PView * getView(int index, PView *view)
Definition: Plugin.cpp:81
PViewDataList::Time
std::vector< double > Time
Definition: PViewDataList.h:25
PViewDataList::VS
std::vector< double > VS
Definition: PViewDataList.h:37
eigenvalue
void eigenvalue(double mat[3][3], double v[3])
Definition: Numeric.cpp:550
PViewDataList::NbTQ
int NbTQ
Definition: PViewDataList.h:32
PViewDataList::SI
std::vector< double > SI
Definition: PViewDataList.h:41
PViewDataList::NbVT
int NbVT
Definition: PViewDataList.h:30
PViewData::getName
virtual std::string getName()
Definition: PViewData.h:70
PViewDataList::VT
std::vector< double > VT
Definition: PViewDataList.h:31
GMSH_Lambda2Plugin
Definition: Lambda2.h:15
PViewDataList::SQ
std::vector< double > SQ
Definition: PViewDataList.h:33
PViewDataList::NbST
int NbST
Definition: PViewDataList.h:30
PViewDataList::NbSY
int NbSY
Definition: PViewDataList.h:42
GMSH_PostPlugin::getDataList
virtual PViewDataList * getDataList(PView *view, bool showError=true)
Definition: Plugin.cpp:107