gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
StreamLines.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 <cmath>
7 #include "GmshConfig.h"
8 #include "StreamLines.h"
9 #include "OctreePost.h"
10 #include "Context.h"
11 #include "PViewOptions.h"
12 
13 #if defined(HAVE_OPENGL)
14 #include "drawContext.h"
15 #endif
16 
17 StringXNumber StreamLinesOptions_Number[] = {
18  {GMSH_FULLRC, "X0", GMSH_StreamLinesPlugin::callbackX0, 0.},
19  {GMSH_FULLRC, "Y0", GMSH_StreamLinesPlugin::callbackY0, 0.},
20  {GMSH_FULLRC, "Z0", GMSH_StreamLinesPlugin::callbackZ0, 0.},
21  {GMSH_FULLRC, "X1", GMSH_StreamLinesPlugin::callbackX1, 1.},
22  {GMSH_FULLRC, "Y1", GMSH_StreamLinesPlugin::callbackY1, 0.},
23  {GMSH_FULLRC, "Z1", GMSH_StreamLinesPlugin::callbackZ1, 0.},
24  {GMSH_FULLRC, "X2", GMSH_StreamLinesPlugin::callbackX2, 0.},
25  {GMSH_FULLRC, "Y2", GMSH_StreamLinesPlugin::callbackY2, 1.},
26  {GMSH_FULLRC, "Z2", GMSH_StreamLinesPlugin::callbackZ2, 0.},
27  {GMSH_FULLRC, "NumPointsU", GMSH_StreamLinesPlugin::callbackU, 10},
28  {GMSH_FULLRC, "NumPointsV", GMSH_StreamLinesPlugin::callbackV, 1},
29  {GMSH_FULLRC, "DT", nullptr, .1},
30  {GMSH_FULLRC, "MaxIter", nullptr, 100},
31  {GMSH_FULLRC, "TimeStep", nullptr, 0},
32  {GMSH_FULLRC, "View", nullptr, -1.},
33  {GMSH_FULLRC, "OtherView", nullptr, -1.}};
34 
35 extern "C" {
36 GMSH_Plugin *GMSH_RegisterStreamLinesPlugin()
37 {
38  return new GMSH_StreamLinesPlugin();
39 }
40 }
41 
42 void GMSH_StreamLinesPlugin::draw(void *context)
43 {
44 #if defined(HAVE_OPENGL)
45  glColor4ubv((GLubyte *)&CTX::instance()->color.fg);
46  drawContext *ctx = (drawContext *)context;
47  double p[3];
48  for(int i = 0; i < getNbU(); ++i) {
49  for(int j = 0; j < getNbV(); ++j) {
50  getPoint(i, j, p);
51  ctx->drawSphere(CTX::instance()->pointSize, p[0], p[1], p[2], 1);
52  }
53  }
54 #endif
55 }
56 
57 double GMSH_StreamLinesPlugin::callback(int num, int action, double value,
58  double *opt, double step, double min,
59  double max)
60 {
61  switch(action) { // configure the input field
62  case 1: return step;
63  case 2: return min;
64  case 3: return max;
65  default: break;
66  }
67  *opt = value;
68  GMSH_Plugin::setDrawFunction(draw);
69  return 0.;
70 }
71 
72 double GMSH_StreamLinesPlugin::callbackX0(int num, int action, double value)
73 {
74  return callback(num, action, value, &StreamLinesOptions_Number[0].def,
75  CTX::instance()->lc / 100., -2 * CTX::instance()->lc,
76  2 * CTX::instance()->lc);
77 }
78 
79 double GMSH_StreamLinesPlugin::callbackY0(int num, int action, double value)
80 {
81  return callback(num, action, value, &StreamLinesOptions_Number[1].def,
82  CTX::instance()->lc / 100., -2 * CTX::instance()->lc,
83  2 * CTX::instance()->lc);
84 }
85 
86 double GMSH_StreamLinesPlugin::callbackZ0(int num, int action, double value)
87 {
88  return callback(num, action, value, &StreamLinesOptions_Number[2].def,
89  CTX::instance()->lc / 100., -2 * CTX::instance()->lc,
90  2 * CTX::instance()->lc);
91 }
92 
93 double GMSH_StreamLinesPlugin::callbackX1(int num, int action, double value)
94 {
95  return callback(num, action, value, &StreamLinesOptions_Number[3].def,
96  CTX::instance()->lc / 100., -2 * CTX::instance()->lc,
97  2 * CTX::instance()->lc);
98 }
99 
100 double GMSH_StreamLinesPlugin::callbackY1(int num, int action, double value)
101 {
102  return callback(num, action, value, &StreamLinesOptions_Number[4].def,
103  CTX::instance()->lc / 100., -2 * CTX::instance()->lc,
104  2 * CTX::instance()->lc);
105 }
106 
107 double GMSH_StreamLinesPlugin::callbackZ1(int num, int action, double value)
108 {
109  return callback(num, action, value, &StreamLinesOptions_Number[5].def,
110  CTX::instance()->lc / 100., -2 * CTX::instance()->lc,
111  2 * CTX::instance()->lc);
112 }
113 
114 double GMSH_StreamLinesPlugin::callbackX2(int num, int action, double value)
115 {
116  return callback(num, action, value, &StreamLinesOptions_Number[6].def,
117  CTX::instance()->lc / 100., -2 * CTX::instance()->lc,
118  2 * CTX::instance()->lc);
119 }
120 
121 double GMSH_StreamLinesPlugin::callbackY2(int num, int action, double value)
122 {
123  return callback(num, action, value, &StreamLinesOptions_Number[7].def,
124  CTX::instance()->lc / 100., -2 * CTX::instance()->lc,
125  2 * CTX::instance()->lc);
126 }
127 
128 double GMSH_StreamLinesPlugin::callbackZ2(int num, int action, double value)
129 {
130  return callback(num, action, value, &StreamLinesOptions_Number[8].def,
131  CTX::instance()->lc / 100., -2 * CTX::instance()->lc,
132  2 * CTX::instance()->lc);
133 }
134 
135 double GMSH_StreamLinesPlugin::callbackU(int num, int action, double value)
136 {
137  return callback(num, action, value, &StreamLinesOptions_Number[9].def, 1, 1,
138  100);
139 }
140 
141 double GMSH_StreamLinesPlugin::callbackV(int num, int action, double value)
142 {
143  return callback(num, action, value, &StreamLinesOptions_Number[10].def, 1, 1,
144  100);
145 }
146 
147 std::string GMSH_StreamLinesPlugin::getHelp() const
148 {
149  return "Plugin(StreamLines) computes stream lines "
150  "from the `TimeStep'-th time step of a vector "
151  "view `View' and optionally interpolates the "
152  "scalar view `OtherView' on the resulting stream "
153  "lines.\n\n"
154  "The plugin takes as input a grid defined by the "
155  "3 points (`X0',`Y0',`Z0') (origin), (`X1',`Y1',`Z1') "
156  "(axis of U) and (`X2',`Y2',`Z2') (axis of V).\n\n"
157  "The number of points along U and V that are to be "
158  "transported is set with the options `NumPointsU' "
159  "and `NumPointsV'. The equation\n\n"
160  "dX(t)/dt = V(x,y,z)\n\n"
161  "is then solved with the initial condition X(t=0) "
162  "chosen as the grid and with V(x,y,z) interpolated "
163  "on the vector view.\n\n"
164  "The time stepping scheme is a RK44 with step size "
165  "`DT' and `MaxIter' maximum number of iterations.\n\n"
166  "If `TimeStep' < 0, the plugin tries to compute "
167  "streamlines of the unsteady flow.\n\n"
168  "If `View' < 0, the plugin is run on the current view.\n\n"
169  "Plugin(StreamLines) creates one new list-based view. This "
170  "view contains multi-step vector points if `OtherView' "
171  "< 0, or single-step scalar lines if `OtherView' >= 0.";
172 }
173 
174 int GMSH_StreamLinesPlugin::getNbOptions() const
175 {
176  return sizeof(StreamLinesOptions_Number) / sizeof(StringXNumber);
177 }
178 
179 StringXNumber *GMSH_StreamLinesPlugin::getOption(int iopt)
180 {
181  return &StreamLinesOptions_Number[iopt];
182 }
183 
184 int GMSH_StreamLinesPlugin::getNbU()
185 {
186  return (int)StreamLinesOptions_Number[9].def;
187 }
188 
189 int GMSH_StreamLinesPlugin::getNbV()
190 {
191  return (int)StreamLinesOptions_Number[10].def;
192 }
193 
194 void GMSH_StreamLinesPlugin::getPoint(int iU, int iV, double *X)
195 {
196  double u = getNbU() > 1 ? (double)iU / (double)(getNbU() - 1.) : 0.;
197  double v = getNbV() > 1 ? (double)iV / (double)(getNbV() - 1.) : 0.;
198  X[0] =
199  StreamLinesOptions_Number[0].def +
200  u * (StreamLinesOptions_Number[3].def - StreamLinesOptions_Number[0].def) +
201  v * (StreamLinesOptions_Number[6].def - StreamLinesOptions_Number[0].def);
202  X[1] =
203  StreamLinesOptions_Number[1].def +
204  u * (StreamLinesOptions_Number[4].def - StreamLinesOptions_Number[1].def) +
205  v * (StreamLinesOptions_Number[7].def - StreamLinesOptions_Number[1].def);
206  X[2] =
207  StreamLinesOptions_Number[2].def +
208  u * (StreamLinesOptions_Number[5].def - StreamLinesOptions_Number[2].def) +
209  v * (StreamLinesOptions_Number[8].def - StreamLinesOptions_Number[2].def);
210 }
211 
212 PView *GMSH_StreamLinesPlugin::execute(PView *v)
213 {
214  double DT = StreamLinesOptions_Number[11].def;
215  int maxIter = (int)StreamLinesOptions_Number[12].def;
216  int timeStep = (int)StreamLinesOptions_Number[13].def;
217  int iView = (int)StreamLinesOptions_Number[14].def;
218  int otherView = (int)StreamLinesOptions_Number[15].def;
219 
220  PView *v1 = getView(iView, v);
221  if(!v1) return v;
222  PViewData *data1 = getPossiblyAdaptiveData(v1);
223 
224  PView *v2 = (otherView < 0) ? nullptr : getView(otherView, v);
225  PViewData *data2 = v2 ? getPossiblyAdaptiveData(v2) : nullptr;
226 
227  // sanity checks
228  if(timeStep > data1->getNumTimeSteps() - 1) {
229  Msg::Error("Invalid time step (%d) in view[%d]", v1->getIndex());
230  return v;
231  }
232 
233  OctreePost o1(v1);
234  double *val2 = nullptr;
235  OctreePost *o2 = nullptr;
236  if(data2) {
237  val2 = new double[data2->getNumTimeSteps()];
238  o2 = new OctreePost(v2);
239  }
240 
241  PView *v3 = new PView();
242  PViewDataList *data3 = getDataList(v3);
243 
244  const double b1 = 1. / 3., b2 = 2. / 3., b3 = 1. / 3., b4 = 1. / 6.;
245  const double a1 = 0.5, a2 = 0.5, a3 = 1., a4 = 1.;
246  double XINIT[3], X[3], DX[3], X1[3], X2[3], X3[3], X4[3];
247 
248  for(int i = 0; i < getNbU(); ++i) {
249  for(int j = 0; j < getNbV(); ++j) {
250  getPoint(i, j, XINIT);
251  getPoint(i, j, X);
252 
253  if(data2) { o2->searchScalar(X[0], X[1], X[2], val2, -1); }
254  else {
255  data3->NbVP++;
256  data3->VP.push_back(X[0]);
257  data3->VP.push_back(X[1]);
258  data3->VP.push_back(X[2]);
259  }
260 
261  int currentTimeStep = 0;
262 
263  for(int iter = 0; iter < maxIter; iter++) {
264  double XPREV[3] = {X[0], X[1], X[2]};
265 
266  if(timeStep < 0) {
267  double T0 = data1->getTime(0);
268  double currentT = T0 + DT * iter;
269  data3->Time.push_back(currentT);
270  for(; currentTimeStep < data1->getNumTimeSteps() - 1 &&
271  currentT > 0.5 * (data1->getTime(currentTimeStep) +
272  data1->getTime(currentTimeStep + 1));
273  currentTimeStep++)
274  ;
275  }
276  else {
277  currentTimeStep = timeStep;
278  }
279 
280  // dX/dt = V
281  // X1 = X + a1 * DT * V(X)
282  // X2 = X + a2 * DT * V(X1)
283  // X3 = X + a3 * DT * V(X2)
284  // X4 = X + a4 * DT * V(X3)
285  // X = X + b1 X1 + b2 X2 + b3 X3 + b4 x4
286  double val[3];
287  o1.searchVector(X[0], X[1], X[2], val, currentTimeStep);
288  for(int k = 0; k < 3; k++) X1[k] = X[k] + DT * val[k] * a1;
289  o1.searchVector(X1[0], X1[1], X1[2], val, currentTimeStep);
290  for(int k = 0; k < 3; k++) X2[k] = X[k] + DT * val[k] * a2;
291  o1.searchVector(X2[0], X2[1], X2[2], val, currentTimeStep);
292  for(int k = 0; k < 3; k++) X3[k] = X[k] + DT * val[k] * a3;
293  o1.searchVector(X3[0], X3[1], X3[2], val, currentTimeStep);
294  for(int k = 0; k < 3; k++) X4[k] = X[k] + DT * val[k] * a4;
295 
296  for(int k = 0; k < 3; k++)
297  X[k] += (b1 * (X1[k] - X[k]) + b2 * (X2[k] - X[k]) +
298  b3 * (X3[k] - X[k]) + b4 * (X4[k] - X[k]));
299  for(int k = 0; k < 3; k++) DX[k] = X[k] - XINIT[k];
300 
301  if(data2) {
302  data3->NbSL++;
303  data3->SL.push_back(XPREV[0]);
304  data3->SL.push_back(X[0]);
305  data3->SL.push_back(XPREV[1]);
306  data3->SL.push_back(X[1]);
307  data3->SL.push_back(XPREV[2]);
308  data3->SL.push_back(X[2]);
309  for(int k = 0; k < data2->getNumTimeSteps(); k++)
310  data3->SL.push_back(val2[k]);
311  o2->searchScalar(X[0], X[1], X[2], val2, -1);
312  for(int k = 0; k < data2->getNumTimeSteps(); k++)
313  data3->SL.push_back(val2[k]);
314  }
315  else {
316  data3->VP.push_back(DX[0]);
317  data3->VP.push_back(DX[1]);
318  data3->VP.push_back(DX[2]);
319  }
320  }
321  }
322  }
323 
324  if(data2) {
325  delete[] val2;
326  delete o2;
327  }
328  else {
329  v3->getOptions()->vectorType = PViewOptions::Displacement;
330  }
331 
332  data3->setName(data1->getName() + "_StreamLines");
333  data3->setFileName(data1->getName() + "_StreamLines.pos");
334  data3->finalize();
335 
336  return v3;
337 }
GMSH_StreamLinesPlugin::execute
PView * execute(PView *)
Definition: StreamLines.cpp:212
GMSH_StreamLinesPlugin::callbackZ2
static double callbackZ2(int, int, double)
Definition: StreamLines.cpp:128
GMSH_StreamLinesPlugin::callbackY1
static double callbackY1(int, int, double)
Definition: StreamLines.cpp:100
PView
Definition: PView.h:27
GMSH_StreamLinesPlugin::callbackY2
static double callbackY2(int, int, double)
Definition: StreamLines.cpp:121
OctreePost::searchScalar
bool searchScalar(double x, double y, double z, double *values, int step=-1, double *size=nullptr, int qn=0, double *qx=nullptr, double *qy=nullptr, double *qz=nullptr, bool grad=false, int dim=-1)
Definition: OctreePost.cpp:578
StreamLines.h
GMSH_Plugin
Definition: Plugin.h:26
StringXNumber::def
double def
Definition: Options.h:922
PViewData::getNumTimeSteps
virtual int getNumTimeSteps()=0
OctreePost.h
PViewDataList
Definition: PViewDataList.h:17
GMSH_Plugin::draw
static void(* draw)(void *)
Definition: Plugin.h:77
a4
#define a4
Definition: GaussQuadratureTet.cpp:11
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
StringXNumber
Definition: Options.h:918
GMSH_StreamLinesPlugin::callbackX2
static double callbackX2(int, int, double)
Definition: StreamLines.cpp:114
OctreePost
Definition: BoundaryLayers.cpp:23
PView::getIndex
int getIndex()
Definition: PView.h:92
GMSH_StreamLinesPlugin::callbackV
static double callbackV(int, int, double)
Definition: StreamLines.cpp:141
PViewData::setFileName
virtual void setFileName(const std::string &val)
Definition: PViewData.h:75
CTX::fg
unsigned int fg
Definition: Context.h:358
CTX::instance
static CTX * instance()
Definition: Context.cpp:122
GMSH_RegisterStreamLinesPlugin
GMSH_Plugin * GMSH_RegisterStreamLinesPlugin()
Definition: StreamLines.cpp:36
PViewData::getTime
virtual double getTime(int step)
Definition: PViewData.h:94
GMSH_StreamLinesPlugin::callbackZ1
static double callbackZ1(int, int, double)
Definition: StreamLines.cpp:107
PViewDataList::VP
std::vector< double > VP
Definition: PViewDataList.h:27
GMSH_FULLRC
#define GMSH_FULLRC
Definition: Options.h:20
a1
const double a1
Definition: GaussQuadratureHex.cpp:10
PViewDataList::SL
std::vector< double > SL
Definition: PViewDataList.h:29
GMSH_StreamLinesPlugin::callbackU
static double callbackU(int, int, double)
Definition: StreamLines.cpp:135
PViewOptions.h
b4
#define b4
Definition: GaussQuadratureTet.cpp:12
drawContext
Definition: drawContext.h:120
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
GMSH_StreamLinesPlugin::callbackX0
static double callbackX0(int, int, double)
Definition: StreamLines.cpp:72
PViewOptions::Displacement
@ Displacement
Definition: PViewOptions.h:25
GMSH_StreamLinesPlugin::callbackY0
static double callbackY0(int, int, double)
Definition: StreamLines.cpp:79
a2
const double a2
Definition: GaussQuadratureHex.cpp:12
PViewDataList::NbVP
int NbVP
Definition: PViewDataList.h:26
PViewData
Definition: PViewData.h:29
GMSH_StreamLinesPlugin::getNbOptions
int getNbOptions() const
Definition: StreamLines.cpp:174
b1
const double b1
Definition: GaussQuadratureHex.cpp:14
GMSH_StreamLinesPlugin
Definition: StreamLines.h:15
GMSH_StreamLinesPlugin::callback
static double callback(int num, int action, double value, double *opt, double step, double min, double max)
Definition: StreamLines.cpp:57
GMSH_StreamLinesPlugin::getPoint
static void getPoint(int iU, int iV, double *X)
Definition: StreamLines.cpp:194
PViewDataList::NbSL
int NbSL
Definition: PViewDataList.h:28
Context.h
GMSH_StreamLinesPlugin::getOption
StringXNumber * getOption(int iopt)
Definition: StreamLines.cpp:179
OctreePost::searchVector
bool searchVector(double x, double y, double z, double *values, int step=-1, double *size=nullptr, int qn=0, double *qx=nullptr, double *qy=nullptr, double *qz=nullptr, bool grad=false, int dim=-1)
Definition: OctreePost.cpp:634
GMSH_Plugin::setDrawFunction
static void setDrawFunction(void(*fct)(void *))
Definition: Plugin.cpp:21
GMSH_PostPlugin::getView
virtual PView * getView(int index, PView *view)
Definition: Plugin.cpp:81
GMSH_PostPlugin::getPossiblyAdaptiveData
virtual PViewData * getPossiblyAdaptiveData(PView *view)
Definition: Plugin.cpp:94
PViewDataList::Time
std::vector< double > Time
Definition: PViewDataList.h:25
PView::getOptions
PViewOptions * getOptions()
Definition: PView.h:81
PViewOptions::vectorType
int vectorType
Definition: PViewOptions.h:60
GMSH_StreamLinesPlugin::callbackX1
static double callbackX1(int, int, double)
Definition: StreamLines.cpp:93
GMSH_StreamLinesPlugin::callbackZ0
static double callbackZ0(int, int, double)
Definition: StreamLines.cpp:86
GMSH_StreamLinesPlugin::getNbV
static int getNbV()
Definition: StreamLines.cpp:189
StreamLinesOptions_Number
StringXNumber StreamLinesOptions_Number[]
Definition: StreamLines.cpp:17
PViewData::getName
virtual std::string getName()
Definition: PViewData.h:70
drawContext::drawSphere
void drawSphere(double R, double x, double y, double z, int n1, int n2, int light)
Definition: drawGlyph.cpp:390
GMSH_StreamLinesPlugin::getNbU
static int getNbU()
Definition: StreamLines.cpp:184
GMSH_PostPlugin::getDataList
virtual PViewDataList * getDataList(PView *view, bool showError=true)
Definition: Plugin.cpp:107
GMSH_StreamLinesPlugin::getHelp
std::string getHelp() const
Definition: StreamLines.cpp:147
DT
Definition: DivideAndConquer.h:50
drawContext.h