gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
InnerVertexPlacement.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 // Contributed by Amaury Johnen
7 
8 #include "InnerVertexPlacement.h"
9 #include "GmshDefines.h"
10 #include "pointsGenerators.h"
11 
12 namespace {
13  typedef std::pair<std::pair<int, int>, int> mytuple;
14 
15  mytuple make_mytuple(int a, int b, int c)
16  {
17  return std::make_pair(std::make_pair(a, b), c);
18  }
19 
20  void double2int(const fullMatrix<double> &doubleMat, fullMatrix<int> &intMat)
21  {
22  intMat.resize(doubleMat.size1(), doubleMat.size2(), false);
23  for(int i = 0; i < intMat.size1(); ++i) {
24  for(int j = 0; j < intMat.size2(); ++j) {
25  intMat(i, j) = static_cast<int>(doubleMat(i, j) + .5);
26  }
27  }
28  }
29 
30  std::map<int, fullMatrix<double> *> storedMatrices[6];
31 } // namespace
32 
33 fullMatrix<double> *getInnerVertexPlacement(int type, int order)
34 {
35  if(type < 3 || type > 8) return nullptr;
36  std::map<int, fullMatrix<double> *>::iterator it;
37  it = storedMatrices[type - 3].find(order);
38  if(it != storedMatrices[type - 3].end()) { return it->second; }
39  else {
40  fullMatrix<double> *matrix = new fullMatrix<double>();
41  switch(type) {
42  case TYPE_TRI:
43  *matrix = gmshGenerateInnerVertexPlacementTriangle(order);
44  break;
45  case TYPE_QUA:
46  *matrix = gmshGenerateInnerVertexPlacementQuadrangle(order);
47  break;
48  case TYPE_TET:
49  *matrix = gmshGenerateInnerVertexPlacementTetrahedron(order);
50  break;
51  case TYPE_PRI:
52  *matrix = gmshGenerateInnerVertexPlacementPrism(order);
53  break;
54  case TYPE_HEX:
55  *matrix = gmshGenerateInnerVertexPlacementHexahedron(order);
56  break;
57  case TYPE_PYR:
58  *matrix = gmshGenerateInnerVertexPlacementPyramid(order);
59  break;
60  }
61  storedMatrices[type - 3][order] = matrix;
62  return matrix;
63  }
64 }
65 
66 fullMatrix<double> gmshGenerateInnerVertexPlacementTriangle(int order)
67 {
68  if(order < 3) return fullMatrix<double>(0, 0);
69 
70  const int szInc = 3 * order;
71  const int szComp = (order + 1) * (order + 2) / 2;
72  const int szInt = szComp - szInc;
73 
74  fullMatrix<double> monomials = gmshGenerateMonomialsTriangle(order, false);
75  fullMatrix<int> coordinates;
76  double2int(monomials, coordinates);
77 
78  std::map<std::pair<int, int>, int> coord2idx;
79  for(int i = 0; i < szInc; ++i) {
80  int u = coordinates(i, 0);
81  int v = coordinates(i, 1);
82  coord2idx[std::make_pair(u, v)] = i;
83  }
84 
85  int &n = order;
86  fullMatrix<double> M(szInt, szInc, true);
87  for(int i = 0; i < szInt; ++i) {
88  int u = coordinates(szInc + i, 0);
89  int v = coordinates(szInc + i, 1);
90  int q = n - u - v;
91  double xi = (double)u / n;
92  double eta = (double)v / n;
93  double rho = (double)q / n;
94 
95  M(i, coord2idx[std::make_pair(u + v, 0)]) += xi;
96  M(i, coord2idx[std::make_pair(u + q, v)]) += xi;
97  M(i, coord2idx[std::make_pair(0, v + u)]) += eta;
98  M(i, coord2idx[std::make_pair(u, v + q)]) += eta;
99  M(i, coord2idx[std::make_pair(0, v)]) += rho;
100  M(i, coord2idx[std::make_pair(u, 0)]) += rho;
101  M(i, coord2idx[std::make_pair(n, 0)]) -= xi;
102  M(i, coord2idx[std::make_pair(0, n)]) -= eta;
103  M(i, coord2idx[std::make_pair(0, 0)]) -= rho;
104  }
105  return M;
106 }
107 
108 fullMatrix<double> gmshGenerateInnerVertexPlacementQuadrangle(int order)
109 {
110  if(order < 2) return fullMatrix<double>(0, 0);
111 
112  const int szInc = 4 * order;
113  const int szComp = (order + 1) * (order + 1);
114  const int szInt = szComp - szInc;
115 
116  fullMatrix<double> monomials = gmshGenerateMonomialsQuadrangle(order, false);
117  fullMatrix<int> coordinates;
118  double2int(monomials, coordinates);
119 
120  std::map<std::pair<int, int>, int> coord2idx;
121  for(int i = 0; i < szInc; ++i) {
122  int u = coordinates(i, 0);
123  int v = coordinates(i, 1);
124  coord2idx[std::make_pair(u, v)] = i;
125  }
126 
127  int &n = order;
128  fullMatrix<double> M(szInt, szInc, true);
129  for(int i = 0; i < szInt; ++i) {
130  int u = coordinates(szInc + i, 0);
131  int v = coordinates(szInc + i, 1);
132  double xi = (double)u / n;
133  double eta = (double)v / n;
134 
135  M(i, coord2idx[std::make_pair(0, v)]) += 1 - xi;
136  M(i, coord2idx[std::make_pair(n, v)]) += xi;
137  M(i, coord2idx[std::make_pair(u, 0)]) += 1 - eta;
138  M(i, coord2idx[std::make_pair(u, n)]) += eta;
139  M(i, coord2idx[std::make_pair(0, 0)]) -= (1 - xi) * (1 - eta);
140  M(i, coord2idx[std::make_pair(n, 0)]) -= xi * (1 - eta);
141  M(i, coord2idx[std::make_pair(n, n)]) -= xi * eta;
142  M(i, coord2idx[std::make_pair(0, n)]) -= (1 - xi) * eta;
143  }
144  return M;
145 }
146 
147 fullMatrix<double> gmshGenerateInnerVertexPlacementTetrahedron(int order)
148 {
149  if(order < 4) return fullMatrix<double>(0, 0);
150 
151  const int szInt = (order - 3) * (order - 2) * (order - 1) / 6;
152  const int szComp = (order + 1) * (order + 2) * (order + 3) / 6;
153  const int szInc = szComp - szInt;
154 
155  fullMatrix<double> monomials = gmshGenerateMonomialsTetrahedron(order, false);
156  fullMatrix<int> coordinates;
157  double2int(monomials, coordinates);
158 
159  std::map<mytuple, int> coord2idx;
160  for(int i = 0; i < szInc; ++i) {
161  int u = coordinates(i, 0);
162  int v = coordinates(i, 1);
163  int w = coordinates(i, 2);
164  coord2idx[make_mytuple(u, v, w)] = i;
165  }
166 
167  int &n = order;
168  fullMatrix<double> M(szInt, szInc, true);
169  for(int i = 0; i < szInt; ++i) {
170  int u = coordinates(szInc + i, 0);
171  int v = coordinates(szInc + i, 1);
172  int w = coordinates(szInc + i, 2);
173  int q = n - u - v - w;
174  double xi = (double)u / n;
175  double eta = (double)v / n;
176  double zeta = (double)w / n;
177  double rho = (double)q / n;
178 
179  M(i, coord2idx[make_mytuple(u + v, 0, w)]) += xi;
180  M(i, coord2idx[make_mytuple(u + w, v, 0)]) += xi;
181  M(i, coord2idx[make_mytuple(u + q, v, w)]) += xi;
182  M(i, coord2idx[make_mytuple(0, v + u, w)]) += eta;
183  M(i, coord2idx[make_mytuple(u, v + w, 0)]) += eta;
184  M(i, coord2idx[make_mytuple(u, v + q, w)]) += eta;
185  M(i, coord2idx[make_mytuple(0, v, w + u)]) += zeta;
186  M(i, coord2idx[make_mytuple(u, 0, w + v)]) += zeta;
187  M(i, coord2idx[make_mytuple(u, v, w + q)]) += zeta;
188  M(i, coord2idx[make_mytuple(0, v, w)]) += rho;
189  M(i, coord2idx[make_mytuple(u, 0, w)]) += rho;
190  M(i, coord2idx[make_mytuple(u, v, 0)]) += rho;
191 
192  M(i, coord2idx[make_mytuple(n - q, 0, 0)]) -= xi;
193  M(i, coord2idx[make_mytuple(n - w, 0, w)]) -= xi;
194  M(i, coord2idx[make_mytuple(n - v, v, 0)]) -= xi;
195  M(i, coord2idx[make_mytuple(0, n - q, 0)]) -= eta;
196  M(i, coord2idx[make_mytuple(0, n - w, w)]) -= eta;
197  M(i, coord2idx[make_mytuple(u, n - u, 0)]) -= eta;
198  M(i, coord2idx[make_mytuple(0, 0, n - q)]) -= zeta;
199  M(i, coord2idx[make_mytuple(0, v, n - v)]) -= zeta;
200  M(i, coord2idx[make_mytuple(u, 0, n - u)]) -= zeta;
201  M(i, coord2idx[make_mytuple(0, 0, w)]) -= rho;
202  M(i, coord2idx[make_mytuple(0, v, 0)]) -= rho;
203  M(i, coord2idx[make_mytuple(u, 0, 0)]) -= rho;
204 
205  M(i, coord2idx[make_mytuple(n, 0, 0)]) += xi;
206  M(i, coord2idx[make_mytuple(0, n, 0)]) += eta;
207  M(i, coord2idx[make_mytuple(0, 0, n)]) += zeta;
208  M(i, coord2idx[make_mytuple(0, 0, 0)]) += rho;
209  }
210  return M;
211 }
212 
213 fullMatrix<double> gmshGenerateInnerVertexPlacementHexahedron(int order)
214 {
215  if(order < 2) return fullMatrix<double>(0, 0);
216 
217  const int szInt = (order - 1) * (order - 1) * (order - 1);
218  const int szComp = (order + 1) * (order + 1) * (order + 1);
219  const int szInc = szComp - szInt;
220 
221  fullMatrix<double> monomials = gmshGenerateMonomialsHexahedron(order, false);
222  fullMatrix<int> coordinates;
223  double2int(monomials, coordinates);
224 
225  std::map<mytuple, int> coord2idx;
226  for(int i = 0; i < szInc; ++i) {
227  int u = coordinates(i, 0);
228  int v = coordinates(i, 1);
229  int w = coordinates(i, 2);
230  coord2idx[make_mytuple(u, v, w)] = i;
231  }
232 
233  int &n = order;
234  fullMatrix<double> M(szInt, szInc, true);
235  for(int i = 0; i < szInt; ++i) {
236  int u = coordinates(szInc + i, 0);
237  int v = coordinates(szInc + i, 1);
238  int w = coordinates(szInc + i, 2);
239  double xi = (double)u / n;
240  double eta = (double)v / n;
241  double zeta = (double)w / n;
242 
243  M(i, coord2idx[make_mytuple(0, v, w)]) += 1 - xi;
244  M(i, coord2idx[make_mytuple(n, v, w)]) += xi;
245  M(i, coord2idx[make_mytuple(u, 0, w)]) += 1 - eta;
246  M(i, coord2idx[make_mytuple(u, n, w)]) += eta;
247  M(i, coord2idx[make_mytuple(u, v, 0)]) += 1 - zeta;
248  M(i, coord2idx[make_mytuple(u, v, n)]) += zeta;
249 
250  M(i, coord2idx[make_mytuple(0, 0, w)]) -= (1 - xi) * (1 - eta);
251  M(i, coord2idx[make_mytuple(0, n, w)]) -= (1 - xi) * eta;
252  M(i, coord2idx[make_mytuple(n, 0, w)]) -= xi * (1 - eta);
253  M(i, coord2idx[make_mytuple(n, n, w)]) -= xi * eta;
254  M(i, coord2idx[make_mytuple(u, 0, 0)]) -= (1 - eta) * (1 - zeta);
255  M(i, coord2idx[make_mytuple(u, 0, n)]) -= (1 - eta) * zeta;
256  M(i, coord2idx[make_mytuple(u, n, 0)]) -= eta * (1 - zeta);
257  M(i, coord2idx[make_mytuple(u, n, n)]) -= eta * zeta;
258  M(i, coord2idx[make_mytuple(0, v, 0)]) -= (1 - zeta) * (1 - xi);
259  M(i, coord2idx[make_mytuple(n, v, 0)]) -= (1 - zeta) * xi;
260  M(i, coord2idx[make_mytuple(0, v, n)]) -= zeta * (1 - xi);
261  M(i, coord2idx[make_mytuple(n, v, n)]) -= zeta * xi;
262 
263  M(i, coord2idx[make_mytuple(0, 0, 0)]) += (1 - xi) * (1 - eta) * (1 - zeta);
264  M(i, coord2idx[make_mytuple(0, 0, n)]) += (1 - xi) * (1 - eta) * zeta;
265  M(i, coord2idx[make_mytuple(0, n, 0)]) += (1 - xi) * eta * (1 - zeta);
266  M(i, coord2idx[make_mytuple(0, n, n)]) += (1 - xi) * eta * zeta;
267  M(i, coord2idx[make_mytuple(n, 0, 0)]) += xi * (1 - eta) * (1 - zeta);
268  M(i, coord2idx[make_mytuple(n, 0, n)]) += xi * (1 - eta) * zeta;
269  M(i, coord2idx[make_mytuple(n, n, 0)]) += xi * eta * (1 - zeta);
270  M(i, coord2idx[make_mytuple(n, n, n)]) += xi * eta * zeta;
271  }
272  return M;
273 }
274 
275 fullMatrix<double> gmshGenerateInnerVertexPlacementPrism(int order)
276 {
277  if(order < 3) return fullMatrix<double>(0, 0);
278 
279  const int szInt = (order - 1) * (order - 2) * (order - 1) / 2;
280  const int szComp = (order + 1) * (order + 1) * (order + 2) / 2;
281  const int szInc = szComp - szInt;
282 
283  fullMatrix<double> monomials = gmshGenerateMonomialsPrism(order, false);
284  fullMatrix<int> coordinates;
285  double2int(monomials, coordinates);
286 
287  std::map<mytuple, int> coord2idx;
288  for(int i = 0; i < szInc; ++i) {
289  int u = coordinates(i, 0);
290  int v = coordinates(i, 1);
291  int w = coordinates(i, 2);
292  coord2idx[make_mytuple(u, v, w)] = i;
293  }
294 
295  int &n = order;
296  fullMatrix<double> M(szInt, szInc, true);
297  for(int i = 0; i < szInt; ++i) {
298  int u = coordinates(szInc + i, 0);
299  int v = coordinates(szInc + i, 1);
300  int w = coordinates(szInc + i, 2);
301  int q = n - u - v;
302  double xi = (double)u / n;
303  double eta = (double)v / n;
304  double zeta = (double)w / n;
305  double rho = (double)q / n;
306 
307  M(i, coord2idx[make_mytuple(u + v, 0, w)]) += xi;
308  M(i, coord2idx[make_mytuple(u + q, v, w)]) += xi;
309  M(i, coord2idx[make_mytuple(0, v + u, w)]) += eta;
310  M(i, coord2idx[make_mytuple(u, v + q, w)]) += eta;
311  M(i, coord2idx[make_mytuple(0, v, w)]) += rho;
312  M(i, coord2idx[make_mytuple(u, 0, w)]) += rho;
313  M(i, coord2idx[make_mytuple(u, v, 0)]) += 1 - zeta;
314  M(i, coord2idx[make_mytuple(u, v, n)]) += zeta;
315 
316  M(i, coord2idx[make_mytuple(n, 0, w)]) -= xi;
317  M(i, coord2idx[make_mytuple(0, n, w)]) -= eta;
318  M(i, coord2idx[make_mytuple(0, 0, w)]) -= rho;
319  M(i, coord2idx[make_mytuple(u + v, 0, 0)]) -= xi * (1 - zeta);
320  M(i, coord2idx[make_mytuple(u + q, v, 0)]) -= xi * (1 - zeta);
321  M(i, coord2idx[make_mytuple(0, v + u, 0)]) -= eta * (1 - zeta);
322  M(i, coord2idx[make_mytuple(u, v + q, 0)]) -= eta * (1 - zeta);
323  M(i, coord2idx[make_mytuple(0, v, 0)]) -= rho * (1 - zeta);
324  M(i, coord2idx[make_mytuple(u, 0, 0)]) -= rho * (1 - zeta);
325  M(i, coord2idx[make_mytuple(u + v, 0, n)]) -= xi * zeta;
326  M(i, coord2idx[make_mytuple(u + q, v, n)]) -= xi * zeta;
327  M(i, coord2idx[make_mytuple(0, v + u, n)]) -= eta * zeta;
328  M(i, coord2idx[make_mytuple(u, v + q, n)]) -= eta * zeta;
329  M(i, coord2idx[make_mytuple(0, v, n)]) -= rho * zeta;
330  M(i, coord2idx[make_mytuple(u, 0, n)]) -= rho * zeta;
331 
332  M(i, coord2idx[make_mytuple(n, 0, 0)]) += xi * (1 - zeta);
333  M(i, coord2idx[make_mytuple(0, n, 0)]) += eta * (1 - zeta);
334  M(i, coord2idx[make_mytuple(0, 0, 0)]) += rho * (1 - zeta);
335  M(i, coord2idx[make_mytuple(n, 0, n)]) += xi * zeta;
336  M(i, coord2idx[make_mytuple(0, n, n)]) += eta * zeta;
337  M(i, coord2idx[make_mytuple(0, 0, n)]) += rho * zeta;
338  }
339  return M;
340 }
341 
342 fullMatrix<double> gmshGenerateInnerVertexPlacementPyramid(int order)
343 {
344  if(order < 3) return fullMatrix<double>(0, 0);
345 
346  const int szInt = (order - 2) * ((order - 2) + 1) * (2 * (order - 2) + 1) / 6;
347  const int szComp =
348  (order + 1) * ((order + 1) + 1) * (2 * (order + 1) + 1) / 6;
349  const int szInc = szComp - szInt;
350 
351  fullMatrix<double> monomials = gmshGenerateMonomialsPyramid(order, false);
352  fullMatrix<int> coordinates;
353  double2int(monomials, coordinates);
354 
355  std::map<mytuple, int> coord2idx;
356  for(int i = 0; i < szInc; ++i) {
357  int u = coordinates(i, 0);
358  int v = coordinates(i, 1);
359  int w = order - coordinates(i, 2);
360  coord2idx[make_mytuple(u, v, w)] = i;
361  }
362 
363  int &n = order;
364  fullMatrix<double> M(szInt, szInc, true);
365  for(int i = 0; i < szInt; ++i) {
366  int u = coordinates(szInc + i, 0);
367  int v = coordinates(szInc + i, 1);
368  int w = order - coordinates(szInc + i, 2);
369  int q = n - u - w;
370  int r = n - v - w;
371  double xi = (double)u / n;
372  double eta = (double)v / n;
373  double rho = (double)q / n;
374  double tau = (double)r / n;
375  double xip = (double)u / (n - w);
376  double etap = (double)v / (n - w);
377 
378  // lateral faces
379  M(i, coord2idx[make_mytuple(0, v, w)]) += 1 - xip;
380  M(i, coord2idx[make_mytuple(u, 0, w)]) += 1 - etap;
381  M(i, coord2idx[make_mytuple(n - w, v, w)]) += xip;
382  M(i, coord2idx[make_mytuple(u, n - w, w)]) += etap;
383  // basis face
384  M(i, coord2idx[make_mytuple(u, v, 0)]) += rho * tau;
385  M(i, coord2idx[make_mytuple(u + w, v, 0)]) += xi * tau;
386  M(i, coord2idx[make_mytuple(u + w, v + w, 0)]) += xi * eta;
387  M(i, coord2idx[make_mytuple(u, v + w, 0)]) += rho * eta;
388  // vertical edges
389  M(i, coord2idx[make_mytuple(0, 0, w)]) -= (1 - xip) * (1 - etap);
390  M(i, coord2idx[make_mytuple(n - w, 0, w)]) -= xip * (1 - etap);
391  M(i, coord2idx[make_mytuple(n - w, n - w, w)]) -= xip * etap;
392  M(i, coord2idx[make_mytuple(0, n - w, w)]) -= (1 - xip) * etap;
393  // basis edges
394  M(i, coord2idx[make_mytuple(0, v, 0)]) -= rho * tau;
395  M(i, coord2idx[make_mytuple(u, 0, 0)]) -= rho * tau;
396  M(i, coord2idx[make_mytuple(u + w, 0, 0)]) -= xi * tau;
397  M(i, coord2idx[make_mytuple(n, v, 0)]) -= xi * tau;
398  M(i, coord2idx[make_mytuple(n, v + w, 0)]) -= xi * eta;
399  M(i, coord2idx[make_mytuple(u + w, n, 0)]) -= xi * eta;
400  M(i, coord2idx[make_mytuple(u, n, 0)]) -= rho * eta;
401  M(i, coord2idx[make_mytuple(0, v + w, 0)]) -= rho * eta;
402  // basis corners
403  M(i, coord2idx[make_mytuple(0, 0, 0)]) += rho * tau;
404  M(i, coord2idx[make_mytuple(n, 0, 0)]) += xi * tau;
405  M(i, coord2idx[make_mytuple(n, n, 0)]) += xi * eta;
406  M(i, coord2idx[make_mytuple(0, n, 0)]) += rho * eta;
407  }
408  return M;
409 }
410 
411 fullMatrix<double> gmshGenerateInnerVertexPlacementTriangleLinear(int order,
412  int dir)
413 {
414  if(order < 3) return fullMatrix<double>(0, 0);
415 
416  const int szInc = 3 * order;
417  const int szComp = (order + 1) * (order + 2) / 2;
418  const int szInt = szComp - szInc;
419 
420  fullMatrix<double> monomials = gmshGenerateMonomialsTriangle(order, false);
421  fullMatrix<int> coordinates;
422  double2int(monomials, coordinates);
423 
424  std::map<std::pair<int, int>, int> coord2idx;
425  for(int i = 0; i < szInc; ++i) {
426  int u = coordinates(i, 0);
427  int v = coordinates(i, 1);
428  coord2idx[std::make_pair(u, v)] = i;
429  }
430 
431  int &n = order;
432  fullMatrix<double> M(szInt, szInc, true);
433  for(int i = 0; i < szInt; ++i) {
434  int u = coordinates(szInc + i, 0);
435  int v = coordinates(szInc + i, 1);
436  double mu;
437 
438  switch(dir) {
439  case 0:
440  mu = (double)u / (n - v);
441  M(i, coord2idx[std::make_pair(0, v)]) += 1 - mu;
442  M(i, coord2idx[std::make_pair(n - v, v)]) += mu;
443  break;
444  case 1:
445  mu = (double)v / (u + v);
446  M(i, coord2idx[std::make_pair(u + v, 0)]) += 1 - mu;
447  M(i, coord2idx[std::make_pair(0, u + v)]) += mu;
448  break;
449  default:
450  case 2:
451  mu = (double)v / (n - u);
452  M(i, coord2idx[std::make_pair(u, 0)]) += 1 - mu;
453  M(i, coord2idx[std::make_pair(u, n - u)]) += mu;
454  break;
455  }
456  }
457 
458  return M;
459 }
460 
461 fullMatrix<double> gmshGenerateInnerVertexPlacementQuadrangleLinear(int order)
462 {
463  if(order < 2) return fullMatrix<double>(0, 0);
464 
465  const int szInc = 4 * order;
466  const int szComp = (order + 1) * (order + 1);
467  const int szInt = szComp - szInc;
468 
469  fullMatrix<double> monomials = gmshGenerateMonomialsQuadrangle(order, false);
470  fullMatrix<int> coordinates;
471  double2int(monomials, coordinates);
472 
473  std::map<std::pair<int, int>, int> coord2idx;
474  for(int i = 0; i < szInc; ++i) {
475  int u = coordinates(i, 0);
476  int v = coordinates(i, 1);
477  coord2idx[std::make_pair(u, v)] = i;
478  }
479 
480  int &n = order;
481  fullMatrix<double> M(szInt, szInc, true);
482  for(int i = 0; i < szInt; ++i) {
483  int u = coordinates(szInc + i, 0);
484  int v = coordinates(szInc + i, 1);
485  double eta = (double)v / n;
486 
487  M(i, coord2idx[std::make_pair(u, 0)]) += 1 - eta;
488  M(i, coord2idx[std::make_pair(u, n)]) += eta;
489  }
490  return M;
491 }
492 
493 fullMatrix<double> gmshGenerateInnerVertexPlacementTetrahedronLinear(int order,
494  int dir)
495 {
496  // 'dir' corresponds to the number of edge
497  if(order < 4) return fullMatrix<double>(0, 0);
498 
499  const int szInt = (order - 3) * (order - 2) * (order - 1) / 6;
500  const int szComp = (order + 1) * (order + 2) * (order + 3) / 6;
501  const int szInc = szComp - szInt;
502 
503  fullMatrix<double> monomials = gmshGenerateMonomialsTetrahedron(order, false);
504  fullMatrix<int> coordinates;
505  double2int(monomials, coordinates);
506 
507  std::map<mytuple, int> coord2idx;
508  for(int i = 0; i < szInc; ++i) {
509  int u = coordinates(i, 0);
510  int v = coordinates(i, 1);
511  int w = coordinates(i, 2);
512  coord2idx[make_mytuple(u, v, w)] = i;
513  }
514 
515  int &n = order;
516  fullMatrix<double> M(szInt, szInc, true);
517  for(int i = 0; i < szInt; ++i) {
518  int u = coordinates(szInc + i, 0);
519  int v = coordinates(szInc + i, 1);
520  int w = coordinates(szInc + i, 2);
521  double mu;
522 
523  switch(dir) {
524  case 0:
525  mu = (double)u / (n - v - w);
526  M(i, coord2idx[make_mytuple(0, v, w)]) += 1 - mu;
527  M(i, coord2idx[make_mytuple(n - v - w, v, w)]) += mu;
528  break;
529  case 2:
530  mu = (double)v / (n - u - w);
531  M(i, coord2idx[make_mytuple(u, 0, w)]) += 1 - mu;
532  M(i, coord2idx[make_mytuple(u, n - u - w, w)]) += mu;
533  break;
534  default:
535  case 3:
536  mu = (double)w / (n - u - v);
537  M(i, coord2idx[make_mytuple(u, v, 0)]) += 1 - mu;
538  M(i, coord2idx[make_mytuple(u, v, n - u - v)]) += mu;
539  break;
540  case 1:
541  mu = (double)v / (u + v);
542  M(i, coord2idx[make_mytuple(u + v, 0, w)]) += 1 - mu;
543  M(i, coord2idx[make_mytuple(0, u + v, w)]) += mu;
544  break;
545  case 4:
546  mu = (double)w / (v + w);
547  M(i, coord2idx[make_mytuple(u, v + w, 0)]) += 1 - mu;
548  M(i, coord2idx[make_mytuple(u, 0, v + w)]) += mu;
549  break;
550  case 5:
551  mu = (double)w / (u + w);
552  M(i, coord2idx[make_mytuple(u + w, v, 0)]) += 1 - mu;
553  M(i, coord2idx[make_mytuple(0, v, u + w)]) += mu;
554  break;
555  }
556  }
557  return M;
558 }
559 
560 fullMatrix<double> gmshGenerateInnerVertexPlacementHexahedronLinear(int order,
561  int dir)
562 {
563  if(order < 2) return fullMatrix<double>(0, 0);
564 
565  const int szInt = (order - 1) * (order - 1) * (order - 1);
566  const int szComp = (order + 1) * (order + 1) * (order + 1);
567  const int szInc = szComp - szInt;
568 
569  fullMatrix<double> monomials = gmshGenerateMonomialsHexahedron(order, false);
570  fullMatrix<int> coordinates;
571  double2int(monomials, coordinates);
572 
573  std::map<mytuple, int> coord2idx;
574  for(int i = 0; i < szInc; ++i) {
575  int u = coordinates(i, 0);
576  int v = coordinates(i, 1);
577  int w = coordinates(i, 2);
578  coord2idx[make_mytuple(u, v, w)] = i;
579  }
580 
581  int &n = order;
582  fullMatrix<double> M(szInt, szInc, true);
583  for(int i = 0; i < szInt; ++i) {
584  int u = coordinates(szInc + i, 0);
585  int v = coordinates(szInc + i, 1);
586  int w = coordinates(szInc + i, 2);
587  double eta;
588 
589  switch(dir) {
590  case 0:
591  eta = (double)u / n;
592  M(i, coord2idx[make_mytuple(0, v, w)]) += 1 - eta;
593  M(i, coord2idx[make_mytuple(n, v, w)]) += eta;
594  break;
595  case 1:
596  eta = (double)v / n;
597  M(i, coord2idx[make_mytuple(u, 0, w)]) += 1 - eta;
598  M(i, coord2idx[make_mytuple(u, n, w)]) += eta;
599  break;
600  default:
601  case 2:
602  eta = (double)w / n;
603  M(i, coord2idx[make_mytuple(u, v, 0)]) += 1 - eta;
604  M(i, coord2idx[make_mytuple(u, v, n)]) += eta;
605  break;
606  }
607  }
608  return M;
609 }
610 
611 fullMatrix<double> gmshGenerateInnerVertexPlacementPrismLinear(int order,
612  int dir)
613 {
614  if(order < 3) return fullMatrix<double>(0, 0);
615 
616  const int szInt = (order - 1) * (order - 2) * (order - 1) / 2;
617  const int szComp = (order + 1) * (order + 1) * (order + 2) / 2;
618  const int szInc = szComp - szInt;
619 
620  fullMatrix<double> monomials = gmshGenerateMonomialsPrism(order, false);
621  fullMatrix<int> coordinates;
622  double2int(monomials, coordinates);
623 
624  std::map<mytuple, int> coord2idx;
625  for(int i = 0; i < szInc; ++i) {
626  int u = coordinates(i, 0);
627  int v = coordinates(i, 1);
628  int w = coordinates(i, 2);
629  coord2idx[make_mytuple(u, v, w)] = i;
630  }
631 
632  int &n = order;
633  fullMatrix<double> M(szInt, szInc, true);
634  for(int i = 0; i < szInt; ++i) {
635  int u = coordinates(szInc + i, 0);
636  int v = coordinates(szInc + i, 1);
637  int w = coordinates(szInc + i, 2);
638  double eta, mu;
639 
640  switch(dir) {
641  case 0:
642  mu = (double)u / (n - v);
643  M(i, coord2idx[make_mytuple(0, v, w)]) += 1 - mu;
644  M(i, coord2idx[make_mytuple(n - v, v, w)]) += mu;
645  break;
646  case 1:
647  mu = (double)v / (u + v);
648  M(i, coord2idx[make_mytuple(u + v, 0, w)]) += 1 - mu;
649  M(i, coord2idx[make_mytuple(0, u + v, w)]) += mu;
650  break;
651  case 2:
652  mu = (double)v / (n - u);
653  M(i, coord2idx[make_mytuple(u, 0, w)]) += 1 - mu;
654  M(i, coord2idx[make_mytuple(u, n - u, w)]) += mu;
655  break;
656  default:
657  case 3:
658  eta = (double)w / n;
659  M(i, coord2idx[make_mytuple(u, v, 0)]) += 1 - eta;
660  M(i, coord2idx[make_mytuple(u, v, n)]) += eta;
661  break;
662  }
663  }
664  return M;
665 }
666 
667 #undef make_mytuple
gmshGenerateInnerVertexPlacementQuadrangle
fullMatrix< double > gmshGenerateInnerVertexPlacementQuadrangle(int order)
Definition: InnerVertexPlacement.cpp:108
gmshGenerateInnerVertexPlacementPrismLinear
fullMatrix< double > gmshGenerateInnerVertexPlacementPrismLinear(int order, int dir)
Definition: InnerVertexPlacement.cpp:611
gmshGenerateInnerVertexPlacementQuadrangleLinear
fullMatrix< double > gmshGenerateInnerVertexPlacementQuadrangleLinear(int order)
Definition: InnerVertexPlacement.cpp:461
gmshGenerateInnerVertexPlacementTriangle
fullMatrix< double > gmshGenerateInnerVertexPlacementTriangle(int order)
Definition: InnerVertexPlacement.cpp:66
gmshGenerateMonomialsTriangle
fullMatrix< double > gmshGenerateMonomialsTriangle(int order, bool serendip)
Definition: pointsGenerators.cpp:182
gmshGenerateInnerVertexPlacementHexahedron
fullMatrix< double > gmshGenerateInnerVertexPlacementHexahedron(int order)
Definition: InnerVertexPlacement.cpp:213
getInnerVertexPlacement
fullMatrix< double > * getInnerVertexPlacement(int type, int order)
Definition: InnerVertexPlacement.cpp:33
c
static double c(int i, int j, fullMatrix< double > &CA, const std::vector< SPoint3 > &P, const std::vector< SPoint3 > &Q)
Definition: discreteFrechetDistance.cpp:15
gmshGenerateInnerVertexPlacementTetrahedronLinear
fullMatrix< double > gmshGenerateInnerVertexPlacementTetrahedronLinear(int order, int dir)
Definition: InnerVertexPlacement.cpp:493
gmshGenerateMonomialsPyramid
fullMatrix< double > gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints)
Definition: pointsGenerators.cpp:737
TYPE_TRI
#define TYPE_TRI
Definition: GmshDefines.h:66
InnerVertexPlacement.h
TYPE_PRI
#define TYPE_PRI
Definition: GmshDefines.h:70
gmshGenerateInnerVertexPlacementPrism
fullMatrix< double > gmshGenerateInnerVertexPlacementPrism(int order)
Definition: InnerVertexPlacement.cpp:275
fullMatrix< double >
gmshGenerateInnerVertexPlacementTriangleLinear
fullMatrix< double > gmshGenerateInnerVertexPlacementTriangleLinear(int order, int dir)
Definition: InnerVertexPlacement.cpp:411
GmshDefines.h
fullMatrix::size2
int size2() const
Definition: fullMatrix.h:275
TYPE_PYR
#define TYPE_PYR
Definition: GmshDefines.h:69
TYPE_QUA
#define TYPE_QUA
Definition: GmshDefines.h:67
fullMatrix::size1
int size1() const
Definition: fullMatrix.h:274
TYPE_HEX
#define TYPE_HEX
Definition: GmshDefines.h:71
TYPE_TET
#define TYPE_TET
Definition: GmshDefines.h:68
gmshGenerateInnerVertexPlacementTetrahedron
fullMatrix< double > gmshGenerateInnerVertexPlacementTetrahedron(int order)
Definition: InnerVertexPlacement.cpp:147
gmshGenerateInnerVertexPlacementPyramid
fullMatrix< double > gmshGenerateInnerVertexPlacementPyramid(int order)
Definition: InnerVertexPlacement.cpp:342
gmshGenerateInnerVertexPlacementHexahedronLinear
fullMatrix< double > gmshGenerateInnerVertexPlacementHexahedronLinear(int order, int dir)
Definition: InnerVertexPlacement.cpp:560
pointsGenerators.h
gmshGenerateMonomialsTetrahedron
fullMatrix< double > gmshGenerateMonomialsTetrahedron(int order, bool serendip)
Definition: pointsGenerators.cpp:317
fullMatrix::resize
bool resize(int r, int c, bool resetValue=true)
Definition: fullMatrix.h:307
gmshGenerateMonomialsPrism
fullMatrix< double > gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
Definition: pointsGenerators.cpp:402
gmshGenerateMonomialsQuadrangle
fullMatrix< double > gmshGenerateMonomialsQuadrangle(int order, bool forSerendipPoints)
Definition: pointsGenerators.cpp:222
gmshGenerateMonomialsHexahedron
fullMatrix< double > gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoints)
Definition: pointsGenerators.cpp:583