gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
adaptiveData.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 <math.h>
7 #include <list>
8 #include <set>
9 #include <algorithm>
10 #include "adaptiveData.h"
11 #include "PViewDataGModel.h"
12 #include "Plugin.h"
13 #include "OS.h"
14 #include "GmshDefines.h"
15 
16 //#define TIMER
17 
18 std::set<adaptiveVertex> adaptivePoint::allVertices;
19 std::set<adaptiveVertex> adaptiveLine::allVertices;
20 std::set<adaptiveVertex> adaptiveTriangle::allVertices;
21 std::set<adaptiveVertex> adaptiveQuadrangle::allVertices;
22 std::set<adaptiveVertex> adaptiveTetrahedron::allVertices;
23 std::set<adaptiveVertex> adaptiveHexahedron::allVertices;
24 std::set<adaptiveVertex> adaptivePrism::allVertices;
25 std::set<adaptiveVertex> adaptivePyramid::allVertices;
26 
27 std::list<adaptivePoint *> adaptivePoint::all;
28 std::list<adaptiveLine *> adaptiveLine::all;
29 std::list<adaptiveTriangle *> adaptiveTriangle::all;
30 std::list<adaptiveQuadrangle *> adaptiveQuadrangle::all;
31 std::list<adaptiveTetrahedron *> adaptiveTetrahedron::all;
32 std::list<adaptiveHexahedron *> adaptiveHexahedron::all;
33 std::list<adaptivePrism *> adaptivePrism::all;
34 std::list<adaptivePyramid *> adaptivePyramid::all;
35 
36 int adaptivePoint::numNodes = 1;
37 int adaptiveLine::numNodes = 2;
38 int adaptiveTriangle::numNodes = 3;
39 int adaptiveQuadrangle::numNodes = 4;
40 int adaptivePrism::numNodes = 6;
41 int adaptiveTetrahedron::numNodes = 4;
42 int adaptiveHexahedron::numNodes = 8;
43 int adaptivePyramid::numNodes = 5;
44 
45 int adaptivePoint::numEdges = 0;
46 int adaptiveLine::numEdges = 1;
47 int adaptiveTriangle::numEdges = 3;
48 int adaptiveQuadrangle::numEdges = 4;
49 int adaptivePrism::numEdges = 9;
50 int adaptiveTetrahedron::numEdges = 6;
51 int adaptiveHexahedron::numEdges = 12;
52 int adaptivePyramid::numEdges = 8;
53 
54 std::vector<vectInt> globalVTKData::vtkGlobalConnectivity;
55 std::vector<int> globalVTKData::vtkGlobalCellType;
56 std::vector<PCoords> globalVTKData::vtkGlobalCoords;
57 std::vector<PValues> globalVTKData::vtkGlobalValues;
58 
59 template <class T> static void cleanElement()
60 {
61  for(auto it = T::all.begin(); it != T::all.end(); ++it) delete *it;
62  T::all.clear();
63  T::allVertices.clear();
64 }
65 
66 static void computeShapeFunctions(fullMatrix<double> *coeffs,
67  fullMatrix<double> *eexps, double u, double v,
68  double w, fullVector<double> *sf,
69  fullVector<double> *tmp)
70 {
71  for(int i = 0; i < eexps->size1(); i++) {
72  (*tmp)(i) = pow(u, (*eexps)(i, 0));
73  if(eexps->size2() > 1) (*tmp)(i) *= pow(v, (*eexps)(i, 1));
74  if(eexps->size2() > 2) (*tmp)(i) *= pow(w, (*eexps)(i, 2));
75  }
76  coeffs->mult(*tmp, *sf);
77 }
78 
90 static void computeShapeFunctionsPyramid(fullMatrix<double> *coeffs,
91  fullMatrix<double> *eexps, double u,
92  double v, double w,
93  fullVector<double> *sf,
94  fullVector<double> *tmp)
95 {
96  double oneMinW = (w == 1) ? 1e-12 : 1 - w;
97  for(int l = 0; l < eexps->size1(); l++) {
98  int i = (*eexps)(l, 0);
99  int j = (*eexps)(l, 1);
100  int k = (*eexps)(l, 2);
101  int m = std::max(i, j);
102  (*tmp)(l) = pow(u, i);
103  (*tmp)(l) *= pow(v, j);
104  (*tmp)(l) *= pow(w, k);
105  (*tmp)(l) *= pow(oneMinW, m - i - j);
106  }
107  coeffs->mult(*tmp, *sf);
108 }
109 
110 adaptiveVertex *adaptiveVertex::add(double x, double y, double z,
111  std::set<adaptiveVertex> &allVertices)
112 {
113  adaptiveVertex p;
114  p.x = x;
115  p.y = y;
116  p.z = z;
117  auto it = allVertices.find(p);
118  if(it == allVertices.end()) {
119  allVertices.insert(p);
120  it = allVertices.find(p);
121  }
122  return (adaptiveVertex *)&(*it);
123 }
124 
125 void adaptivePoint::create(int maxlevel)
126 {
127  cleanElement<adaptivePoint>();
128  adaptiveVertex *p1 = adaptiveVertex::add(0, 0, 0, allVertices);
129  adaptivePoint *t = new adaptivePoint(p1);
130  recurCreate(t, maxlevel, 0);
131 }
132 
133 void adaptivePoint::recurCreate(adaptivePoint *e, int maxlevel, int level)
134 {
135  all.push_back(e);
136 }
137 
138 void adaptivePoint::error(double AVG, double tol)
139 {
140  adaptivePoint *e = *all.begin();
141  recurError(e, AVG, tol);
142 }
143 
144 void adaptivePoint::recurError(adaptivePoint *e, double AVG, double tol)
145 {
146  e->visible = true;
147 }
148 
149 void adaptiveLine::create(int maxlevel)
150 {
151  cleanElement<adaptiveLine>();
152  adaptiveVertex *p1 = adaptiveVertex::add(-1, 0, 0, allVertices);
153  adaptiveVertex *p2 = adaptiveVertex::add(1, 0, 0, allVertices);
154  adaptiveLine *t = new adaptiveLine(p1, p2);
155  recurCreate(t, maxlevel, 0);
156 }
157 
158 void adaptiveLine::recurCreate(adaptiveLine *e, int maxlevel, int level)
159 {
160  all.push_back(e);
161  if(level++ >= maxlevel) return;
162 
163  // p1 p12 p2
164  adaptiveVertex *p1 = e->p[0];
165  adaptiveVertex *p2 = e->p[1];
166  adaptiveVertex *p12 =
167  adaptiveVertex::add((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5,
168  (p1->z + p2->z) * 0.5, allVertices);
169  adaptiveLine *e1 = new adaptiveLine(p1, p12);
170  recurCreate(e1, maxlevel, level);
171  adaptiveLine *e2 = new adaptiveLine(p12, p2);
172  recurCreate(e2, maxlevel, level);
173  e->e[0] = e1;
174  e->e[1] = e2;
175 }
176 
177 void adaptiveLine::error(double AVG, double tol)
178 {
179  adaptiveLine *e = *all.begin();
180  recurError(e, AVG, tol);
181 }
182 
183 void adaptiveLine::recurError(adaptiveLine *e, double AVG, double tol)
184 {
185  if(!e->e[0])
186  e->visible = true;
187  else {
188  double vr;
189  if(!e->e[0]->e[0]) {
190  double v1 = e->e[0]->V();
191  double v2 = e->e[1]->V();
192  vr = (v1 + v2) / 2.;
193  double v = e->V();
194  if(fabs(v - vr) > AVG * tol) {
195  e->visible = false;
196  recurError(e->e[0], AVG, tol);
197  recurError(e->e[1], AVG, tol);
198  }
199  else
200  e->visible = true;
201  }
202  else {
203  double v11 = e->e[0]->e[0]->V();
204  double v12 = e->e[0]->e[1]->V();
205  double v21 = e->e[1]->e[0]->V();
206  double v22 = e->e[1]->e[1]->V();
207  double vr1 = (v11 + v12) / 2.;
208  double vr2 = (v21 + v22) / 2.;
209  vr = (vr1 + vr2) / 2.;
210  if(fabs(e->e[0]->V() - vr1) > AVG * tol ||
211  fabs(e->e[1]->V() - vr2) > AVG * tol ||
212  fabs(e->V() - vr) > AVG * tol) {
213  e->visible = false;
214  recurError(e->e[0], AVG, tol);
215  recurError(e->e[1], AVG, tol);
216  }
217  else
218  e->visible = true;
219  }
220  }
221 }
222 
223 void adaptiveTriangle::create(int maxlevel)
224 {
225  cleanElement<adaptiveTriangle>();
226  adaptiveVertex *p1 = adaptiveVertex::add(0, 0, 0, allVertices);
227  adaptiveVertex *p2 = adaptiveVertex::add(0, 1, 0, allVertices);
228  adaptiveVertex *p3 = adaptiveVertex::add(1, 0, 0, allVertices);
229  adaptiveTriangle *t = new adaptiveTriangle(p1, p2, p3);
230  recurCreate(t, maxlevel, 0);
231 }
232 
233 void adaptiveTriangle::recurCreate(adaptiveTriangle *t, int maxlevel, int level)
234 {
235  all.push_back(t);
236  if(level++ >= maxlevel) return;
237 
238  // p3
239  // p13 p23
240  // p1 p12 p2
241  adaptiveVertex *p1 = t->p[0];
242  adaptiveVertex *p2 = t->p[1];
243  adaptiveVertex *p3 = t->p[2];
244  adaptiveVertex *p12 =
245  adaptiveVertex::add((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5,
246  (p1->z + p2->z) * 0.5, allVertices);
247  adaptiveVertex *p13 =
248  adaptiveVertex::add((p1->x + p3->x) * 0.5, (p1->y + p3->y) * 0.5,
249  (p1->z + p3->z) * 0.5, allVertices);
250  adaptiveVertex *p23 =
251  adaptiveVertex::add((p3->x + p2->x) * 0.5, (p3->y + p2->y) * 0.5,
252  (p3->z + p2->z) * 0.5, allVertices);
253  adaptiveTriangle *t1 = new adaptiveTriangle(p1, p12, p13);
254  recurCreate(t1, maxlevel, level);
255  adaptiveTriangle *t2 = new adaptiveTriangle(p2, p23, p12);
256  recurCreate(t2, maxlevel, level);
257  adaptiveTriangle *t3 = new adaptiveTriangle(p3, p13, p23);
258  recurCreate(t3, maxlevel, level);
259  adaptiveTriangle *t4 = new adaptiveTriangle(p12, p23, p13);
260  recurCreate(t4, maxlevel, level);
261  t->e[0] = t1;
262  t->e[1] = t2;
263  t->e[2] = t3;
264  t->e[3] = t4;
265 }
266 
267 void adaptiveTriangle::error(double AVG, double tol)
268 {
269  adaptiveTriangle *t = *all.begin();
270  recurError(t, AVG, tol);
271 }
272 
273 void adaptiveTriangle::recurError(adaptiveTriangle *t, double AVG, double tol)
274 {
275  if(!t->e[0])
276  t->visible = true;
277  else {
278  double vr;
279  if(!t->e[0]->e[0]) {
280  double v1 = t->e[0]->V();
281  double v2 = t->e[1]->V();
282  double v3 = t->e[2]->V();
283  double v4 = t->e[3]->V();
284  vr = (2 * v1 + 2 * v2 + 2 * v3 + v4) / 7.;
285  double v = t->V();
286  if(fabs(v - vr) > AVG * tol) {
287  t->visible = false;
288  recurError(t->e[0], AVG, tol);
289  recurError(t->e[1], AVG, tol);
290  recurError(t->e[2], AVG, tol);
291  recurError(t->e[3], AVG, tol);
292  }
293  else
294  t->visible = true;
295  }
296  else {
297  double v11 = t->e[0]->e[0]->V();
298  double v12 = t->e[0]->e[1]->V();
299  double v13 = t->e[0]->e[2]->V();
300  double v14 = t->e[0]->e[3]->V();
301  double v21 = t->e[1]->e[0]->V();
302  double v22 = t->e[1]->e[1]->V();
303  double v23 = t->e[1]->e[2]->V();
304  double v24 = t->e[1]->e[3]->V();
305  double v31 = t->e[2]->e[0]->V();
306  double v32 = t->e[2]->e[1]->V();
307  double v33 = t->e[2]->e[2]->V();
308  double v34 = t->e[2]->e[3]->V();
309  double v41 = t->e[3]->e[0]->V();
310  double v42 = t->e[3]->e[1]->V();
311  double v43 = t->e[3]->e[2]->V();
312  double v44 = t->e[3]->e[3]->V();
313  double vr1 = (2 * v11 + 2 * v12 + 2 * v13 + v14) / 7.;
314  double vr2 = (2 * v21 + 2 * v22 + 2 * v23 + v24) / 7.;
315  double vr3 = (2 * v31 + 2 * v32 + 2 * v33 + v34) / 7.;
316  double vr4 = (2 * v41 + 2 * v42 + 2 * v43 + v44) / 7.;
317  vr = (2 * vr1 + 2 * vr2 + 2 * vr3 + vr4) / 7.;
318  if(fabs(t->e[0]->V() - vr1) > AVG * tol ||
319  fabs(t->e[1]->V() - vr2) > AVG * tol ||
320  fabs(t->e[2]->V() - vr3) > AVG * tol ||
321  fabs(t->e[3]->V() - vr4) > AVG * tol ||
322  fabs(t->V() - vr) > AVG * tol) {
323  t->visible = false;
324  recurError(t->e[0], AVG, tol);
325  recurError(t->e[1], AVG, tol);
326  recurError(t->e[2], AVG, tol);
327  recurError(t->e[3], AVG, tol);
328  }
329  else
330  t->visible = true;
331  }
332  }
333 }
334 
335 void adaptiveQuadrangle::create(int maxlevel)
336 {
337  cleanElement<adaptiveQuadrangle>();
338  adaptiveVertex *p1 = adaptiveVertex::add(-1, -1, 0, allVertices);
339  adaptiveVertex *p2 = adaptiveVertex::add(1, -1, 0, allVertices);
340  adaptiveVertex *p3 = adaptiveVertex::add(1, 1, 0, allVertices);
341  adaptiveVertex *p4 = adaptiveVertex::add(-1, 1, 0, allVertices);
342  adaptiveQuadrangle *q = new adaptiveQuadrangle(p1, p2, p3, p4);
343  recurCreate(q, maxlevel, 0);
344 }
345 
346 void adaptiveQuadrangle::recurCreate(adaptiveQuadrangle *q, int maxlevel,
347  int level)
348 {
349  all.push_back(q);
350  if(level++ >= maxlevel) return;
351 
352  // p4 p34 p3
353  // p14 pc p23
354  // p1 p12 p2
355  adaptiveVertex *p1 = q->p[0];
356  adaptiveVertex *p2 = q->p[1];
357  adaptiveVertex *p3 = q->p[2];
358  adaptiveVertex *p4 = q->p[3];
359  adaptiveVertex *p12 =
360  adaptiveVertex::add((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5,
361  (p1->z + p2->z) * 0.5, allVertices);
362  adaptiveVertex *p23 =
363  adaptiveVertex::add((p2->x + p3->x) * 0.5, (p2->y + p3->y) * 0.5,
364  (p2->z + p3->z) * 0.5, allVertices);
365  adaptiveVertex *p34 =
366  adaptiveVertex::add((p3->x + p4->x) * 0.5, (p3->y + p4->y) * 0.5,
367  (p3->z + p4->z) * 0.5, allVertices);
368  adaptiveVertex *p14 =
369  adaptiveVertex::add((p1->x + p4->x) * 0.5, (p1->y + p4->y) * 0.5,
370  (p1->z + p4->z) * 0.5, allVertices);
371  adaptiveVertex *pc =
372  adaptiveVertex::add((p1->x + p2->x + p3->x + p4->x) * 0.25,
373  (p1->y + p2->y + p3->y + p4->y) * 0.25,
374  (p1->z + p2->z + p3->z + p4->z) * 0.25, allVertices);
375  adaptiveQuadrangle *q1 = new adaptiveQuadrangle(p1, p12, pc, p14);
376  recurCreate(q1, maxlevel, level);
377  adaptiveQuadrangle *q2 = new adaptiveQuadrangle(p2, p23, pc, p12);
378  recurCreate(q2, maxlevel, level);
379  adaptiveQuadrangle *q3 = new adaptiveQuadrangle(p3, p34, pc, p23);
380  recurCreate(q3, maxlevel, level);
381  adaptiveQuadrangle *q4 = new adaptiveQuadrangle(p4, p14, pc, p34);
382  recurCreate(q4, maxlevel, level);
383  q->e[0] = q1;
384  q->e[1] = q2;
385  q->e[2] = q3;
386  q->e[3] = q4;
387 }
388 
389 void adaptiveQuadrangle::error(double AVG, double tol)
390 {
391  adaptiveQuadrangle *q = *all.begin();
392  recurError(q, AVG, tol);
393 }
394 
395 void adaptiveQuadrangle::recurError(adaptiveQuadrangle *q, double AVG,
396  double tol)
397 {
398  if(!q->e[0])
399  q->visible = true;
400  else {
401  double vr;
402  double vd = (q->p[0]->val + q->p[2]->val) / 2.;
403  if(!q->e[0]->e[0]) {
404  double v1 = q->e[0]->V();
405  double v2 = q->e[1]->V();
406  double v3 = q->e[2]->V();
407  double v4 = q->e[3]->V();
408  vr = (v1 + v2 + v3 + v4) / 4.;
409  double v = q->V();
410  if(fabs(v - vr) > AVG * tol || fabs(vd - vr) > AVG * tol) {
411  q->visible = false;
412  recurError(q->e[0], AVG, tol);
413  recurError(q->e[1], AVG, tol);
414  recurError(q->e[2], AVG, tol);
415  recurError(q->e[3], AVG, tol);
416  }
417  else
418  q->visible = true;
419  }
420  else {
421  double v11 = q->e[0]->e[0]->V();
422  double v12 = q->e[0]->e[1]->V();
423  double v13 = q->e[0]->e[2]->V();
424  double v14 = q->e[0]->e[3]->V();
425  double v21 = q->e[1]->e[0]->V();
426  double v22 = q->e[1]->e[1]->V();
427  double v23 = q->e[1]->e[2]->V();
428  double v24 = q->e[1]->e[3]->V();
429  double v31 = q->e[2]->e[0]->V();
430  double v32 = q->e[2]->e[1]->V();
431  double v33 = q->e[2]->e[2]->V();
432  double v34 = q->e[2]->e[3]->V();
433  double v41 = q->e[3]->e[0]->V();
434  double v42 = q->e[3]->e[1]->V();
435  double v43 = q->e[3]->e[2]->V();
436  double v44 = q->e[3]->e[3]->V();
437  double vr1 = (v11 + v12 + v13 + v14) / 4.;
438  double vr2 = (v21 + v22 + v23 + v24) / 4.;
439  double vr3 = (v31 + v32 + v33 + v34) / 4.;
440  double vr4 = (v41 + v42 + v43 + v44) / 4.;
441  vr = (vr1 + vr2 + vr3 + vr4) / 4.;
442  if(fabs(q->e[0]->V() - vr1) > AVG * tol ||
443  fabs(q->e[1]->V() - vr2) > AVG * tol ||
444  fabs(q->e[2]->V() - vr3) > AVG * tol ||
445  fabs(q->e[3]->V() - vr4) > AVG * tol ||
446  fabs(q->V() - vr) > AVG * tol || fabs(vd - vr) > AVG * tol) {
447  q->visible = false;
448  recurError(q->e[0], AVG, tol);
449  recurError(q->e[1], AVG, tol);
450  recurError(q->e[2], AVG, tol);
451  recurError(q->e[3], AVG, tol);
452  }
453  else
454  q->visible = true;
455  }
456  }
457 }
458 
459 void adaptiveTetrahedron::create(int maxlevel)
460 {
461  cleanElement<adaptiveTetrahedron>();
462  adaptiveVertex *p1 = adaptiveVertex::add(0, 0, 0, allVertices);
463  adaptiveVertex *p2 = adaptiveVertex::add(0, 1, 0, allVertices);
464  adaptiveVertex *p3 = adaptiveVertex::add(1, 0, 0, allVertices);
465  adaptiveVertex *p4 = adaptiveVertex::add(0, 0, 1, allVertices);
466  adaptiveTetrahedron *t = new adaptiveTetrahedron(p1, p2, p3, p4);
467  recurCreate(t, maxlevel, 0);
468 }
469 
470 void adaptiveTetrahedron::recurCreate(adaptiveTetrahedron *t, int maxlevel,
471  int level)
472 {
473  all.push_back(t);
474  if(level++ >= maxlevel) return;
475 
476  adaptiveVertex *p0 = t->p[0];
477  adaptiveVertex *p1 = t->p[1];
478  adaptiveVertex *p2 = t->p[2];
479  adaptiveVertex *p3 = t->p[3];
480  adaptiveVertex *pe0 =
481  adaptiveVertex::add((p0->x + p1->x) * 0.5, (p0->y + p1->y) * 0.5,
482  (p0->z + p1->z) * 0.5, allVertices);
483  adaptiveVertex *pe1 =
484  adaptiveVertex::add((p0->x + p2->x) * 0.5, (p0->y + p2->y) * 0.5,
485  (p0->z + p2->z) * 0.5, allVertices);
486  adaptiveVertex *pe2 =
487  adaptiveVertex::add((p0->x + p3->x) * 0.5, (p0->y + p3->y) * 0.5,
488  (p0->z + p3->z) * 0.5, allVertices);
489  adaptiveVertex *pe3 =
490  adaptiveVertex::add((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5,
491  (p1->z + p2->z) * 0.5, allVertices);
492  adaptiveVertex *pe4 =
493  adaptiveVertex::add((p1->x + p3->x) * 0.5, (p1->y + p3->y) * 0.5,
494  (p1->z + p3->z) * 0.5, allVertices);
495  adaptiveVertex *pe5 =
496  adaptiveVertex::add((p2->x + p3->x) * 0.5, (p2->y + p3->y) * 0.5,
497  (p2->z + p3->z) * 0.5, allVertices);
498  adaptiveTetrahedron *t1 = new adaptiveTetrahedron(p0, pe0, pe1, pe2);
499  recurCreate(t1, maxlevel, level);
500  adaptiveTetrahedron *t2 = new adaptiveTetrahedron(pe0, p1, pe3, pe4);
501  recurCreate(t2, maxlevel, level);
502  adaptiveTetrahedron *t3 = new adaptiveTetrahedron(pe1, pe3, p2, pe5);
503  recurCreate(t3, maxlevel, level);
504  adaptiveTetrahedron *t4 = new adaptiveTetrahedron(pe2, pe4, pe5, p3);
505  recurCreate(t4, maxlevel, level);
506  adaptiveTetrahedron *t5 = new adaptiveTetrahedron(pe3, pe5, pe2, pe4);
507  recurCreate(t5, maxlevel, level);
508  adaptiveTetrahedron *t6 = new adaptiveTetrahedron(pe3, pe2, pe0, pe4);
509  recurCreate(t6, maxlevel, level);
510  adaptiveTetrahedron *t7 = new adaptiveTetrahedron(pe2, pe5, pe3, pe1);
511  recurCreate(t7, maxlevel, level);
512  adaptiveTetrahedron *t8 = new adaptiveTetrahedron(pe0, pe2, pe3, pe1);
513  recurCreate(t8, maxlevel, level);
514  t->e[0] = t1;
515  t->e[1] = t2;
516  t->e[2] = t3;
517  t->e[3] = t4;
518  t->e[4] = t5;
519  t->e[5] = t6;
520  t->e[6] = t7;
521  t->e[7] = t8;
522 }
523 
524 void adaptiveTetrahedron::error(double AVG, double tol)
525 {
526  adaptiveTetrahedron *t = *all.begin();
527  recurError(t, AVG, tol);
528 }
529 
530 void adaptiveTetrahedron::recurError(adaptiveTetrahedron *t, double AVG,
531  double tol)
532 {
533  if(!t->e[0])
534  t->visible = true;
535  else {
536  const double v1 = t->e[0]->V();
537  const double v2 = t->e[1]->V();
538  const double v3 = t->e[2]->V();
539  const double v4 = t->e[3]->V();
540  const double v5 = t->e[4]->V();
541  const double v6 = t->e[5]->V();
542  const double v7 = t->e[6]->V();
543  const double v8 = t->e[7]->V();
544  const double vr = (v1 + v2 + v3 + v4 + v5 + v6 + v7 + v8) * .125;
545  const double v = t->V();
546  if(!t->e[0]->e[0]) {
547  if(fabs(v - vr) > AVG * tol) {
548  t->visible = false;
549  recurError(t->e[0], AVG, tol);
550  recurError(t->e[1], AVG, tol);
551  recurError(t->e[2], AVG, tol);
552  recurError(t->e[3], AVG, tol);
553  recurError(t->e[4], AVG, tol);
554  recurError(t->e[5], AVG, tol);
555  recurError(t->e[6], AVG, tol);
556  recurError(t->e[7], AVG, tol);
557  }
558  else
559  t->visible = true;
560  }
561  else {
562  double vi[8][8];
563  for(int k = 0; k < 8; k++)
564  for(int l = 0; l < 8; l++) vi[k][l] = t->e[k]->e[l]->V();
565  double vri[8];
566  for(int k = 0; k < 8; k++) {
567  vri[k] = 0.0;
568  for(int l = 0; l < 8; l++) { vri[k] += vi[k][l]; }
569  vri[k] /= 8.0;
570  }
571  if(fabs(t->e[0]->V() - vri[0]) > AVG * tol ||
572  fabs(t->e[1]->V() - vri[1]) > AVG * tol ||
573  fabs(t->e[2]->V() - vri[2]) > AVG * tol ||
574  fabs(t->e[3]->V() - vri[3]) > AVG * tol ||
575  fabs(t->e[4]->V() - vri[4]) > AVG * tol ||
576  fabs(t->e[5]->V() - vri[5]) > AVG * tol ||
577  fabs(t->e[6]->V() - vri[6]) > AVG * tol ||
578  fabs(t->e[7]->V() - vri[7]) > AVG * tol || fabs(v - vr) > AVG * tol) {
579  t->visible = false;
580  recurError(t->e[0], AVG, tol);
581  recurError(t->e[1], AVG, tol);
582  recurError(t->e[2], AVG, tol);
583  recurError(t->e[3], AVG, tol);
584  recurError(t->e[4], AVG, tol);
585  recurError(t->e[5], AVG, tol);
586  recurError(t->e[6], AVG, tol);
587  recurError(t->e[7], AVG, tol);
588  }
589  else
590  t->visible = true;
591  }
592  }
593 }
594 
595 void adaptiveHexahedron::create(int maxlevel)
596 {
597  cleanElement<adaptiveHexahedron>();
598  adaptiveVertex *p1 = adaptiveVertex::add(-1, -1, -1, allVertices);
599  adaptiveVertex *p2 = adaptiveVertex::add(-1, 1, -1, allVertices);
600  adaptiveVertex *p3 = adaptiveVertex::add(1, 1, -1, allVertices);
601  adaptiveVertex *p4 = adaptiveVertex::add(1, -1, -1, allVertices);
602  adaptiveVertex *p11 = adaptiveVertex::add(-1, -1, 1, allVertices);
603  adaptiveVertex *p21 = adaptiveVertex::add(-1, 1, 1, allVertices);
604  adaptiveVertex *p31 = adaptiveVertex::add(1, 1, 1, allVertices);
605  adaptiveVertex *p41 = adaptiveVertex::add(1, -1, 1, allVertices);
606  adaptiveHexahedron *h =
607  new adaptiveHexahedron(p1, p2, p3, p4, p11, p21, p31, p41);
608  recurCreate(h, maxlevel, 0);
609 }
610 
611 void adaptiveHexahedron::recurCreate(adaptiveHexahedron *h, int maxlevel,
612  int level)
613 {
614  all.push_back(h);
615  if(level++ >= maxlevel) return;
616 
617  adaptiveVertex *p0 = h->p[0];
618  adaptiveVertex *p1 = h->p[1];
619  adaptiveVertex *p2 = h->p[2];
620  adaptiveVertex *p3 = h->p[3];
621  adaptiveVertex *p4 = h->p[4];
622  adaptiveVertex *p5 = h->p[5];
623  adaptiveVertex *p6 = h->p[6];
624  adaptiveVertex *p7 = h->p[7];
625  adaptiveVertex *p01 =
626  adaptiveVertex::add((p0->x + p1->x) * 0.5, (p0->y + p1->y) * 0.5,
627  (p0->z + p1->z) * 0.5, allVertices);
628  adaptiveVertex *p12 =
629  adaptiveVertex::add((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5,
630  (p1->z + p2->z) * 0.5, allVertices);
631  adaptiveVertex *p23 =
632  adaptiveVertex::add((p2->x + p3->x) * 0.5, (p2->y + p3->y) * 0.5,
633  (p2->z + p3->z) * 0.5, allVertices);
634  adaptiveVertex *p03 =
635  adaptiveVertex::add((p3->x + p0->x) * 0.5, (p3->y + p0->y) * 0.5,
636  (p3->z + p0->z) * 0.5, allVertices);
637  adaptiveVertex *p45 =
638  adaptiveVertex::add((p4->x + p5->x) * 0.5, (p4->y + p5->y) * 0.5,
639  (p4->z + p5->z) * 0.5, allVertices);
640  adaptiveVertex *p56 =
641  adaptiveVertex::add((p5->x + p6->x) * 0.5, (p5->y + p6->y) * 0.5,
642  (p5->z + p6->z) * 0.5, allVertices);
643  adaptiveVertex *p67 =
644  adaptiveVertex::add((p6->x + p7->x) * 0.5, (p6->y + p7->y) * 0.5,
645  (p6->z + p7->z) * 0.5, allVertices);
646  adaptiveVertex *p47 =
647  adaptiveVertex::add((p7->x + p4->x) * 0.5, (p7->y + p4->y) * 0.5,
648  (p7->z + p4->z) * 0.5, allVertices);
649  adaptiveVertex *p04 =
650  adaptiveVertex::add((p4->x + p0->x) * 0.5, (p4->y + p0->y) * 0.5,
651  (p4->z + p0->z) * 0.5, allVertices);
652  adaptiveVertex *p15 =
653  adaptiveVertex::add((p5->x + p1->x) * 0.5, (p5->y + p1->y) * 0.5,
654  (p5->z + p1->z) * 0.5, allVertices);
655  adaptiveVertex *p26 =
656  adaptiveVertex::add((p6->x + p2->x) * 0.5, (p6->y + p2->y) * 0.5,
657  (p6->z + p2->z) * 0.5, allVertices);
658  adaptiveVertex *p37 =
659  adaptiveVertex::add((p7->x + p3->x) * 0.5, (p7->y + p3->y) * 0.5,
660  (p7->z + p3->z) * 0.5, allVertices);
661  adaptiveVertex *p0145 =
662  adaptiveVertex::add((p45->x + p01->x) * 0.5, (p45->y + p01->y) * 0.5,
663  (p45->z + p01->z) * 0.5, allVertices);
664  adaptiveVertex *p1256 =
665  adaptiveVertex::add((p12->x + p56->x) * 0.5, (p12->y + p56->y) * 0.5,
666  (p12->z + p56->z) * 0.5, allVertices);
667  adaptiveVertex *p2367 =
668  adaptiveVertex::add((p23->x + p67->x) * 0.5, (p23->y + p67->y) * 0.5,
669  (p23->z + p67->z) * 0.5, allVertices);
670  adaptiveVertex *p0347 =
671  adaptiveVertex::add((p03->x + p47->x) * 0.5, (p03->y + p47->y) * 0.5,
672  (p03->z + p47->z) * 0.5, allVertices);
673  adaptiveVertex *p4756 =
674  adaptiveVertex::add((p47->x + p56->x) * 0.5, (p47->y + p56->y) * 0.5,
675  (p47->z + p56->z) * 0.5, allVertices);
676  adaptiveVertex *p0312 =
677  adaptiveVertex::add((p03->x + p12->x) * 0.5, (p03->y + p12->y) * 0.5,
678  (p03->z + p12->z) * 0.5, allVertices);
679  adaptiveVertex *pc = adaptiveVertex::add(
680  (p0->x + p1->x + p2->x + p3->x + p4->x + p5->x + p6->x + p7->x) * 0.125,
681  (p0->y + p1->y + p2->y + p3->y + p4->y + p5->y + p6->y + p7->y) * 0.125,
682  (p0->z + p1->z + p2->z + p3->z + p4->z + p5->z + p6->z + p7->z) * 0.125,
683  allVertices);
684 
685  adaptiveHexahedron *h1 =
686  new adaptiveHexahedron(p0, p01, p0312, p03, p04, p0145, pc, p0347); // p0
687  recurCreate(h1, maxlevel, level);
688  adaptiveHexahedron *h2 =
689  new adaptiveHexahedron(p01, p0145, p15, p1, p0312, pc, p1256, p12); // p1
690  recurCreate(h2, maxlevel, level);
691  adaptiveHexahedron *h3 =
692  new adaptiveHexahedron(p04, p4, p45, p0145, p0347, p47, p4756, pc); // p4
693  recurCreate(h3, maxlevel, level);
694  adaptiveHexahedron *h4 =
695  new adaptiveHexahedron(p0145, p45, p5, p15, pc, p4756, p56, p1256); // p5
696  recurCreate(h4, maxlevel, level);
697  adaptiveHexahedron *h5 =
698  new adaptiveHexahedron(p0347, p47, p4756, pc, p37, p7, p67, p2367); // p7
699  recurCreate(h5, maxlevel, level);
700  adaptiveHexahedron *h6 =
701  new adaptiveHexahedron(pc, p4756, p56, p1256, p2367, p67, p6, p26); // p6
702  recurCreate(h6, maxlevel, level);
703  adaptiveHexahedron *h7 =
704  new adaptiveHexahedron(p03, p0347, pc, p0312, p3, p37, p2367, p23); // p3
705  recurCreate(h7, maxlevel, level);
706  adaptiveHexahedron *h8 =
707  new adaptiveHexahedron(p0312, pc, p1256, p12, p23, p2367, p26, p2); // p2
708  recurCreate(h8, maxlevel, level);
709  h->e[0] = h1;
710  h->e[1] = h2;
711  h->e[2] = h3;
712  h->e[3] = h4;
713  h->e[4] = h5;
714  h->e[5] = h6;
715  h->e[6] = h7;
716  h->e[7] = h8;
717 }
718 
719 void adaptiveHexahedron::error(double AVG, double tol)
720 {
721  adaptiveHexahedron *h = *all.begin();
722  recurError(h, AVG, tol);
723 }
724 
725 void adaptiveHexahedron::recurError(adaptiveHexahedron *h, double AVG,
726  double tol)
727 {
728  if(!h->e[0])
729  h->visible = true;
730  else {
731  const double v1 = h->e[0]->V();
732  const double v2 = h->e[1]->V();
733  const double v3 = h->e[2]->V();
734  const double v4 = h->e[3]->V();
735  const double v5 = h->e[4]->V();
736  const double v6 = h->e[5]->V();
737  const double v7 = h->e[6]->V();
738  const double v8 = h->e[7]->V();
739  const double vr = (v1 + v2 + v3 + v4 + v5 + v6 + v7 + v8) * .125;
740  const double v = h->V();
741  // we use diagonal 1-7 because it's the one used for drawing
742  const double vd = (h->p[1]->val + h->p[7]->val) / 2.;
743  if(!h->e[0]->e[0]) {
744  if(fabs(v - vr) > AVG * tol || fabs(vd - vr) > AVG * tol) {
745  h->visible = false;
746  recurError(h->e[0], AVG, tol);
747  recurError(h->e[1], AVG, tol);
748  recurError(h->e[2], AVG, tol);
749  recurError(h->e[3], AVG, tol);
750  recurError(h->e[4], AVG, tol);
751  recurError(h->e[5], AVG, tol);
752  recurError(h->e[6], AVG, tol);
753  recurError(h->e[7], AVG, tol);
754  }
755  else
756  h->visible = true;
757  }
758  else {
759  double vii[8][8];
760  for(int i = 0; i < 8; i++)
761  for(int j = 0; j < 8; j++) vii[i][j] = h->e[i]->e[j]->V();
762  double vri[8];
763  for(int k = 0; k < 8; k++) {
764  vri[k] = 0.0;
765  for(int l = 0; l < 8; l++) { vri[k] += vii[k][l]; }
766  vri[k] /= 8.0;
767  }
768  if(fabs(h->e[0]->V() - vri[0]) > AVG * tol ||
769  fabs(h->e[1]->V() - vri[1]) > AVG * tol ||
770  fabs(h->e[2]->V() - vri[2]) > AVG * tol ||
771  fabs(h->e[3]->V() - vri[3]) > AVG * tol ||
772  fabs(h->e[4]->V() - vri[4]) > AVG * tol ||
773  fabs(h->e[5]->V() - vri[5]) > AVG * tol ||
774  fabs(h->e[6]->V() - vri[6]) > AVG * tol ||
775  fabs(h->e[7]->V() - vri[7]) > AVG * tol || fabs(v - vr) > AVG * tol ||
776  fabs(vd - vr) > AVG * tol) {
777  h->visible = false;
778  recurError(h->e[0], AVG, tol);
779  recurError(h->e[1], AVG, tol);
780  recurError(h->e[2], AVG, tol);
781  recurError(h->e[3], AVG, tol);
782  recurError(h->e[4], AVG, tol);
783  recurError(h->e[5], AVG, tol);
784  recurError(h->e[6], AVG, tol);
785  recurError(h->e[7], AVG, tol);
786  }
787  else
788  h->visible = true;
789  }
790  }
791 }
792 
793 void adaptivePrism::create(int maxlevel)
794 {
795  cleanElement<adaptivePrism>();
796  adaptiveVertex *p1 = adaptiveVertex::add(0, 0, -1, allVertices);
797  adaptiveVertex *p2 = adaptiveVertex::add(1, 0, -1, allVertices);
798  adaptiveVertex *p3 = adaptiveVertex::add(0, 1, -1, allVertices);
799  adaptiveVertex *p4 = adaptiveVertex::add(0, 0, 1, allVertices);
800  adaptiveVertex *p5 = adaptiveVertex::add(1, 0, 1, allVertices);
801  adaptiveVertex *p6 = adaptiveVertex::add(0, 1, 1, allVertices);
802  adaptivePrism *p = new adaptivePrism(p1, p2, p3, p4, p5, p6);
803  recurCreate(p, maxlevel, 0);
804 }
805 
806 void adaptivePrism::recurCreate(adaptivePrism *p, int maxlevel, int level)
807 {
808  all.push_back(p);
809  if(level++ >= maxlevel) return;
810 
811  // p4 p34 p3
812  // p14 pc p23
813  // p1 p12 p2
814  adaptiveVertex *p1 = p->p[0];
815  adaptiveVertex *p2 = p->p[1];
816  adaptiveVertex *p3 = p->p[2];
817  adaptiveVertex *p4 = p->p[3];
818  adaptiveVertex *p5 = p->p[4];
819  adaptiveVertex *p6 = p->p[5];
820  adaptiveVertex *p14 =
821  adaptiveVertex::add((p1->x + p4->x) * 0.5, (p1->y + p4->y) * 0.5,
822  (p1->z + p4->z) * 0.5, allVertices);
823  adaptiveVertex *p25 =
824  adaptiveVertex::add((p2->x + p5->x) * 0.5, (p2->y + p5->y) * 0.5,
825  (p2->z + p5->z) * 0.5, allVertices);
826  adaptiveVertex *p36 =
827  adaptiveVertex::add((p3->x + p6->x) * 0.5, (p3->y + p6->y) * 0.5,
828  (p3->z + p6->z) * 0.5, allVertices);
829  adaptiveVertex *p12 =
830  adaptiveVertex::add((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5,
831  (p1->z + p2->z) * 0.5, allVertices);
832  adaptiveVertex *p23 =
833  adaptiveVertex::add((p2->x + p3->x) * 0.5, (p2->y + p3->y) * 0.5,
834  (p2->z + p3->z) * 0.5, allVertices);
835  adaptiveVertex *p31 =
836  adaptiveVertex::add((p3->x + p1->x) * 0.5, (p3->y + p1->y) * 0.5,
837  (p3->z + p1->z) * 0.5, allVertices);
838  adaptiveVertex *p1425 =
839  adaptiveVertex::add((p14->x + p25->x) * 0.5, (p14->y + p25->y) * 0.5,
840  (p14->z + p25->z) * 0.5, allVertices);
841  adaptiveVertex *p2536 =
842  adaptiveVertex::add((p25->x + p36->x) * 0.5, (p25->y + p36->y) * 0.5,
843  (p25->z + p36->z) * 0.5, allVertices);
844  adaptiveVertex *p3614 =
845  adaptiveVertex::add((p36->x + p14->x) * 0.5, (p36->y + p14->y) * 0.5,
846  (p36->z + p14->z) * 0.5, allVertices);
847  adaptiveVertex *p45 =
848  adaptiveVertex::add((p4->x + p5->x) * 0.5, (p4->y + p5->y) * 0.5,
849  (p4->z + p5->z) * 0.5, allVertices);
850  adaptiveVertex *p56 =
851  adaptiveVertex::add((p5->x + p6->x) * 0.5, (p5->y + p6->y) * 0.5,
852  (p5->z + p6->z) * 0.5, allVertices);
853  adaptiveVertex *p64 =
854  adaptiveVertex::add((p6->x + p4->x) * 0.5, (p6->y + p4->y) * 0.5,
855  (p6->z + p4->z) * 0.5, allVertices);
856  p->e[0] = new adaptivePrism(p1, p12, p31, p14, p1425, p3614);
857  recurCreate(p->e[0], maxlevel, level);
858  p->e[1] = new adaptivePrism(p2, p23, p12, p25, p2536, p1425);
859  recurCreate(p->e[1], maxlevel, level);
860  p->e[2] = new adaptivePrism(p3, p31, p23, p36, p3614, p2536);
861  recurCreate(p->e[2], maxlevel, level);
862  p->e[3] = new adaptivePrism(p12, p23, p31, p1425, p2536, p3614);
863  recurCreate(p->e[3], maxlevel, level);
864  p->e[4] = new adaptivePrism(p14, p1425, p3614, p4, p45, p64);
865  recurCreate(p->e[4], maxlevel, level);
866  p->e[5] = new adaptivePrism(p25, p2536, p1425, p5, p56, p45);
867  recurCreate(p->e[5], maxlevel, level);
868  p->e[6] = new adaptivePrism(p36, p3614, p2536, p6, p64, p56);
869  recurCreate(p->e[6], maxlevel, level);
870  p->e[7] = new adaptivePrism(p1425, p2536, p3614, p45, p56, p64);
871  recurCreate(p->e[7], maxlevel, level);
872 }
873 
874 void adaptivePrism::error(double AVG, double tol)
875 {
876  adaptivePrism *p = *all.begin();
877  recurError(p, AVG, tol);
878 }
879 
880 void adaptivePrism::recurError(adaptivePrism *p, double AVG, double tol)
881 {
882  if(!p->e[0])
883  p->visible = true;
884  else {
885  double vi[8];
886  for(int i = 0; i < 8; i++) vi[i] = p->e[i]->V();
887  const double vr =
888  (vi[0] + vi[1] + vi[2] + vi[3] / 2 + vi[4] + vi[5] + vi[6] + vi[7] / 2) /
889  7;
890  const double v = p->V();
891  if(!p->e[0]->e[0]) {
892  if(fabs(v - vr) > AVG * tol) {
893  p->visible = false;
894  recurError(p->e[0], AVG, tol);
895  recurError(p->e[1], AVG, tol);
896  recurError(p->e[2], AVG, tol);
897  recurError(p->e[3], AVG, tol);
898  recurError(p->e[4], AVG, tol);
899  recurError(p->e[5], AVG, tol);
900  recurError(p->e[6], AVG, tol);
901  recurError(p->e[7], AVG, tol);
902  }
903  else
904  p->visible = true;
905  }
906  else {
907  bool err = false;
908  for(int i = 0; i < 8; i++) {
909  double vi1 = p->e[i]->e[0]->V();
910  double vi2 = p->e[i]->e[1]->V();
911  double vi3 = p->e[i]->e[2]->V();
912  double vi4 = p->e[i]->e[3]->V();
913  double vi5 = p->e[i]->e[4]->V();
914  double vi6 = p->e[i]->e[5]->V();
915  double vi7 = p->e[i]->e[6]->V();
916  double vi8 = p->e[i]->e[7]->V();
917  double vri =
918  (vi1 + vi2 + vi3 + vi4 / 2 + vi5 + vi6 + vi7 + vi8 / 2) / 7;
919  err |= (fabs((vi[i] - vri)) > AVG * tol);
920  }
921  err |= (fabs((v - vr)) > AVG * tol);
922  if(err) {
923  p->visible = false;
924  for(int i = 0; i < 8; i++) recurError(p->e[i], AVG, tol);
925  }
926  else
927  p->visible = true;
928  }
929  }
930 }
931 
932 void adaptivePyramid::create(int maxlevel)
933 {
934  cleanElement<adaptivePyramid>();
935  adaptiveVertex *p1 = adaptiveVertex::add(-1, -1, 0, allVertices);
936  adaptiveVertex *p2 = adaptiveVertex::add(1, -1, 0, allVertices);
937  adaptiveVertex *p3 = adaptiveVertex::add(1, 1, 0, allVertices);
938  adaptiveVertex *p4 = adaptiveVertex::add(-1, 1, 0, allVertices);
939  adaptiveVertex *p5 = adaptiveVertex::add(0, 0, 1, allVertices);
940  adaptivePyramid *p = new adaptivePyramid(p1, p2, p3, p4, p5);
941  recurCreate(p, maxlevel, 0);
942 }
943 
944 void adaptivePyramid::recurCreate(adaptivePyramid *p, int maxlevel, int level)
945 {
946  all.push_back(p);
947  if(level++ >= maxlevel) return;
948 
949  // quad points
950  adaptiveVertex *p1 = p->p[0];
951  adaptiveVertex *p2 = p->p[1];
952  adaptiveVertex *p3 = p->p[2];
953  adaptiveVertex *p4 = p->p[3];
954 
955  // apex
956  adaptiveVertex *p5 = p->p[4];
957 
958  // center of the quad
959 
960  adaptiveVertex *p1234 =
961  adaptiveVertex::add((p1->x + p2->x + p3->x + p4->x) * 0.25,
962  (p1->y + p2->y + p3->y + p4->y) * 0.25,
963  (p1->z + p2->z + p3->z + p4->z) * 0.25, allVertices);
964 
965  // quad edge points
966 
967  adaptiveVertex *p12 =
968  adaptiveVertex::add((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5,
969  (p1->z + p2->z) * 0.5, allVertices);
970 
971  adaptiveVertex *p23 =
972  adaptiveVertex::add((p2->x + p3->x) * 0.5, (p2->y + p3->y) * 0.5,
973  (p2->z + p3->z) * 0.5, allVertices);
974 
975  adaptiveVertex *p34 =
976  adaptiveVertex::add((p3->x + p4->x) * 0.5, (p3->y + p4->y) * 0.5,
977  (p3->z + p4->z) * 0.5, allVertices);
978 
979  adaptiveVertex *p41 =
980  adaptiveVertex::add((p4->x + p1->x) * 0.5, (p4->y + p1->y) * 0.5,
981  (p4->z + p1->z) * 0.5, allVertices);
982 
983  // quad vertex to apex edge points
984 
985  adaptiveVertex *p15 =
986  adaptiveVertex::add((p1->x + p5->x) * 0.5, (p1->y + p5->y) * 0.5,
987  (p1->z + p5->z) * 0.5, allVertices);
988 
989  adaptiveVertex *p25 =
990  adaptiveVertex::add((p2->x + p5->x) * 0.5, (p2->y + p5->y) * 0.5,
991  (p2->z + p5->z) * 0.5, allVertices);
992 
993  adaptiveVertex *p35 =
994  adaptiveVertex::add((p3->x + p5->x) * 0.5, (p3->y + p5->y) * 0.5,
995  (p3->z + p5->z) * 0.5, allVertices);
996 
997  adaptiveVertex *p45 =
998  adaptiveVertex::add((p4->x + p5->x) * 0.5, (p4->y + p5->y) * 0.5,
999  (p4->z + p5->z) * 0.5, allVertices);
1000 
1001  // four base pyramids on the quad base
1002 
1003  p->e[0] = new adaptivePyramid(p1, p12, p1234, p41, p15);
1004  recurCreate(p->e[0], maxlevel, level);
1005  p->e[1] = new adaptivePyramid(p2, p23, p1234, p12, p25);
1006  recurCreate(p->e[1], maxlevel, level);
1007  p->e[2] = new adaptivePyramid(p3, p34, p1234, p23, p35);
1008  recurCreate(p->e[2], maxlevel, level);
1009  p->e[3] = new adaptivePyramid(p4, p41, p1234, p34, p45);
1010  recurCreate(p->e[3], maxlevel, level);
1011 
1012  // top pyramids
1013 
1014  p->e[4] = new adaptivePyramid(p15, p25, p35, p45, p5);
1015  recurCreate(p->e[4], maxlevel, level);
1016  p->e[5] = new adaptivePyramid(p15, p45, p35, p25, p1234);
1017  recurCreate(p->e[5], maxlevel, level);
1018 
1019  // degenerated pyramids to replace the remaining tetrahedral holes
1020  // degenerated quad in the interior of the element, apices on the quad edges
1021 
1022  p->e[6] = new adaptivePyramid(p1234, p25, p15, p1234, p12);
1023  recurCreate(p->e[6], maxlevel, level);
1024  p->e[7] = new adaptivePyramid(p1234, p35, p25, p1234, p23);
1025  recurCreate(p->e[7], maxlevel, level);
1026  p->e[8] = new adaptivePyramid(p1234, p45, p35, p1234, p34);
1027  recurCreate(p->e[8], maxlevel, level);
1028  p->e[9] = new adaptivePyramid(p1234, p15, p45, p1234, p41);
1029  recurCreate(p->e[9], maxlevel, level);
1030 }
1031 
1032 void adaptivePyramid::error(double AVG, double tol)
1033 {
1034  adaptivePyramid *p = *all.begin();
1035  recurError(p, AVG, tol);
1036 }
1037 
1038 void adaptivePyramid::recurError(adaptivePyramid *p, double AVG, double tol)
1039 {
1040  if(!p->e[0])
1041  p->visible = true;
1042  else {
1043  double vi[10];
1044  for(int i = 0; i < 10; i++) vi[i] = p->e[i]->V();
1045  double vr = 0;
1046  for(int i = 0; i < 6; i++) vr += vi[i]; // pyramids have volume V/8
1047  for(int i = 6; i < 10; i++)
1048  vr += vi[i] * 0.5; // tetrahedra have volume V/16
1049  vr /= 8.;
1050  const double v = p->V();
1051  if(!p->e[0]->e[0]) {
1052  if(fabs(v - vr) > AVG * tol) {
1053  p->visible = false;
1054  for(int i = 0; i < 10; i++) recurError(p->e[i], AVG, tol);
1055  }
1056  else
1057  p->visible = true;
1058  }
1059  else {
1060  bool err = false;
1061  for(int i = 0; i < 10; i++) {
1062  double vj[10];
1063  for(int j = 0; j < 10; j++) vj[j] = p->e[i]->e[j]->V();
1064  double vri = 0;
1065  for(int j = 0; j < 6; j++) vri += vj[j];
1066  for(int j = 6; j < 10; j++) vri += vj[j] * 0.5;
1067  vri /= 8.;
1068  err |= (fabs((vi[i] - vri)) > AVG * tol);
1069  }
1070  err |= (fabs((v - vr)) > AVG * tol);
1071  if(err) {
1072  p->visible = false;
1073  for(int i = 0; i < 10; i++) recurError(p->e[i], AVG, tol);
1074  }
1075  else
1076  p->visible = true;
1077  }
1078  }
1079 }
1080 
1081 template <class T>
1082 adaptiveElements<T>::adaptiveElements(std::vector<fullMatrix<double> *> &p)
1083  : _coeffsVal(nullptr), _eexpsVal(nullptr), _interpolVal(nullptr),
1084  _coeffsGeom(nullptr), _eexpsGeom(nullptr), _interpolGeom(nullptr)
1085 {
1086  if(p.size() >= 2) {
1087  _coeffsVal = p[0];
1088  _eexpsVal = p[1];
1089  }
1090  if(p.size() == 4) {
1091  _coeffsGeom = p[2];
1092  _eexpsGeom = p[3];
1093  }
1094 }
1095 
1096 template <class T> adaptiveElements<T>::~adaptiveElements()
1097 {
1098  if(_interpolVal) delete _interpolVal;
1099  if(_interpolGeom) delete _interpolGeom;
1100  cleanElement<T>();
1101 }
1102 
1103 template <class T> void adaptiveElements<T>::init(int level)
1104 {
1105 #ifdef TIMER
1106  double t1 = TimeOfDay();
1107 #endif
1108 
1109  T::create(level);
1110  int numVals = _coeffsVal ? _coeffsVal->size1() : T::numNodes;
1111  int numNodes = _coeffsGeom ? _coeffsGeom->size1() : T::numNodes;
1112 
1113  if(_interpolVal) delete _interpolVal;
1114  _interpolVal = new fullMatrix<double>(T::allVertices.size(), numVals);
1115 
1116  if(_interpolGeom) delete _interpolGeom;
1117  _interpolGeom = new fullMatrix<double>(T::allVertices.size(), numNodes);
1118 
1119  fullVector<double> sfv(numVals), *tmpv = nullptr;
1120  fullVector<double> sfg(numNodes), *tmpg = nullptr;
1121  if(_eexpsVal) tmpv = new fullVector<double>(_eexpsVal->size1());
1122  if(_eexpsGeom) tmpg = new fullVector<double>(_eexpsGeom->size1());
1123 
1124  int i = 0;
1125  for(auto it = T::allVertices.begin(); it != T::allVertices.end(); ++it) {
1126  if(_coeffsVal && _eexpsVal)
1127  computeShapeFunctions(_coeffsVal, _eexpsVal, it->x, it->y, it->z, &sfv,
1128  tmpv);
1129  else
1130  T::GSF(it->x, it->y, it->z, sfv);
1131  for(int j = 0; j < numVals; j++) (*_interpolVal)(i, j) = sfv(j);
1132 
1133  if(_coeffsGeom && _eexpsGeom)
1134  computeShapeFunctions(_coeffsGeom, _eexpsGeom, it->x, it->y, it->z, &sfg,
1135  tmpg);
1136  else
1137  T::GSF(it->x, it->y, it->z, sfg);
1138  for(int j = 0; j < numNodes; j++) (*_interpolGeom)(i, j) = sfg(j);
1139 
1140  i++;
1141  }
1142 
1143  if(tmpv) delete tmpv;
1144  if(tmpg) delete tmpg;
1145 
1146 #ifdef TIMER
1147  adaptiveData::timerInit += TimeOfDay() - t1;
1148  return;
1149 #endif
1150 }
1151 
1152 template <> void adaptiveElements<adaptivePyramid>::init(int level)
1153 {
1154 #ifdef TIMER
1155  double t1 = TimeOfDay();
1156 #endif
1157 
1158  adaptivePyramid::create(level);
1159  int numVals = _coeffsVal ? _coeffsVal->size1() : adaptivePyramid::numNodes;
1160  int numNodes = _coeffsGeom ? _coeffsGeom->size1() : adaptivePyramid::numNodes;
1161 
1162  if(_interpolVal) delete _interpolVal;
1163  _interpolVal =
1164  new fullMatrix<double>(adaptivePyramid::allVertices.size(), numVals);
1165 
1166  if(_interpolGeom) delete _interpolGeom;
1167  _interpolGeom =
1168  new fullMatrix<double>(adaptivePyramid::allVertices.size(), numNodes);
1169 
1170  fullVector<double> sfv(numVals), *tmpv = nullptr;
1171  fullVector<double> sfg(numNodes), *tmpg = nullptr;
1172  if(_eexpsVal) tmpv = new fullVector<double>(_eexpsVal->size1());
1173  if(_eexpsGeom) tmpg = new fullVector<double>(_eexpsGeom->size1());
1174 
1175  int i = 0;
1176  for(auto it = adaptivePyramid::allVertices.begin();
1177  it != adaptivePyramid::allVertices.end(); ++it) {
1178  if(_coeffsVal && _eexpsVal)
1179  computeShapeFunctionsPyramid(_coeffsVal, _eexpsVal, it->x, it->y, it->z,
1180  &sfv, tmpv);
1181  else
1182  adaptivePyramid::GSF(it->x, it->y, it->z, sfv);
1183 
1184  for(int j = 0; j < numVals; j++) (*_interpolVal)(i, j) = sfv(j);
1185 
1186  if(_coeffsGeom && _eexpsGeom)
1187  computeShapeFunctionsPyramid(_coeffsGeom, _eexpsGeom, it->x, it->y, it->z,
1188  &sfg, tmpg);
1189  else
1190  adaptivePyramid::GSF(it->x, it->y, it->z, sfg);
1191  for(int j = 0; j < numNodes; j++) (*_interpolGeom)(i, j) = sfg(j);
1192 
1193  i++;
1194  }
1195 
1196  if(tmpv) delete tmpv;
1197  if(tmpg) delete tmpg;
1198 
1199 #ifdef TIMER
1200  adaptiveData::timerInit += TimeOfDay() - t1;
1201  return;
1202 #endif
1203 }
1204 
1205 template <class T>
1206 bool adaptiveElements<T>::adapt(double tol, int numComp,
1207  std::vector<PCoords> &coords,
1208  std::vector<PValues> &values, double &minVal,
1209  double &maxVal, GMSH_PostPlugin *plug,
1210  bool onlyComputeMinMax)
1211 {
1212  int numVertices = T::allVertices.size();
1213 
1214  if(!numVertices) {
1215  Msg::Warning("No adapted vertices to interpolate");
1216  return false;
1217  }
1218 
1219  int numVals = _coeffsVal ? _coeffsVal->size1() : T::numNodes;
1220 
1221  if(numVals != (int)values.size()) {
1222  Msg::Warning("Wrong number of values in adaptation %d != %i", numVals,
1223  values.size());
1224  return false;
1225  }
1226 
1227 #ifdef TIMER
1228  double t1 = TimeOfDay();
1229 #endif
1230 
1231  fullVector<double> val(numVals), res(numVertices);
1232  switch(numComp) {
1233  case 1: {
1234  for(int i = 0; i < numVals; i++) val(i) = values[i].v[0];
1235  break;
1236  }
1237  case 3:
1238  case 9: {
1239  for(int i = 0; i < numVals; i++) {
1240  val(i) = 0;
1241  for(int k = 0; k < numComp; k++)
1242  val(i) += values[i].v[k] * values[i].v[k];
1243  }
1244  break;
1245  }
1246  default: {
1247  Msg::Error("Can only adapt scalar, vector or tensor data");
1248  return false;
1249  }
1250  }
1251 
1252  _interpolVal->mult(val, res);
1253 
1254  // minVal = VAL_INF;
1255  // maxVal = -VAL_INF;
1256  for(int i = 0; i < numVertices; i++) {
1257  minVal = std::min(minVal, res(i));
1258  maxVal = std::max(maxVal, res(i));
1259  }
1260  if(onlyComputeMinMax) return true;
1261 
1262  fullMatrix<double> *resxyz = nullptr;
1263  if(numComp == 3 || numComp == 9) {
1264  fullMatrix<double> valxyz(numVals, numComp);
1265  resxyz = new fullMatrix<double>(numVertices, numComp);
1266  for(int i = 0; i < numVals; i++) {
1267  for(int k = 0; k < numComp; k++) { valxyz(i, k) = values[i].v[k]; }
1268  }
1269  _interpolVal->mult(valxyz, *resxyz);
1270  }
1271 
1272  int numNodes = _coeffsGeom ? _coeffsGeom->size1() : T::numNodes;
1273  if(numNodes != (int)coords.size()) {
1274  Msg::Error("Wrong number of nodes in adaptation %d != %i", numNodes,
1275  coords.size());
1276  if(resxyz) delete resxyz;
1277  return false;
1278  }
1279 
1280  fullMatrix<double> xyz(numNodes, 3), XYZ(numVertices, 3);
1281  for(int i = 0; i < numNodes; i++) {
1282  xyz(i, 0) = coords[i].c[0];
1283  xyz(i, 1) = coords[i].c[1];
1284  xyz(i, 2) = coords[i].c[2];
1285  }
1286  _interpolGeom->mult(xyz, XYZ);
1287 
1288 #ifdef TIMER
1289  adaptiveData::timerAdapt += TimeOfDay() - t1;
1290  return true;
1291 #endif
1292 
1293  int i = 0;
1294  for(auto it = T::allVertices.begin(); it != T::allVertices.end(); ++it) {
1295  // ok because we know this will not change the set ordering
1296  adaptiveVertex *p = (adaptiveVertex *)&(*it);
1297  p->val = res(i);
1298  if(resxyz) {
1299  p->val = (*resxyz)(i, 0);
1300  p->valy = (*resxyz)(i, 1);
1301  p->valz = (*resxyz)(i, 2);
1302  if(numComp == 9) {
1303  p->valyx = (*resxyz)(i, 3);
1304  p->valyy = (*resxyz)(i, 4);
1305  p->valyz = (*resxyz)(i, 5);
1306  p->valzx = (*resxyz)(i, 6);
1307  p->valzy = (*resxyz)(i, 7);
1308  p->valzz = (*resxyz)(i, 8);
1309  }
1310  }
1311  p->X = XYZ(i, 0);
1312  p->Y = XYZ(i, 1);
1313  p->Z = XYZ(i, 2);
1314  i++;
1315  }
1316 
1317  if(resxyz) delete resxyz;
1318 
1319  for(auto it = T::all.begin(); it != T::all.end(); it++)
1320  (*it)->visible = false;
1321 
1322  if(!plug || tol != 0.) {
1323  double avg = fabs(maxVal - minVal);
1324  if(tol < 0) avg = 1.; // force visibility to the smallest subdivision
1325  T::error(avg, tol);
1326  }
1327 
1328  if(plug) plug->assignSpecificVisibility();
1329 
1330  coords.clear();
1331  values.clear();
1332  for(auto it = T::all.begin(); it != T::all.end(); it++) {
1333  if((*it)->visible) {
1334  adaptiveVertex **p = (*it)->p;
1335  for(int i = 0; i < T::numNodes; i++) {
1336  coords.push_back(PCoords(p[i]->X, p[i]->Y, p[i]->Z));
1337  switch(numComp) {
1338  case 1: values.push_back(PValues(p[i]->val)); break;
1339  case 3:
1340  values.push_back(PValues(p[i]->val, p[i]->valy, p[i]->valz));
1341  break;
1342  case 9:
1343  values.push_back(PValues(p[i]->val, p[i]->valy, p[i]->valz,
1344  p[i]->valyx, p[i]->valyy, p[i]->valyz,
1345  p[i]->valzx, p[i]->valzy, p[i]->valzz));
1346  break;
1347  }
1348  }
1349  }
1350  }
1351 
1352  return true;
1353 }
1354 
1355 template <class T>
1356 void adaptiveElements<T>::addInView(double tol, int step, PViewData *in,
1357  PViewDataList *out, GMSH_PostPlugin *plug)
1358 {
1359  int numComp = in->getNumComponents(0, 0, 0);
1360  if(numComp != 1 && numComp != 3 && numComp != 9) return;
1361 
1362  int numEle = 0, *outNb = nullptr;
1363  std::vector<double> *outList = nullptr;
1364  switch(T::numEdges) {
1365  case 0:
1366  numEle = in->getNumPoints();
1367  outNb =
1368  (numComp == 1) ? &out->NbSP : ((numComp == 3) ? &out->NbVP : &out->NbTP);
1369  outList =
1370  (numComp == 1) ? &out->SP : ((numComp == 3) ? &out->VP : &out->TP);
1371  break;
1372  case 1:
1373  numEle = in->getNumLines();
1374  outNb =
1375  (numComp == 1) ? &out->NbSL : ((numComp == 3) ? &out->NbVL : &out->NbTL);
1376  outList =
1377  (numComp == 1) ? &out->SL : ((numComp == 3) ? &out->VL : &out->TL);
1378  break;
1379  case 3:
1380  numEle = in->getNumTriangles();
1381  outNb =
1382  (numComp == 1) ? &out->NbST : ((numComp == 3) ? &out->NbVT : &out->NbTT);
1383  outList =
1384  (numComp == 1) ? &out->ST : ((numComp == 3) ? &out->VT : &out->TT);
1385  break;
1386  case 4:
1387  numEle = in->getNumQuadrangles();
1388  outNb =
1389  (numComp == 1) ? &out->NbSQ : ((numComp == 3) ? &out->NbVQ : &out->NbTQ);
1390  outList =
1391  (numComp == 1) ? &out->SQ : ((numComp == 3) ? &out->VQ : &out->TQ);
1392  break;
1393  case 6:
1394  numEle = in->getNumTetrahedra();
1395  outNb =
1396  (numComp == 1) ? &out->NbSS : ((numComp == 3) ? &out->NbVS : &out->NbTS);
1397  outList =
1398  (numComp == 1) ? &out->SS : ((numComp == 3) ? &out->VS : &out->TS);
1399  break;
1400  case 9:
1401  numEle = in->getNumPrisms();
1402  outNb =
1403  (numComp == 1) ? &out->NbSI : ((numComp == 3) ? &out->NbVI : &out->NbTI);
1404  outList =
1405  (numComp == 1) ? &out->SI : ((numComp == 3) ? &out->VI : &out->TI);
1406  break;
1407  case 8:
1408  numEle = in->getNumPyramids();
1409  outNb =
1410  (numComp == 1) ? &out->NbSY : ((numComp == 3) ? &out->NbVY : &out->NbTY);
1411  outList =
1412  (numComp == 1) ? &out->SY : ((numComp == 3) ? &out->VY : &out->TY);
1413  break;
1414  case 12:
1415  numEle = in->getNumHexahedra();
1416  outNb =
1417  (numComp == 1) ? &out->NbSH : ((numComp == 3) ? &out->NbVH : &out->NbTH);
1418  outList =
1419  (numComp == 1) ? &out->SH : ((numComp == 3) ? &out->VH : &out->TH);
1420  break;
1421  }
1422  if(!numEle) return;
1423 
1424  outList->clear();
1425  *outNb = 0;
1426 
1427  for(int ent = 0; ent < in->getNumEntities(step); ent++) {
1428  for(int ele = 0; ele < in->getNumElements(step, ent); ele++) {
1429  if(in->skipElement(step, ent, ele) ||
1430  in->getNumEdges(step, ent, ele) != T::numEdges)
1431  continue;
1432  int numNodes = in->getNumNodes(step, ent, ele);
1433  std::vector<PCoords> coords;
1434  for(int i = 0; i < numNodes; i++) {
1435  double x, y, z;
1436  in->getNode(step, ent, ele, i, x, y, z);
1437  coords.push_back(PCoords(x, y, z));
1438  }
1439  int numVal = in->getNumValues(step, ent, ele);
1440  std::vector<PValues> values;
1441 
1442  switch(numComp) {
1443  case 1:
1444  for(int i = 0; i < numVal; i++) {
1445  double val;
1446  in->getValue(step, ent, ele, i, val);
1447  values.push_back(PValues(val));
1448  }
1449  break;
1450  case 3: {
1451  for(int i = 0; i < numVal / 3; i++) {
1452  double vx, vy, vz;
1453  in->getValue(step, ent, ele, 3 * i + 0, vx);
1454  in->getValue(step, ent, ele, 3 * i + 1, vy);
1455  in->getValue(step, ent, ele, 3 * i + 2, vz);
1456  values.push_back(PValues(vx, vy, vz));
1457  }
1458  break;
1459  }
1460  case 9: {
1461  for(int i = 0; i < numVal / 9; i++) {
1462  double vxx, vxy, vxz, vyx, vyy, vyz, vzx, vzy, vzz;
1463  in->getValue(step, ent, ele, 9 * i + 0, vxx);
1464  in->getValue(step, ent, ele, 9 * i + 1, vxy);
1465  in->getValue(step, ent, ele, 9 * i + 2, vxz);
1466  in->getValue(step, ent, ele, 9 * i + 3, vyx);
1467  in->getValue(step, ent, ele, 9 * i + 4, vyy);
1468  in->getValue(step, ent, ele, 9 * i + 5, vyz);
1469  in->getValue(step, ent, ele, 9 * i + 6, vzx);
1470  in->getValue(step, ent, ele, 9 * i + 7, vzy);
1471  in->getValue(step, ent, ele, 9 * i + 8, vzz);
1472  values.push_back(
1473  PValues(vxx, vxy, vxz, vyx, vyy, vyz, vzx, vzy, vzz));
1474  }
1475  break;
1476  }
1477  }
1478  if(adapt(tol, numComp, coords, values, out->Min, out->Max, plug)) {
1479  *outNb += coords.size() / T::numNodes;
1480  for(std::size_t i = 0; i < coords.size() / T::numNodes; i++) {
1481  for(int k = 0; k < T::numNodes; ++k)
1482  outList->push_back(coords[T::numNodes * i + k].c[0]);
1483  for(int k = 0; k < T::numNodes; ++k)
1484  outList->push_back(coords[T::numNodes * i + k].c[1]);
1485  for(int k = 0; k < T::numNodes; ++k)
1486  outList->push_back(coords[T::numNodes * i + k].c[2]);
1487  for(int k = 0; k < T::numNodes; ++k)
1488  for(int l = 0; l < numComp; ++l)
1489  outList->push_back(values[T::numNodes * i + k].v[l]);
1490  }
1491  }
1492  }
1493  }
1494 }
1495 
1496 adaptiveData::adaptiveData(PViewData *data, bool outDataInit)
1497  : _step(-1), _level(-1), _tol(-1.), _inData(data), _points(nullptr),
1498  _lines(nullptr), _triangles(nullptr), _quadrangles(nullptr),
1499  _tetrahedra(nullptr), _hexahedra(nullptr), _prisms(nullptr),
1500  _pyramids(nullptr)
1501 {
1502  if(outDataInit ==
1503  true) { // For visualization of the adapted view in GMSH GUI only
1504  _outData = new PViewDataList(true);
1505  _outData->setName(data->getName() + "_adapted");
1506  }
1507  else {
1508  _outData = nullptr; // For external used
1509  }
1510  std::vector<fullMatrix<double> *> p;
1511  if(_inData->getNumPoints()) {
1512  _inData->getInterpolationMatrices(TYPE_PNT, p);
1513  _points = new adaptiveElements<adaptivePoint>(p);
1514  }
1515  if(_inData->getNumLines()) {
1516  _inData->getInterpolationMatrices(TYPE_LIN, p);
1517  _lines = new adaptiveElements<adaptiveLine>(p);
1518  }
1519  if(_inData->getNumTriangles()) {
1520  _inData->getInterpolationMatrices(TYPE_TRI, p);
1521  _triangles = new adaptiveElements<adaptiveTriangle>(p);
1522  }
1523  if(_inData->getNumQuadrangles()) {
1524  _inData->getInterpolationMatrices(TYPE_QUA, p);
1525  _quadrangles = new adaptiveElements<adaptiveQuadrangle>(p);
1526  }
1527  if(_inData->getNumTetrahedra()) {
1528  _inData->getInterpolationMatrices(TYPE_TET, p);
1529  _tetrahedra = new adaptiveElements<adaptiveTetrahedron>(p);
1530  }
1531  if(_inData->getNumPrisms()) {
1532  _inData->getInterpolationMatrices(TYPE_PRI, p);
1533  _prisms = new adaptiveElements<adaptivePrism>(p);
1534  }
1535  if(_inData->getNumHexahedra()) {
1536  _inData->getInterpolationMatrices(TYPE_HEX, p);
1537  _hexahedra = new adaptiveElements<adaptiveHexahedron>(p);
1538  }
1539  if(_inData->getNumPyramids()) {
1540  _inData->getInterpolationMatrices(TYPE_PYR, p);
1541  _pyramids = new adaptiveElements<adaptivePyramid>(p);
1542  }
1543  upWriteVTK(true); // By default, write VTK data if called...
1544  upBuildStaticData(false); // ... and do not generated global static data
1545  // structure (only useful for ParaView plugin).
1546 }
1547 
1548 adaptiveData::~adaptiveData()
1549 {
1550  if(_points) delete _points;
1551  if(_lines) delete _lines;
1552  if(_triangles) delete _triangles;
1553  if(_quadrangles) delete _quadrangles;
1554  if(_tetrahedra) delete _tetrahedra;
1555  if(_prisms) delete _prisms;
1556  if(_hexahedra) delete _hexahedra;
1557  if(_pyramids) delete _pyramids;
1558  if(_outData) delete _outData;
1559 }
1560 
1561 double adaptiveData::timerInit = 0.;
1562 double adaptiveData::timerAdapt = 0.;
1563 
1564 void adaptiveData::changeResolution(int step, int level, double tol,
1565  GMSH_PostPlugin *plug)
1566 {
1567  timerInit = timerAdapt = 0.;
1568 
1569  if(_level != level) {
1570  if(_points) _points->init(level);
1571  if(_lines) _lines->init(level);
1572  if(_triangles) _triangles->init(level);
1573  if(_quadrangles) _quadrangles->init(level);
1574  if(_tetrahedra) _tetrahedra->init(level);
1575  if(_prisms) _prisms->init(level);
1576  if(_hexahedra) _hexahedra->init(level);
1577  if(_pyramids) _pyramids->init(level);
1578  }
1579  if(plug || _step != step || _level != level || _tol != tol) {
1580  _outData->setDirty(true);
1581  if(_points) _points->addInView(tol, step, _inData, _outData, plug);
1582  if(_lines) _lines->addInView(tol, step, _inData, _outData, plug);
1583  if(_triangles) _triangles->addInView(tol, step, _inData, _outData, plug);
1584  if(_quadrangles)
1585  _quadrangles->addInView(tol, step, _inData, _outData, plug);
1586  if(_tetrahedra) _tetrahedra->addInView(tol, step, _inData, _outData, plug);
1587  if(_prisms) _prisms->addInView(tol, step, _inData, _outData, plug);
1588  if(_hexahedra) _hexahedra->addInView(tol, step, _inData, _outData, plug);
1589  if(_pyramids) _pyramids->addInView(tol, step, _inData, _outData, plug);
1590  _outData->finalize();
1591  }
1592  _step = step;
1593  _level = level;
1594  _tol = tol;
1595 
1596 #ifdef TIMER
1597  printf("init time = %g\n", timerInit);
1598  printf("adapt time = %g\n", timerAdapt);
1599 #endif
1600 }
1601 
1602 bool VTKData::isLittleEndian()
1603 {
1604  int num = 1;
1605  if(*(char *)&num == 1)
1606  return true; // Little Endian
1607  else
1608  return false; // Big Endian
1609 }
1610 
1611 void VTKData::SwapArrayByteOrder(void *array, int nbytes, int nItems)
1612 {
1613  // This swaps the byte order for the array of nItems each of size nbytes
1614  int i, j;
1615  unsigned char *ucDst = (unsigned char *)array;
1616 
1617  for(i = 0; i < nItems; i++) {
1618  for(j = 0; j < (nbytes / 2); j++)
1619  std::swap(ucDst[j], ucDst[(nbytes - 1) - j]);
1620  ucDst += nbytes;
1621  }
1622 }
1623 
1624 void VTKData::writeVTKElmData()
1625 {
1626  // This routine writes vtu files (ascii or binary) from a elemental data base
1627  // of nodes coordinates, cell connectivity, type and offset, and point data
1628  // (either scalar or vector field)
1629 
1630  // Format choice
1631  if(vtkFormat == "vtu") {
1632  if(vtkCountTotElmLev0 <= numPartMinElm * minElmPerPart) {
1633  if((vtkCountTotElmLev0 - 1) % minElmPerPart == 0) { // new filename
1634  vtkCountFile = (vtkCountTotElmLev0 - 1) / minElmPerPart;
1635  initVTKFile();
1636  }
1637  }
1638  else {
1639  if((vtkCountTotElmLev0 - 1 - numPartMinElm * minElmPerPart) %
1640  maxElmPerPart ==
1641  0) {
1642  // new filename
1643  vtkCountFile = numPartMinElm + (vtkCountTotElmLev0 - 1 -
1644  numPartMinElm * minElmPerPart) /
1645  maxElmPerPart;
1646  initVTKFile();
1647  }
1648  }
1649 
1650  if(vtkIsBinary == true) { // Use appended format for raw binary
1651 
1652  // Write raw binary data to separate files first. Text headers will be
1653  // added later, as wall as raw data size (needs to know the size before)
1654 
1655  int counter;
1656  uint64_t *i64array;
1657  uint8_t *i8array;
1658  double *darray;
1659 
1660  // Node value
1661  counter = 0;
1662  darray = new double[vtkNumComp * vtkLocalValues.size()];
1663  for(auto it = vtkLocalValues.begin(); it != vtkLocalValues.end(); ++it) {
1664  for(int i = 0; i < vtkNumComp; i++) { darray[counter + i] = it->v[i]; }
1665  counter += vtkNumComp;
1666  vtkCountTotVal += vtkNumComp;
1667  }
1668  assert(counter == vtkNumComp * (int)vtkLocalValues.size());
1669  fwrite(darray, sizeof(double), vtkNumComp * vtkLocalValues.size(),
1670  vtkFileNodVal);
1671  delete[] darray;
1672 
1673  // Points
1674  int sizeArray = (int)vtkLocalCoords.size();
1675  darray = new double[3 * sizeArray];
1676  counter = 0;
1677  for(auto it = vtkLocalCoords.begin(); it != vtkLocalCoords.end(); ++it) {
1678  for(int i = 0; i < 3; i++) { darray[counter + i] = (*it).c[i]; }
1679  counter += 3;
1680  vtkCountCoord += 3;
1681  }
1682  fwrite(darray, sizeof(double), 3 * sizeArray, vtkFileCoord);
1683  delete[] darray;
1684 
1685  // Cells
1686 
1687  // First count the number of integers that described the cell data in
1688  // vtkConnectivity See
1689  // http://www.vtk.org/wp-content/uploads/2015/04/file-formats.pdf (page 4)
1690  int cellSizeData = 0;
1691  for(auto it = vtkLocalConnectivity.begin();
1692  it != vtkLocalConnectivity.end(); ++it) {
1693  // Contrary to vtk format, no +1 required for the number of nodes in the
1694  // element
1695  cellSizeData += (int)it->size();
1696  }
1697 
1698  // Connectivity (and build offset at the same time)
1699  counter = 0;
1700  int cellcounter = 0;
1701  i64array = new uint64_t[cellSizeData];
1702  uint64_t *cellOffset = new uint64_t[vtkLocalConnectivity.size()];
1703  for(auto it = vtkLocalConnectivity.begin();
1704  it != vtkLocalConnectivity.end(); ++it) {
1705  for(auto jt = it->begin(); jt != it->end(); ++jt) {
1706  i64array[counter] = *jt;
1707  counter++;
1708  vtkCountTotNodConnect++;
1709  }
1710  cellOffset[cellcounter] = vtkCountTotNodConnect; // build the offset
1711  cellcounter++;
1712  }
1713  fwrite(i64array, sizeof(uint64_t), cellSizeData, vtkFileConnect);
1714  delete[] i64array;
1715 
1716  // Cell offset
1717  fwrite(cellOffset, sizeof(uint64_t), vtkLocalConnectivity.size(),
1718  vtkFileCellOffset);
1719  delete[] cellOffset;
1720 
1721  // Cell type
1722  counter = 0;
1723  i8array = new uint8_t[vtkLocalConnectivity.size()];
1724  for(auto it = vtkLocalCellType.begin(); it != vtkLocalCellType.end();
1725  it++) {
1726  i8array[counter] = *it;
1727  counter++;
1728  }
1729  fwrite(i8array, sizeof(uint8_t), vtkLocalConnectivity.size(),
1730  vtkFileCellType);
1731  delete[] i8array;
1732  }
1733  else { // ascii
1734 
1735  // Node values
1736  for(auto it = vtkLocalValues.begin(); it != vtkLocalValues.end(); ++it) {
1737  for(int i = 0; i < vtkNumComp; i++) {
1738  fprintf(vtkFileNodVal, "%23.16e ", (*it).v[i]);
1739  vtkCountTotVal++;
1740  if(vtkCountTotVal % 6 == 0) fprintf(vtkFileNodVal, "\n");
1741  }
1742  }
1743 
1744  for(auto it = vtkLocalCoords.begin(); it != vtkLocalCoords.end(); it++) {
1745  fprintf(vtkFileCoord, "%23.16e %23.16e %23.16e ", (*it).c[0],
1746  (*it).c[1], (*it).c[2]);
1747  vtkCountCoord += 3;
1748  if(vtkCountCoord % 6 == 0) fprintf(vtkFileCoord, "\n");
1749  }
1750 
1751  // Cells
1752  // Connectivity
1753  int *cellOffset = new int[vtkLocalConnectivity.size()];
1754  int cellcounter = 0;
1755  for(auto it = vtkLocalConnectivity.begin();
1756  it != vtkLocalConnectivity.end(); ++it) {
1757  for(auto jt = it->begin(); jt != it->end(); ++jt) {
1758  fprintf(vtkFileConnect, "%d ", *jt);
1759  vtkCountTotNodConnect++;
1760  if(vtkCountTotNodConnect % 6 == 0) fprintf(vtkFileConnect, "\n");
1761  }
1762  cellOffset[cellcounter] = vtkCountTotNodConnect; // build the offset
1763  cellcounter++;
1764  }
1765 
1766  // Cell offset
1767  for(uint64_t i = 0; i < vtkLocalConnectivity.size(); i++) {
1768  fprintf(vtkFileCellOffset, "%d ", cellOffset[i]);
1769  vtkCountCellOffset++;
1770  if(vtkCountCellOffset % 6 == 0) fprintf(vtkFileCellOffset, "\n");
1771  }
1772  delete[] cellOffset;
1773 
1774  // Cell type
1775  for(auto it = vtkLocalCellType.begin(); it != vtkLocalCellType.end();
1776  it++) {
1777  fprintf(vtkFileCellType, "%d ", *it);
1778  vtkCountCellType++;
1779  if(vtkCountCellType % 6 == 0) fprintf(vtkFileCellType, "\n");
1780  }
1781 
1782  } // if ascii
1783 
1784  // finalize and close current vtu file
1785  if(vtkCountTotElmLev0 <= numPartMinElm * minElmPerPart) {
1786  if(vtkCountTotElmLev0 % minElmPerPart == 0) { finalizeVTKFile(); }
1787  }
1788  else {
1789  if((vtkCountTotElmLev0 - numPartMinElm * minElmPerPart) % maxElmPerPart ==
1790  0) {
1791  finalizeVTKFile();
1792  }
1793  }
1794 
1795  } // vtu format
1796  else
1797  Msg::Error("Unknown format");
1798 }
1799 
1800 void VTKData::initVTKFile()
1801 {
1802  // Temporary files
1803  vtkFileCoord = fopen("vtkCoords.vtu", "wb");
1804  vtkFileConnect = fopen("vtkConnectivity.vtu", "wb");
1805  vtkFileCellOffset = fopen("vtkCellOffset.vtu", "wb");
1806  vtkFileCellType = fopen("vtkCellType.vtu", "wb");
1807  vtkFileNodVal = fopen("vtkNodeValue.vtu", "wb");
1808 
1809  if(vtkCountFile == 0) {
1810  // write the pvtu file and create the corresponding directory for vtu files
1811 
1812  if(vtkUseDefaultName == 1) {
1813  vtkDirName = vtkFieldName + "_step" + ToString<int>(vtkStep) + "_level" +
1814  ToString<int>(vtkLevel) + "_tol" + ToString<double>(vtkTol) +
1815  "_npart" + ToString<int>(vtkNpart);
1816  }
1817  else {
1818  // Remove existing extension here to avoid duplicate
1819  std::size_t found = vtkFileName.find_last_of('.');
1820  // remove extension
1821  if(found != std::string::npos) vtkFileName = vtkFileName.substr(0, found);
1822  vtkDirName = vtkFileName;
1823  }
1824 
1825  CreateSingleDir(vtkDirName);
1826 
1827  vtkFileName =
1828  vtkDirName + ".p" + vtkFormat; // add pvtu extension to file name
1829  vtkFile = fopen(vtkFileName.c_str(), "w");
1830 
1831  bool littleEndian = isLittleEndian(); // Determine endianess
1832  if(littleEndian == true)
1833  fprintf(vtkFile, "<VTKFile type=\"PUnstructuredGrid\" version=\"1.0\" "
1834  "byte_order=\"LittleEndian\">\n");
1835  else
1836  fprintf(vtkFile, "<VTKFile type=\"PUnstructuredGrid\" version=\"1.0\" "
1837  "byte_order=\"BigEndian\">\n");
1838 
1839  fprintf(vtkFile, "<PUnstructuredGrid GhostLevel=\"0\">\n");
1840  fprintf(vtkFile, "<PPoints>\n");
1841  fprintf(vtkFile, "<DataArray type=\"Float64\" Name=\"Points\" "
1842  "NumberOfComponents=\"3\"/>\n");
1843  fprintf(vtkFile, "</PPoints>\n");
1844 
1845  fprintf(vtkFile, "<PCells>\n");
1846  fprintf(vtkFile, "<PDataArray type=\"Int64\" Name=\"connectivity\" "
1847  "NumberOfComponents=\"1\"/>\n");
1848  fprintf(vtkFile, "<PDataArray type=\"Int64\" Name=\"offsets\" "
1849  "NumberOfComponents=\"1\"/>\n");
1850  fprintf(
1851  vtkFile,
1852  "<PDataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\"/>\n");
1853  fprintf(vtkFile, "</PCells>\n");
1854 
1855  fprintf(vtkFile, "<PPointData>\n");
1856  fprintf(
1857  vtkFile,
1858  "<PDataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"%d\"/>\n",
1859  vtkFieldName.c_str(), vtkNumComp);
1860  fprintf(vtkFile, "</PPointData>\n");
1861 
1862  fprintf(vtkFile, "<PCellData>\n");
1863  fprintf(vtkFile, "</PCellData>\n");
1864 
1865  for(int i = 0; i < vtkNpart; i++)
1866  fprintf(vtkFile, "<Piece Source=\"%s/data%d.vtu\"/>\n",
1867  vtkDirName.c_str(), i);
1868  fprintf(vtkFile, "</PUnstructuredGrid>\n");
1869  fprintf(vtkFile, "</VTKFile>\n");
1870  fclose(vtkFile);
1871  }
1872 }
1873 
1874 void VTKData::finalizeVTKFile()
1875 {
1876  // This routine writes vtu files (ascii or binary) from a complete data base
1877  // of nodes coordinates, cell connectivity, type and offset, and point data
1878  // (either scalar or vector field)
1879 
1880  // Close first temporary files. Todo: Avoid multiple open/close actions and
1881  // keep the file open to write all information
1882  fclose(vtkFileCoord);
1883  fclose(vtkFileConnect);
1884  fclose(vtkFileCellOffset);
1885  fclose(vtkFileCellType);
1886  fclose(vtkFileNodVal);
1887 
1888  bool littleEndian = isLittleEndian(); // Determine endianess
1889 
1890  // Open final file
1891  std::string filename;
1892  filename = vtkDirName + "/data" + ToString(vtkCountFile) + "." + vtkFormat;
1893 
1894  Msg::StatusBar(true,
1895  "Writing VTK data in %s: fieldname = %s - numElm = %d - "
1896  "numNod = %d nodes\n",
1897  filename.c_str(), vtkFieldName.c_str(), vtkCountTotElm,
1898  vtkCountTotNod);
1899 
1900  assert(vtkCountTotNod == vtkCountCoord / 3);
1901 
1902  // Now concatenate headers with data files
1903  if(vtkFormat == "vtu") { // Format choice
1904 
1905  if(vtkIsBinary == true) { // Binary or ascii
1906 
1907  vtkFile = fopen(filename.c_str(), "wb");
1908  if(vtkFile == nullptr) {
1909  printf("Could not open file %s\n", filename.c_str());
1910  return;
1911  }
1912 
1913  uint64_t byteoffset = 0;
1914 
1915  // Headers first
1916 
1917  if(littleEndian == true)
1918  fprintf(vtkFile, "<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" "
1919  "byte_order=\"LittleEndian\" "
1920  "header_type=\"UInt64\">\n");
1921  else
1922  fprintf(vtkFile, "<VTKFile type=\"PUnstructuredGrid\" version=\"1.0\" "
1923  "byte_order=\"BigEndian\" header_type=\"UInt64\">\n");
1924  fprintf(vtkFile, "<UnstructuredGrid>\n");
1925  fprintf(vtkFile, "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",
1926  vtkCountTotNod, vtkCountTotElm);
1927 
1928  // Node value
1929  fprintf(vtkFile, "<PointData>\n");
1930  fprintf(vtkFile,
1931  "<DataArray type=\"Float64\" Name=\"%s\" "
1932  "NumberOfComponents=\"%d\" format=\"appended\" offset=\"%" PRIu64
1933  "\"/>\n",
1934  vtkFieldName.c_str(), vtkNumComp, byteoffset);
1935  fprintf(vtkFile, "</PointData>\n");
1936  byteoffset = byteoffset + (vtkCountTotNod * vtkNumComp + 1) *
1937  sizeof(double); // +1 for datasize in bytes
1938 
1939  // Cell values (none here but may change)
1940  fprintf(vtkFile, "<CellData>\n");
1941  fprintf(vtkFile,
1942  "</CellData>\n"); // no offset here because empty cell data
1943 
1944  // Nodes
1945  fprintf(vtkFile, "<Points>\n");
1946  fprintf(vtkFile,
1947  "<DataArray type=\"Float64\" Name=\"Points\" "
1948  "NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PRIu64
1949  "\"/>\n",
1950  byteoffset);
1951  fprintf(vtkFile, "</Points>\n");
1952  byteoffset = byteoffset + (vtkCountCoord + 1) *
1953  sizeof(double); // +1 for datasize in bytes
1954 
1955  // Cells
1956  fprintf(vtkFile, "<Cells>\n");
1957  fprintf(vtkFile,
1958  "<DataArray type=\"Int64\" Name=\"connectivity\" "
1959  "format=\"appended\" offset=\"%" PRIu64 "\"/>\n",
1960  byteoffset);
1961  byteoffset = byteoffset + (vtkCountTotNodConnect + 1) * sizeof(uint64_t);
1962  fprintf(vtkFile,
1963  "<DataArray type=\"Int64\" Name=\"offsets\" format=\"appended\" "
1964  "offset=\"%" PRIu64 "\"/>\n",
1965  byteoffset);
1966  byteoffset = byteoffset + (vtkCountTotElm + 1) * sizeof(uint64_t);
1967  fprintf(vtkFile,
1968  "<DataArray type=\"UInt8\" Name=\"types\" format=\"appended\" "
1969  "offset=\"%" PRIu64 "\"/>\n",
1970  byteoffset);
1971  byteoffset = byteoffset + (vtkCountTotElm + 1) * sizeof(uint8_t);
1972  fprintf(vtkFile, "</Cells>\n");
1973 
1974  fprintf(vtkFile, "</Piece>\n");
1975  fprintf(vtkFile, "</UnstructuredGrid>\n");
1976 
1977  fprintf(vtkFile, "<AppendedData encoding=\"raw\">\n");
1978  fprintf(vtkFile, "_");
1979 
1980  uint64_t datasize;
1981 
1982  // Node values
1983  datasize = vtkNumComp * vtkCountTotNod * sizeof(double);
1984  fwrite(&datasize, sizeof(uint64_t), 1, vtkFile);
1985  fclose(vtkFile);
1986 
1987  std::ifstream if_vtkNodeValue("vtkNodeValue.vtu", std::ios_base::binary);
1988  std::ofstream of_vtkfile(filename.c_str(),
1989  std::ios_base::binary | std::ios_base::app);
1990  of_vtkfile << if_vtkNodeValue.rdbuf();
1991  if_vtkNodeValue.close();
1992  of_vtkfile.close();
1993 
1994  // Points
1995  vtkFile = fopen(filename.c_str(), "ab");
1996  datasize = vtkCountTotNod * 3 * sizeof(double);
1997  fwrite(&datasize, sizeof(uint64_t), 1, vtkFile);
1998  fclose(vtkFile);
1999 
2000  std::ifstream if_vtkCoords("vtkCoords.vtu", std::ios_base::binary);
2001  of_vtkfile.open(filename.c_str(),
2002  std::ios_base::binary | std::ios_base::app);
2003  of_vtkfile << if_vtkCoords.rdbuf();
2004  if_vtkCoords.close();
2005  of_vtkfile.close();
2006 
2007  // Cells
2008  // Connectivity
2009  vtkFile = fopen(filename.c_str(), "ab");
2010  datasize = vtkCountTotNodConnect * sizeof(uint64_t);
2011  fwrite(&datasize, sizeof(uint64_t), 1, vtkFile);
2012  fclose(vtkFile);
2013 
2014  std::ifstream if_vtkConnectivity("vtkConnectivity.vtu",
2015  std::ios_base::binary);
2016  of_vtkfile.open(filename.c_str(),
2017  std::ios_base::binary | std::ios_base::app);
2018  of_vtkfile << if_vtkConnectivity.rdbuf();
2019  if_vtkConnectivity.close();
2020  of_vtkfile.close();
2021 
2022  // Cell offset
2023  vtkFile = fopen(filename.c_str(), "ab");
2024  datasize = vtkCountTotElm * sizeof(uint64_t);
2025  fwrite(&datasize, sizeof(uint64_t), 1, vtkFile);
2026  fclose(vtkFile);
2027 
2028  std::ifstream if_vtkCellOffset("vtkCellOffset.vtu",
2029  std::ios_base::binary);
2030  of_vtkfile.open(filename.c_str(),
2031  std::ios_base::binary | std::ios_base::app);
2032  of_vtkfile << if_vtkCellOffset.rdbuf();
2033  if_vtkCellOffset.close();
2034  of_vtkfile.close();
2035 
2036  // Cell type
2037  vtkFile = fopen(filename.c_str(), "ab");
2038  datasize = vtkCountTotElm * sizeof(uint8_t);
2039  fwrite(&datasize, sizeof(uint64_t), 1, vtkFile);
2040  fclose(vtkFile);
2041 
2042  std::ifstream if_vtkCellType("vtkCellType.vtu", std::ios_base::binary);
2043  of_vtkfile.open(filename.c_str(),
2044  std::ios_base::binary | std::ios_base::app);
2045  of_vtkfile << if_vtkCellType.rdbuf();
2046  if_vtkCellType.close();
2047  of_vtkfile.close();
2048 
2049  vtkFile = fopen(filename.c_str(), "ab");
2050  fprintf(vtkFile, "\n");
2051  fprintf(vtkFile, "</AppendedData>\n");
2052  fprintf(vtkFile, "</VTKFile>\n"); // for both binary and ascii
2053  fclose(vtkFile);
2054  }
2055  else { // ascii
2056 
2057  vtkFile = fopen(filename.c_str(), "w");
2058  if(vtkFile == nullptr) {
2059  printf("Could not open file %s\n", filename.c_str());
2060  return;
2061  }
2062 
2063  if(littleEndian == true)
2064  fprintf(vtkFile, "<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" "
2065  "byte_order=\"LittleEndian\" "
2066  "header_type=\"UInt64\">\n");
2067  else
2068  fprintf(vtkFile, "<VTKFile type=\"PUnstructuredGrid\" version=\"1.0\" "
2069  "byte_order=\"BigEndian\" header_type=\"UInt64\">\n");
2070  fprintf(vtkFile, "<UnstructuredGrid>\n");
2071  fprintf(vtkFile, "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",
2072  vtkCountTotNod, vtkCountTotElm);
2073 
2074  // Node values
2075  fprintf(vtkFile, "<PointData>\n");
2076  fprintf(vtkFile,
2077  "<DataArray type=\"Float64\" Name=\"%s\" "
2078  "NumberOfComponents=\"%d\" format=\"ascii\">\n",
2079  vtkFieldName.c_str(), vtkNumComp);
2080  fclose(vtkFile); // close file for binary concatenation
2081 
2082  std::ifstream if_vtkNodeValue("vtkNodeValue.vtu", std::ios_base::binary);
2083  std::ofstream of_vtkfile(filename.c_str(),
2084  std::ios_base::binary | std::ios_base::app);
2085  of_vtkfile << if_vtkNodeValue.rdbuf();
2086  if_vtkNodeValue.close();
2087  of_vtkfile.close();
2088 
2089  vtkFile = fopen(filename.c_str(), "a");
2090  fprintf(vtkFile, "</DataArray>\n");
2091  fprintf(vtkFile, "</PointData>\n");
2092 
2093  // Cell values
2094  fprintf(vtkFile, "<CellData>\n");
2095  fprintf(vtkFile, "</CellData>\n");
2096 
2097  // Nodes
2098  fprintf(vtkFile, "<Points>\n");
2099  fprintf(vtkFile, "<DataArray type=\"Float64\" Name=\"Points\" "
2100  "NumberOfComponents=\"3\" format=\"ascii\">\n");
2101  fclose(vtkFile); // close file for binary concatenation
2102 
2103  of_vtkfile.open(filename.c_str(),
2104  std::ios_base::binary | std::ios_base::app);
2105  std::ifstream if_vtkCoords("vtkCoords.vtu", std::ios_base::binary);
2106  of_vtkfile << if_vtkCoords.rdbuf();
2107  if_vtkCoords.close();
2108  of_vtkfile.close();
2109 
2110  vtkFile = fopen(filename.c_str(), "a");
2111  fprintf(vtkFile, "</DataArray>\n");
2112  fprintf(vtkFile, "</Points>\n");
2113 
2114  // Cells
2115  fprintf(vtkFile, "<Cells>\n");
2116  fprintf(
2117  vtkFile,
2118  "<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">\n");
2119  fclose(vtkFile); // close file for binary concatenation
2120 
2121  // Connectivity
2122  of_vtkfile.open(filename.c_str(),
2123  std::ios_base::binary | std::ios_base::app);
2124  std::ifstream if_vtkConnectivity("vtkConnectivity.vtu",
2125  std::ios_base::binary);
2126  of_vtkfile << if_vtkConnectivity.rdbuf();
2127  if_vtkConnectivity.close();
2128  of_vtkfile.close();
2129 
2130  vtkFile = fopen(filename.c_str(), "a");
2131  fprintf(vtkFile, "</DataArray>\n");
2132 
2133  // Cell offset
2134  fprintf(vtkFile,
2135  "<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">\n");
2136  fclose(vtkFile); // close file for binary concatenation
2137 
2138  of_vtkfile.open(filename.c_str(),
2139  std::ios_base::binary | std::ios_base::app);
2140  std::ifstream if_vtkCellOffset("vtkCellOffset.vtu",
2141  std::ios_base::binary);
2142  of_vtkfile << if_vtkCellOffset.rdbuf();
2143  if_vtkCellOffset.close();
2144  of_vtkfile.close();
2145 
2146  vtkFile = fopen(filename.c_str(), "a");
2147  fprintf(vtkFile, "</DataArray>\n");
2148 
2149  // Cell type
2150  fprintf(vtkFile,
2151  "<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n");
2152  fclose(vtkFile); // close file for binary concatenation
2153 
2154  of_vtkfile.open(filename.c_str(),
2155  std::ios_base::binary | std::ios_base::app);
2156  std::ifstream if_vtkCellType("vtkCellType.vtu", std::ios_base::binary);
2157  of_vtkfile << if_vtkCellType.rdbuf();
2158  if_vtkCellType.close();
2159  of_vtkfile.close();
2160 
2161  vtkFile = fopen(filename.c_str(), "a");
2162  fprintf(vtkFile, "</DataArray>\n");
2163  fprintf(vtkFile, "</Cells>\n");
2164 
2165  fprintf(vtkFile, "</Piece>\n");
2166  fprintf(vtkFile, "</UnstructuredGrid>\n");
2167 
2168  fprintf(vtkFile, "</VTKFile>\n"); // for both binary and ascii
2169  fclose(vtkFile);
2170  } // if binary/ascii
2171 
2172  // Remove temporary files now
2173  if(remove("vtkCoords.vtu") != 0)
2174  printf("ERROR: Could not remove vtkCoords.vtu\n");
2175  if(remove("vtkConnectivity.vtu") != 0)
2176  printf("ERROR: Could not remove vtkConnectivity.vtu\n");
2177  if(remove("vtkCellOffset.vtu") != 0)
2178  printf("ERROR: Could not remove vtkCellOffset.vtu\n");
2179  if(remove("vtkCellType.vtu") != 0)
2180  printf("ERROR: Could not remove vtkCellType.vtu\n");
2181  if(remove("vtkNodeValue.vtu") != 0)
2182  printf("ERROR: Could not remove vtkNodeValue.vtu\n");
2183 
2184  // Reset counters for next file
2185  vtkCountTotNod = 0;
2186  vtkCountTotElm = 0;
2187  vtkCountCoord = 0;
2188  vtkCountTotNodConnect = 0;
2189  vtkCountTotVal = 0;
2190  vtkCountCellOffset = 0;
2191  vtkCountCellType = 0;
2192  }
2193 
2194  else
2195  Msg::Error("File format unknown: %s", vtkFormat.c_str());
2196 }
2197 
2198 int VTKData::getPVCellType(int numEdges)
2199 {
2200  int cellType; // Convention for cell types in ParaView
2201  switch(numEdges) {
2202  case 0:
2203  printf(
2204  "WARNING: Trying to write a node to the ParaView data base and file\n");
2205  cellType = -1;
2206  break;
2207  case 1:
2208  printf(
2209  "WARNING: Trying to write a node to the ParaView data base and file\n");
2210  cellType = -2;
2211  break;
2212  case 3:
2213  cellType = 5; // 2D VTK triangle
2214  break;
2215  case 4:
2216  cellType = 9; // 2D VTK quadrangle
2217  break;
2218  case 6:
2219  cellType = 10; // 3D VTK tetrahedron
2220  break;
2221  case 9:
2222  cellType = 13; // 3D VTK prism/wedge
2223  break;
2224  case 8:
2225  cellType = 14; // 3D VTK pyramid
2226  break;
2227  case 12:
2228  cellType = 12; // 3D VTK hexahedron
2229  break;
2230  default:
2231  printf("ERROR: No cell type was detected\n");
2232  cellType = -1;
2233  break;
2234  }
2235 
2236  return cellType;
2237 }
2238 
2239 template <class T>
2240 void adaptiveElements<T>::adaptForVTK(double tol, int numComp,
2241  std::vector<PCoords> &coords,
2242  std::vector<PValues> &values,
2243  double &minVal, double &maxVal)
2244 {
2245  int numVertices = T::allVertices.size();
2246 
2247  if(!numVertices) {
2248  Msg::Error("No adapted vertices to interpolate");
2249  return;
2250  }
2251 
2252  int numVals = _coeffsVal ? _coeffsVal->size1() : T::numNodes;
2253  if(numVals != (int)values.size()) {
2254  Msg::Error("Wrong number of values in adaptation %d != %i", numVals,
2255  values.size());
2256  return;
2257  }
2258 
2259 #ifdef TIMER
2260  double t1 = TimeOfDay();
2261 #endif
2262 
2263  fullVector<double> val(numVals), res(numVertices);
2264  switch(numComp) {
2265  case 1: {
2266  for(int i = 0; i < numVals; i++) val(i) = values[i].v[0];
2267  break;
2268  }
2269  case 3:
2270  case 9: {
2271  for(int i = 0; i < numVals; i++) {
2272  val(i) = 0;
2273  for(int k = 0; k < numComp; k++)
2274  val(i) += values[i].v[k] * values[i].v[k];
2275  }
2276  break;
2277  }
2278  default: {
2279  Msg::Error("Can only adapt scalar, vector or tensor data");
2280  return;
2281  }
2282  }
2283 
2284  _interpolVal->mult(val, res);
2285 
2286  for(int i = 0; i < numVertices; i++) {
2287  minVal = std::min(minVal, res(i));
2288  maxVal = std::max(maxVal, res(i));
2289  }
2290 
2291  fullMatrix<double> *resxyz = nullptr;
2292  if(numComp == 3 || numComp == 9) {
2293  fullMatrix<double> valxyz(numVals, numComp);
2294  resxyz = new fullMatrix<double>(numVertices, numComp);
2295  for(int i = 0; i < numVals; i++) {
2296  for(int k = 0; k < numComp; k++) { valxyz(i, k) = values[i].v[k]; }
2297  }
2298  _interpolVal->mult(valxyz, *resxyz);
2299  }
2300 
2301  int numNodes = _coeffsGeom ? _coeffsGeom->size1() : T::numNodes;
2302  if(numNodes != (int)coords.size()) {
2303  Msg::Error("Wrong number of nodes in adaptation %d != %i", numNodes,
2304  coords.size());
2305  if(resxyz) delete resxyz;
2306  return;
2307  }
2308 
2309  fullMatrix<double> xyz(numNodes, 3), XYZ(numVertices, 3);
2310  for(int i = 0; i < numNodes; i++) {
2311  xyz(i, 0) = coords[i].c[0];
2312  xyz(i, 1) = coords[i].c[1];
2313  xyz(i, 2) = coords[i].c[2];
2314  }
2315  _interpolGeom->mult(xyz, XYZ);
2316 
2317 #ifdef TIMER
2318  adaptiveData::timerAdapt += TimeOfDay() - t1;
2319  return;
2320 #endif
2321 
2322  int i = 0;
2323  for(auto it = T::allVertices.begin(); it != T::allVertices.end(); ++it) {
2324  // ok because we know this will not change the set ordering
2325  adaptiveVertex *p = (adaptiveVertex *)&(*it);
2326  p->val = res(i);
2327  if(resxyz) {
2328  p->val = (*resxyz)(i, 0);
2329  p->valy = (*resxyz)(i, 1);
2330  p->valz = (*resxyz)(i, 2);
2331  if(numComp == 9) {
2332  p->valyx = (*resxyz)(i, 3);
2333  p->valyy = (*resxyz)(i, 4);
2334  p->valyz = (*resxyz)(i, 5);
2335  p->valzx = (*resxyz)(i, 6);
2336  p->valzy = (*resxyz)(i, 7);
2337  p->valzz = (*resxyz)(i, 8);
2338  }
2339  }
2340  p->X = XYZ(i, 0);
2341  p->Y = XYZ(i, 1);
2342  p->Z = XYZ(i, 2);
2343  i++;
2344  }
2345 
2346  if(resxyz) delete resxyz;
2347 
2348  for(auto it = T::all.begin(); it != T::all.end(); it++)
2349  (*it)->visible = false;
2350 
2351  if(tol != 0.) {
2352  double avg = fabs(maxVal - minVal);
2353  if(tol < 0) avg = 1.; // force visibility to the smallest subdivision
2354  T::error(avg, tol);
2355  }
2356 
2357  coords.clear();
2358  values.clear();
2359  for(auto it = T::all.begin(); it != T::all.end(); it++) {
2360  if((*it)->visible) {
2361  adaptiveVertex **p = (*it)->p;
2362  for(int i = 0; i < T::numNodes; i++) {
2363  coords.push_back(PCoords(p[i]->X, p[i]->Y, p[i]->Z));
2364  switch(numComp) {
2365  case 1: values.push_back(PValues(p[i]->val)); break;
2366  case 3:
2367  values.push_back(PValues(p[i]->val, p[i]->valy, p[i]->valz));
2368  break;
2369  case 9:
2370  values.push_back(PValues(p[i]->val, p[i]->valy, p[i]->valz,
2371  p[i]->valyx, p[i]->valyy, p[i]->valyz,
2372  p[i]->valzx, p[i]->valzy, p[i]->valzz));
2373  break;
2374  }
2375  }
2376  }
2377  }
2378 }
2379 
2380 template <class T>
2381 void adaptiveElements<T>::buildMapping(nodMap<T> &myNodMap, double tol,
2382  int &numNodInsert)
2383 {
2384  if(tol > 0.0 || myNodMap.getSize() == 0) {
2385  // Either this is not a uniform refinement and we need to rebuild the whole
2386  // mapping for each canonical element, or this is the first time we try to
2387  // build the mapping
2388 
2389  myNodMap
2390  .cleanMapping(); // Required if tol > 0 (local error based adaptation)
2391 
2392  for(auto itleaf = T::all.begin(); itleaf != T::all.end(); itleaf++) {
2393  // Visit all the leaves of the refined canonical element
2394 
2395  if((*itleaf)->visible == true) {
2396  // Find the leaves that are flagged for visibility
2397 
2398  for(int i = 0; i < T::numNodes; i++) {
2399  // Visit each nodes of the leaf (3 for triangles, 4 for quadrangle,
2400  // etc)
2401  adaptiveVertex pquery;
2402  pquery.x = (*itleaf)->p[i]->x;
2403  pquery.y = (*itleaf)->p[i]->y;
2404  pquery.z = (*itleaf)->p[i]->z;
2405  auto it = T::allVertices.find(pquery);
2406  if(it == T::allVertices.end()) {
2407  Msg::Error("Could not find adaptive Vertex in "
2408  "adaptiveElements<T>::buildMapping %f %f %f",
2409  pquery.x, pquery.y, pquery.z);
2410  }
2411  else {
2412  // Compute the distance in the list to get the mapping for
2413  // the canonical element (note std:distance returns long int
2414  int dist = (int)std::distance(T::allVertices.begin(), it);
2415  myNodMap.mapping.push_back(dist);
2416  }
2417  // quit properly if vertex not found - Should not happen though
2418  assert(it != T::allVertices.end());
2419  } // for
2420  } // if
2421  } // for
2422 
2423  if(myNodMap.mapping.size() == 0) {
2424  Msg::Error("Node mapping in buildMapping has zero size");
2425  }
2426 
2427  // Count number of unique nodes from the mapping
2428  // Use an ordered set for efficiency
2429  // This set is also used in case of partiel refinement
2430  std::set<int> uniqueNod;
2431  for(auto it = myNodMap.mapping.begin(); it != myNodMap.mapping.end();
2432  it++) {
2433  uniqueNod.insert(*it);
2434  }
2435  numNodInsert = (int)uniqueNod.size();
2436 
2437  // Renumber the elm in the mapping in case of partial refinement (when vis
2438  // tolerance > 0) so that we have a continuous numbering starting from 0
2439  // with no missing node id in the connectivity This require a new local and
2440  // temporary mapping, based on uniqueNod already generated above
2441  if(tol > 0.0) {
2442  for(auto it = myNodMap.mapping.begin(); it != myNodMap.mapping.end();
2443  ++it) {
2444  auto jt = uniqueNod.find(*it);
2445  *it = std::distance(uniqueNod.begin(), jt);
2446  }
2447  }
2448  }
2449 }
2450 
2451 template <class T>
2452 void adaptiveElements<T>::addInViewForVTK(int step, PViewData *in,
2453  VTKData &myVTKData, bool writeVTK,
2454  bool buildStaticData)
2455 {
2456  int numComp = in->getNumComponents(0, 0, 0);
2457  if(numComp != 1 && numComp != 3 && numComp != 9) return;
2458 
2459  int numEle = 0;
2460  switch(T::numEdges) {
2461  case 0: numEle = in->getNumPoints(); break;
2462  case 1: numEle = in->getNumLines(); break;
2463  case 3: numEle = in->getNumTriangles(); break;
2464  case 4: numEle = in->getNumQuadrangles(); break;
2465  case 6: numEle = in->getNumTetrahedra(); break;
2466  case 9: numEle = in->getNumPrisms(); break;
2467  case 8: numEle = in->getNumPyramids(); break;
2468  case 12: numEle = in->getNumHexahedra(); break;
2469  }
2470  if(!numEle) return;
2471 
2472  // New variables for high order visualiztion through vtk files
2473  int numNodInsert = 0;
2474  nodMap<T> myNodMap;
2475 
2476  double minVal = in->getMin(step);
2477  double maxVal = in->getMax(step);
2478 
2479  for(int ent = 0; ent < in->getNumEntities(step); ent++) {
2480  for(int ele = 0; ele < in->getNumElements(step, ent); ele++) {
2481  if(in->skipElement(step, ent, ele) ||
2482  in->getNumEdges(step, ent, ele) != T::numEdges)
2483  continue;
2484  int numNodes = in->getNumNodes(step, ent, ele);
2485  std::vector<PCoords> coords;
2486  for(int i = 0; i < numNodes; i++) {
2487  double x, y, z;
2488  in->getNode(step, ent, ele, i, x, y, z);
2489  coords.push_back(PCoords(x, y, z));
2490  }
2491  int numVal = in->getNumValues(step, ent, ele);
2492  std::vector<PValues> values;
2493 
2494  switch(numComp) {
2495  case 1:
2496  for(int i = 0; i < numVal; i++) {
2497  double val;
2498  in->getValue(step, ent, ele, i, val);
2499  values.push_back(PValues(val));
2500  }
2501  break;
2502  case 3:
2503  for(int i = 0; i < numVal / 3; i++) {
2504  double vx, vy, vz;
2505  in->getValue(step, ent, ele, 3 * i, vx);
2506  in->getValue(step, ent, ele, 3 * i + 1, vy);
2507  in->getValue(step, ent, ele, 3 * i + 2, vz);
2508  values.push_back(PValues(vx, vy, vz));
2509  }
2510  break;
2511  case 9:
2512  for(int i = 0; i < numVal / 9; i++) {
2513  double vxx, vxy, vxz, vyx, vyy, vyz, vzx, vzy, vzz;
2514  in->getValue(step, ent, ele, 9 * i + 0, vxx);
2515  in->getValue(step, ent, ele, 9 * i + 1, vxy);
2516  in->getValue(step, ent, ele, 9 * i + 2, vxz);
2517  in->getValue(step, ent, ele, 9 * i + 3, vyx);
2518  in->getValue(step, ent, ele, 9 * i + 4, vyy);
2519  in->getValue(step, ent, ele, 9 * i + 5, vyz);
2520  in->getValue(step, ent, ele, 9 * i + 6, vzx);
2521  in->getValue(step, ent, ele, 9 * i + 7, vzy);
2522  in->getValue(step, ent, ele, 9 * i + 8, vzz);
2523  values.push_back(
2524  PValues(vxx, vxy, vxz, vyx, vyy, vyz, vzx, vzy, vzz));
2525  }
2526  break;
2527  }
2528 
2529  adaptForVTK(myVTKData.vtkTol, numComp, coords, values, minVal,
2530  maxVal); // ,plug);
2531 
2532  // Inside initial element, after adapt() has been called
2533 
2534  // Build the mapping of the canonical element,
2535  // or recycle existing one in case of uniform refinement
2536  buildMapping(myNodMap, myVTKData.vtkTol, numNodInsert);
2537 
2538  // Pre-allocate some space for the local coordinates and connectivity
2539  // in order to write to any component of the vector later through vec[i]
2540 
2541  myVTKData.vtkLocalCoords.resize(numNodInsert, PCoords(0.0, 0.0, 0.0));
2542  myVTKData.vtkLocalValues.resize(numNodInsert, PValues(numComp));
2543 
2544  for(std::size_t i = 0; i < coords.size() / T::numNodes; i++) {
2545  // Loop over
2546  // - all refined elements if refinement level > 0
2547  // - single initial element when refinement box is checked for the
2548  // first time (ref level =0)
2549 
2550  // local connectivity for the considered sub triangle
2551  vectInt vtkElmConnectivity;
2552 
2553  for(int k = 0; k < T::numNodes; ++k) {
2554  // Connectivity of the considered sub-element
2555  int countTotNodloc = T::numNodes * i + k; // Nodes are duplicate here
2556  int vtkNodeId =
2557  myVTKData.vtkCountTotNod + myNodMap.mapping[countTotNodloc];
2558  vtkElmConnectivity.push_back(vtkNodeId);
2559 
2560  // Coordinates of the nodes of the considered sub-element
2561  double px, py, pz;
2562  px = coords[T::numNodes * i + k].c[0];
2563  py = coords[T::numNodes * i + k].c[1];
2564  pz = coords[T::numNodes * i + k].c[2];
2565  PCoords tmpCoords = PCoords(px, py, pz);
2566  myVTKData.vtkLocalCoords[myNodMap.mapping[countTotNodloc]] =
2567  tmpCoords;
2568 
2569  // Value associated with each nodes of the sub-element
2570  myVTKData.vtkLocalValues[myNodMap.mapping[countTotNodloc]] =
2571  values[T::numNodes * i + k];
2572  }
2573 
2574  // Add elm connectivity to vector
2575  myVTKData.vtkLocalConnectivity.push_back(vtkElmConnectivity);
2576 
2577  // Increment global elm number
2578  myVTKData.incrementTotElm(1);
2579 
2580  // Save element type
2581  myVTKData.vtkLocalCellType.push_back(
2582  myVTKData.getPVCellType(T::numEdges));
2583 
2584  // Global variables
2585  if(buildStaticData == true) {
2586  globalVTKData::vtkGlobalConnectivity.push_back(vtkElmConnectivity);
2587  globalVTKData::vtkGlobalCellType.push_back(
2588  myVTKData.getPVCellType(T::numEdges));
2589  }
2590 
2591  // Clear existing structure (safer)
2592  vtkElmConnectivity.clear();
2593  }
2594 
2595  // Increment global node and elm-lev0 number
2596  myVTKData.incrementTotNod(numNodInsert);
2597  myVTKData.incrementTotElmLev0(1);
2598 
2599  // Write the VTK data structure of the consider element to vtu file
2600 
2601  if(writeVTK == true) { myVTKData.writeVTKElmData(); }
2602 
2603  if(buildStaticData == true) {
2604  for(int i = 0; i < numNodInsert; i++) {
2605  globalVTKData::vtkGlobalCoords.push_back(myVTKData.vtkLocalCoords[i]);
2606  }
2607 
2608  for(int i = 0; i < numNodInsert; i++) {
2609  globalVTKData::vtkGlobalValues.push_back(myVTKData.vtkLocalValues[i]);
2610  }
2611  }
2612 
2613  myVTKData.clearLocalData();
2614 
2615  } // loop over mesh element
2616  }
2617 }
2618 
2619 template <class T>
2620 int adaptiveElements<T>::countElmLev0(int step, PViewData *in)
2621 {
2622  int sum = 0;
2623  for(int ent = 0; ent < in->getNumEntities(step); ent++) {
2624  for(int ele = 0; ele < in->getNumElements(step, ent); ele++) {
2625  if(in->skipElement(step, ent, ele) ||
2626  in->getNumEdges(step, ent, ele) != T::numEdges)
2627  continue;
2628  else
2629  sum++;
2630  }
2631  }
2632  return sum;
2633 }
2634 
2635 int adaptiveData::countTotElmLev0(int step, PViewData *in)
2636 {
2637  int sumElm = 0;
2638 
2639  if(_triangles) sumElm += _triangles->countElmLev0(step, in);
2640  if(_quadrangles) sumElm += _quadrangles->countElmLev0(step, in);
2641  if(_tetrahedra) sumElm += _tetrahedra->countElmLev0(step, in);
2642  if(_prisms) sumElm += _prisms->countElmLev0(step, in);
2643  if(_hexahedra) sumElm += _hexahedra->countElmLev0(step, in);
2644  if(_pyramids) sumElm += _pyramids->countElmLev0(step, in);
2645 
2646  return sumElm;
2647 }
2648 
2649 void adaptiveData::changeResolutionForVTK(int step, int level, double tol,
2650  int npart, bool isBinary,
2651  const std::string &guiFileName,
2652  int useDefaultName)
2653 {
2654  // clean global VTK data structure before (re)generating it
2655  if(buildStaticData == true) globalVTKData::clearGlobalData();
2656 
2657  VTKData myVTKData(_inData->getName(), _inData->getNumComponents(0, 0, 0),
2658  step, level, tol, guiFileName, useDefaultName, npart,
2659  isBinary);
2660  myVTKData.vtkTotNumElmLev0 = countTotElmLev0(step, _inData);
2661  myVTKData.setFileDistribution();
2662 
2663  // Views of 2D and 3D elements only supported for VTK. _points and _lines are
2664  // currently ignored.
2665  if(_triangles) _triangles->init(myVTKData.vtkLevel);
2666  if(_quadrangles) _quadrangles->init(myVTKData.vtkLevel);
2667  if(_tetrahedra) _tetrahedra->init(myVTKData.vtkLevel);
2668  if(_prisms) _prisms->init(myVTKData.vtkLevel);
2669  if(_hexahedra) _hexahedra->init(myVTKData.vtkLevel);
2670  if(_pyramids) _pyramids->init(myVTKData.vtkLevel);
2671 
2672  if(_triangles)
2673  _triangles->addInViewForVTK(step, _inData, myVTKData, writeVTK,
2674  buildStaticData);
2675  if(_quadrangles)
2676  _quadrangles->addInViewForVTK(step, _inData, myVTKData, writeVTK,
2677  buildStaticData);
2678  if(_tetrahedra)
2679  _tetrahedra->addInViewForVTK(step, _inData, myVTKData, writeVTK,
2680  buildStaticData);
2681  if(_prisms)
2682  _prisms->addInViewForVTK(step, _inData, myVTKData, writeVTK,
2683  buildStaticData);
2684  if(_hexahedra)
2685  _hexahedra->addInViewForVTK(step, _inData, myVTKData, writeVTK,
2686  buildStaticData);
2687  if(_pyramids)
2688  _pyramids->addInViewForVTK(step, _inData, myVTKData, writeVTK,
2689  buildStaticData);
2690 
2691  Msg::StatusBar(true, "Done writing VTK data");
2692 }
adaptiveData::_triangles
adaptiveElements< adaptiveTriangle > * _triangles
Definition: adaptiveData.h:648
PViewData::getNumEdges
virtual int getNumEdges(int step, int ent, int ele)
Definition: PViewData.h:180
PViewDataList::VQ
std::vector< double > VQ
Definition: PViewDataList.h:33
PViewDataList::TS
std::vector< double > TS
Definition: PViewDataList.h:37
adaptiveVertex::valzx
double valzx
Definition: adaptiveData.h:52
adaptiveLine::all
static std::list< adaptiveLine * > all
Definition: adaptiveData.h:109
TimeOfDay
double TimeOfDay()
Definition: OS.cpp:399
adaptiveElements::_coeffsVal
fullMatrix< double > * _coeffsVal
Definition: adaptiveData.h:598
PViewData::getNumLines
virtual int getNumLines(int step=-1)
Definition: PViewData.h:115
adaptiveData::~adaptiveData
~adaptiveData()
Definition: adaptiveData.cpp:1548
adaptiveData::_tol
double _tol
Definition: adaptiveData.h:643
adaptivePyramid::create
static void create(int maxlevel)
Definition: adaptiveData.cpp:932
adaptiveQuadrangle::allVertices
static std::set< adaptiveVertex > allVertices
Definition: adaptiveData.h:169
PViewDataList::NbVY
int NbVY
Definition: PViewDataList.h:42
adaptiveLine::adaptiveLine
adaptiveLine(adaptiveVertex *p1, adaptiveVertex *p2)
Definition: adaptiveData.h:114
VTKData::initVTKFile
void initVTKFile()
Definition: adaptiveData.cpp:1800
PViewDataList::TQ
std::vector< double > TQ
Definition: PViewDataList.h:33
PViewData::getNumPyramids
virtual int getNumPyramids(int step=-1)
Definition: PViewData.h:122
adaptiveData::countTotElmLev0
int countTotElmLev0(int step, PViewData *in)
Definition: adaptiveData.cpp:2635
PViewDataList::Min
double Min
Definition: PViewDataList.h:22
adaptivePyramid::numEdges
static int numEdges
Definition: adaptiveData.h:340
adaptiveHexahedron::allVertices
static std::set< adaptiveVertex > allVertices
Definition: adaptiveData.h:288
std::swap
void swap(picojson::value &x, picojson::value &y)
Definition: picojson.h:1136
adaptivePyramid::p
adaptiveVertex * p[5]
Definition: adaptiveData.h:336
PViewData::skipElement
virtual bool skipElement(int step, int ent, int ele, bool checkVisibility=false, int samplingRate=1)
Definition: PViewData.cpp:90
adaptivePyramid
Definition: adaptiveData.h:333
VTKData::vtkCountTotElm
int vtkCountTotElm
Definition: adaptiveData.h:516
adaptiveTriangle::V
double V() const
Definition: adaptiveData.h:150
Plugin.h
distance
double distance(MVertex *v1, MVertex *v2)
Definition: MVertex.h:245
adaptiveHexahedron::recurCreate
static void recurCreate(adaptiveHexahedron *h, int maxlevel, int level)
Definition: adaptiveData.cpp:611
adaptiveElements::countElmLev0
int countElmLev0(int step, PViewData *in)
Definition: adaptiveData.cpp:2620
VTKData::getPVCellType
int getPVCellType(int numEdges)
Definition: adaptiveData.cpp:2198
PViewDataList::NbTL
int NbTL
Definition: PViewDataList.h:28
VTKData::setFileDistribution
void setFileDistribution()
Definition: adaptiveData.h:580
PViewDataList::SS
std::vector< double > SS
Definition: PViewDataList.h:37
adaptivePyramid::allVertices
static std::set< adaptiveVertex > allVertices
Definition: adaptiveData.h:339
fullVector< double >
TYPE_LIN
#define TYPE_LIN
Definition: GmshDefines.h:65
adaptivePoint::adaptivePoint
adaptivePoint(adaptiveVertex *p1)
Definition: adaptiveData.h:88
PViewDataList::TL
std::vector< double > TL
Definition: PViewDataList.h:29
PCoords
Definition: adaptiveData.h:374
adaptiveTetrahedron::numNodes
static int numNodes
Definition: adaptiveData.h:251
adaptivePrism::create
static void create(int maxlevel)
Definition: adaptiveData.cpp:793
adaptivePyramid::all
static std::list< adaptivePyramid * > all
Definition: adaptiveData.h:338
adaptiveElements::adaptiveElements
adaptiveElements(std::vector< fullMatrix< double > * > &interpolationMatrices)
Definition: adaptiveData.cpp:1082
OS.h
adaptiveHexahedron::adaptiveHexahedron
adaptiveHexahedron(adaptiveVertex *p1, adaptiveVertex *p2, adaptiveVertex *p3, adaptiveVertex *p4, adaptiveVertex *p5, adaptiveVertex *p6, adaptiveVertex *p7, adaptiveVertex *p8)
Definition: adaptiveData.h:292
PViewDataList::NbTP
int NbTP
Definition: PViewDataList.h:26
PViewDataList::NbTT
int NbTT
Definition: PViewDataList.h:30
PViewDataList::Max
double Max
Definition: PViewDataList.h:22
adaptivePrism::numEdges
static int numEdges
Definition: adaptiveData.h:207
adaptiveHexahedron::e
adaptiveHexahedron * e[8]
Definition: adaptiveData.h:286
PViewDataList
Definition: PViewDataList.h:17
adaptiveTriangle::all
static std::list< adaptiveTriangle * > all
Definition: adaptiveData.h:137
Msg::Warning
static void Warning(const char *fmt,...)
Definition: GmshMessage.cpp:543
PViewData::getNode
virtual int getNode(int step, int ent, int ele, int nod, double &x, double &y, double &z)
Definition: PViewData.h:141
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
PViewDataList::NbSP
int NbSP
Definition: PViewDataList.h:26
Msg::StatusBar
static void StatusBar(bool log, const char *fmt,...)
Definition: GmshMessage.cpp:686
PViewData::getInterpolationMatrices
int getInterpolationMatrices(int type, std::vector< fullMatrix< double > * > &p)
Definition: PViewData.cpp:176
adaptiveElements::_eexpsGeom
fullMatrix< double > * _eexpsGeom
Definition: adaptiveData.h:599
adaptivePoint::allVertices
static std::set< adaptiveVertex > allVertices
Definition: adaptiveData.h:84
VTKData::vtkIsBinary
bool vtkIsBinary
Definition: adaptiveData.h:500
PViewData::getValue
virtual void getValue(int step, int ent, int ele, int idx, double &val)
Definition: PViewData.h:159
PViewDataList::TI
std::vector< double > TI
Definition: PViewDataList.h:41
PViewDataList::NbSQ
int NbSQ
Definition: PViewDataList.h:32
TYPE_PNT
#define TYPE_PNT
Definition: GmshDefines.h:64
adaptivePoint::recurCreate
static void recurCreate(adaptivePoint *e, int maxlevel, int level)
Definition: adaptiveData.cpp:133
VTKData::finalizeVTKFile
void finalizeVTKFile()
Definition: adaptiveData.cpp:1874
PViewDataList::SY
std::vector< double > SY
Definition: PViewDataList.h:43
adaptiveVertex::valzz
double valzz
Definition: adaptiveData.h:52
PViewData::getNumPoints
virtual int getNumPoints(int step=-1)
Definition: PViewData.h:114
PViewDataList::NbSH
int NbSH
Definition: PViewDataList.h:38
TYPE_TRI
#define TYPE_TRI
Definition: GmshDefines.h:66
adaptiveVertex::val
double val
Definition: adaptiveData.h:50
adaptiveTriangle::numNodes
static int numNodes
Definition: adaptiveData.h:139
PViewDataList::VL
std::vector< double > VL
Definition: PViewDataList.h:29
adaptiveElements::buildMapping
void buildMapping(nodMap< T > &myNodMap, double tol, int &numNodInsert)
Definition: adaptiveData.cpp:2381
PViewData::getNumValues
virtual int getNumValues(int step, int ent, int ele)
Definition: PViewData.h:156
adaptivePyramid::numNodes
static int numNodes
Definition: adaptiveData.h:340
adaptiveData::_step
int _step
Definition: adaptiveData.h:642
globalVTKData::vtkGlobalValues
static std::vector< PValues > vtkGlobalValues
Definition: adaptiveData.h:449
adaptiveTetrahedron
Definition: adaptiveData.h:244
ToString
std::string ToString(const T &val)
Definition: adaptiveData.h:39
globalVTKData::vtkGlobalCellType
static std::vector< int > vtkGlobalCellType
Definition: adaptiveData.h:446
adaptivePrism::error
static void error(double AVG, double tol)
Definition: adaptiveData.cpp:874
adaptivePoint::all
static std::list< adaptivePoint * > all
Definition: adaptiveData.h:83
adaptiveQuadrangle
Definition: adaptiveData.h:163
adaptiveHexahedron::p
adaptiveVertex * p[8]
Definition: adaptiveData.h:285
adaptiveHexahedron::create
static void create(int maxlevel)
Definition: adaptiveData.cpp:595
PViewDataList::VY
std::vector< double > VY
Definition: PViewDataList.h:43
adaptiveQuadrangle::recurError
static void recurError(adaptiveQuadrangle *q, double AVG, double tol)
Definition: adaptiveData.cpp:395
adaptiveData::changeResolution
void changeResolution(int step, int level, double tol, GMSH_PostPlugin *plug=nullptr)
Definition: adaptiveData.cpp:1564
adaptiveHexahedron::numEdges
static int numEdges
Definition: adaptiveData.h:289
adaptivePoint::error
static void error(double AVG, double tol)
Definition: adaptiveData.cpp:138
PViewData::getNumEntities
virtual int getNumEntities(int step=-1)
Definition: PViewData.h:127
computeShapeFunctionsPyramid
static void computeShapeFunctionsPyramid(fullMatrix< double > *coeffs, fullMatrix< double > *eexps, double u, double v, double w, fullVector< double > *sf, fullVector< double > *tmp)
Definition: adaptiveData.cpp:90
TYPE_PRI
#define TYPE_PRI
Definition: GmshDefines.h:70
adaptiveLine::create
static void create(int maxlevel)
Definition: adaptiveData.cpp:149
PViewDataList::NbVH
int NbVH
Definition: PViewDataList.h:38
adaptivePoint::create
static void create(int maxlevel)
Definition: adaptiveData.cpp:125
adaptiveVertex::y
float y
Definition: adaptiveData.h:48
adaptiveTetrahedron::create
static void create(int maxlevel)
Definition: adaptiveData.cpp:459
adaptiveVertex::valyz
double valyz
Definition: adaptiveData.h:51
adaptiveTriangle::allVertices
static std::set< adaptiveVertex > allVertices
Definition: adaptiveData.h:138
PViewDataList::VH
std::vector< double > VH
Definition: PViewDataList.h:39
adaptiveHexahedron
Definition: adaptiveData.h:282
VTKData::vtkFormat
std::string vtkFormat
Definition: adaptiveData.h:491
adaptiveData::changeResolutionForVTK
void changeResolutionForVTK(int step, int level, double tol, int npart=1, bool isBinary=true, const std::string &guifileName="unknown", int useDefaultName=1)
Definition: adaptiveData.cpp:2649
adaptiveQuadrangle::p
adaptiveVertex * p[4]
Definition: adaptiveData.h:166
adaptiveTriangle::adaptiveTriangle
adaptiveTriangle(adaptiveVertex *p1, adaptiveVertex *p2, adaptiveVertex *p3)
Definition: adaptiveData.h:142
adaptiveTriangle::create
static void create(int maxlevel)
Definition: adaptiveData.cpp:223
adaptiveHexahedron::recurError
static void recurError(adaptiveHexahedron *h, double AVG, double tol)
Definition: adaptiveData.cpp:725
VTKData::vtkFileCellType
FILE * vtkFileCellType
Definition: adaptiveData.h:509
nodMap::getSize
int getSize()
Definition: adaptiveData.h:75
adaptiveLine::V
double V() const
Definition: adaptiveData.h:120
adaptiveData.h
VTKData::vtkUseDefaultName
int vtkUseDefaultName
Definition: adaptiveData.h:501
globalVTKData::clearGlobalData
static void clearGlobalData()
Definition: adaptiveData.h:475
adaptiveElements::adapt
bool adapt(double tol, int numComp, std::vector< PCoords > &coords, std::vector< PValues > &values, double &minVal, double &maxVal, GMSH_PostPlugin *plug=nullptr, bool onlyComputeMinMax=false)
Definition: adaptiveData.cpp:1206
fullMatrix< double >
PViewDataList::TY
std::vector< double > TY
Definition: PViewDataList.h:43
adaptiveData::timerInit
static double timerInit
Definition: adaptiveData.h:671
adaptiveTetrahedron::visible
bool visible
Definition: adaptiveData.h:246
PViewDataList::TP
std::vector< double > TP
Definition: PViewDataList.h:27
adaptiveLine::allVertices
static std::set< adaptiveVertex > allVertices
Definition: adaptiveData.h:110
adaptiveElements< adaptivePoint >
adaptiveHexahedron::visible
bool visible
Definition: adaptiveData.h:284
PViewData::getMax
virtual double getMax(int step=-1, bool onlyVisible=false, int tensorRep=0, int forceNumComponents=0, int componentMap[9]=nullptr)=0
VTKData::vtkCountTotNodConnect
int vtkCountTotNodConnect
Definition: adaptiveData.h:518
VTKData::vtkCountCellOffset
int vtkCountCellOffset
Definition: adaptiveData.h:520
VTKData::vtkFieldName
std::string vtkFieldName
Definition: adaptiveData.h:489
adaptiveTriangle::error
static void error(double AVG, double tol)
Definition: adaptiveData.cpp:267
adaptiveVertex::valz
double valz
maximal three values
Definition: adaptiveData.h:50
PViewData::getNumHexahedra
virtual int getNumHexahedra(int step=-1)
Definition: PViewData.h:120
pe2
const double pe2
Definition: GaussQuadratureQuad.cpp:76
adaptiveQuadrangle::error
static void error(double AVG, double tol)
Definition: adaptiveData.cpp:389
adaptiveQuadrangle::V
double V() const
Definition: adaptiveData.h:183
PViewDataList::VP
std::vector< double > VP
Definition: PViewDataList.h:27
adaptiveLine::p
adaptiveVertex * p[2]
Definition: adaptiveData.h:107
PViewDataList::NbTH
int NbTH
Definition: PViewDataList.h:38
CreateSingleDir
int CreateSingleDir(const std::string &dirName)
Definition: OS.cpp:502
PViewDataList::SP
std::vector< double > SP
Definition: PViewDataList.h:27
VTKData::clearLocalData
void clearLocalData()
Definition: adaptiveData.h:557
adaptivePrism::adaptivePrism
adaptivePrism(adaptiveVertex *p1, adaptiveVertex *p2, adaptiveVertex *p3, adaptiveVertex *p4, adaptiveVertex *p5, adaptiveVertex *p6)
Definition: adaptiveData.h:210
adaptiveTriangle::recurCreate
static void recurCreate(adaptiveTriangle *t, int maxlevel, int level)
Definition: adaptiveData.cpp:233
VTKData::vtkCountCellType
int vtkCountCellType
Definition: adaptiveData.h:521
adaptivePyramid::GSF
static void GSF(double u, double v, double w, fullVector< double > &sf)
Definition: adaptiveData.h:359
adaptiveElements::~adaptiveElements
~adaptiveElements()
Definition: adaptiveData.cpp:1096
VTKData::isLittleEndian
bool isLittleEndian()
Definition: adaptiveData.cpp:1602
VTKData::vtkCountTotNod
int vtkCountTotNod
Definition: adaptiveData.h:515
adaptiveTetrahedron::V
double V() const
Definition: adaptiveData.h:265
adaptiveTetrahedron::numEdges
static int numEdges
Definition: adaptiveData.h:251
adaptiveData::buildStaticData
bool buildStaticData
Definition: adaptiveData.h:660
adaptiveData::_lines
adaptiveElements< adaptiveLine > * _lines
Definition: adaptiveData.h:647
PViewData::getNumQuadrangles
virtual int getNumQuadrangles(int step=-1)
Definition: PViewData.h:117
PViewDataList::SL
std::vector< double > SL
Definition: PViewDataList.h:29
VTKData::vtkNpart
int vtkNpart
Definition: adaptiveData.h:498
adaptiveData::_outData
PViewDataList * _outData
Definition: adaptiveData.h:645
adaptiveQuadrangle::recurCreate
static void recurCreate(adaptiveQuadrangle *q, int maxlevel, int level)
Definition: adaptiveData.cpp:346
PViewDataList::NbVS
int NbVS
Definition: PViewDataList.h:36
VTKData::vtkFileCoord
FILE * vtkFileCoord
Definition: adaptiveData.h:506
VTKData::vtkFileCellOffset
FILE * vtkFileCellOffset
Definition: adaptiveData.h:508
VTKData::vtkLocalConnectivity
std::vector< vectInt > vtkLocalConnectivity
Definition: adaptiveData.h:523
adaptiveData::_tetrahedra
adaptiveElements< adaptiveTetrahedron > * _tetrahedra
Definition: adaptiveData.h:650
adaptiveTriangle::recurError
static void recurError(adaptiveTriangle *t, double AVG, double tol)
Definition: adaptiveData.cpp:273
PViewData::getNumNodes
virtual int getNumNodes(int step, int ent, int ele)
Definition: PViewData.h:137
GmshDefines.h
adaptiveQuadrangle::all
static std::list< adaptiveQuadrangle * > all
Definition: adaptiveData.h:168
adaptiveElements::_eexpsVal
fullMatrix< double > * _eexpsVal
Definition: adaptiveData.h:598
adaptiveTetrahedron::error
static void error(double AVG, double tol)
Definition: adaptiveData.cpp:524
adaptiveVertex::z
float z
parametric coordinates
Definition: adaptiveData.h:48
adaptiveTriangle
Definition: adaptiveData.h:132
VTKData
Definition: adaptiveData.h:485
adaptiveTetrahedron::allVertices
static std::set< adaptiveVertex > allVertices
Definition: adaptiveData.h:250
VTKData::vtkTol
double vtkTol
Definition: adaptiveData.h:497
adaptivePoint::e
adaptivePoint * e[1]
Definition: adaptiveData.h:82
adaptiveData::timerAdapt
static double timerAdapt
Definition: adaptiveData.h:671
adaptiveVertex::valzy
double valzy
Definition: adaptiveData.h:52
adaptiveElements::_interpolVal
fullMatrix< double > * _interpolVal
Definition: adaptiveData.h:598
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
adaptivePyramid::recurError
static void recurError(adaptivePyramid *h, double AVG, double tol)
Definition: adaptiveData.cpp:1038
PViewDataGModel.h
adaptiveHexahedron::numNodes
static int numNodes
Definition: adaptiveData.h:289
adaptivePrism::p
adaptiveVertex * p[6]
Definition: adaptiveData.h:203
adaptiveQuadrangle::visible
bool visible
Definition: adaptiveData.h:165
VTKData::SwapArrayByteOrder
void SwapArrayByteOrder(void *array, int nbytes, int nItems)
Definition: adaptiveData.cpp:1611
adaptiveLine::visible
bool visible
Definition: adaptiveData.h:106
PViewDataList::ST
std::vector< double > ST
Definition: PViewDataList.h:31
adaptivePoint::numNodes
static int numNodes
Definition: adaptiveData.h:85
adaptiveTetrahedron::all
static std::list< adaptiveTetrahedron * > all
Definition: adaptiveData.h:249
adaptiveVertex::Z
double Z
cartesian coordinates
Definition: adaptiveData.h:49
PViewDataList::TT
std::vector< double > TT
Definition: PViewDataList.h:31
adaptiveLine::recurError
static void recurError(adaptiveLine *e, double AVG, double tol)
Definition: adaptiveData.cpp:183
GMSH_PostPlugin::assignSpecificVisibility
virtual void assignSpecificVisibility() const
Definition: Plugin.h:111
VTKData::vtkLevel
int vtkLevel
Definition: adaptiveData.h:495
VTKData::vtkDirName
std::string vtkDirName
Definition: adaptiveData.h:492
PViewDataList::NbVP
int NbVP
Definition: PViewDataList.h:26
PViewDataList::NbTS
int NbTS
Definition: PViewDataList.h:36
PViewData::setDirty
virtual void setDirty(bool val)
Definition: PViewData.h:63
adaptiveLine::e
adaptiveLine * e[2]
Definition: adaptiveData.h:108
adaptiveTriangle::p
adaptiveVertex * p[3]
Definition: adaptiveData.h:135
fullMatrix::size2
int size2() const
Definition: fullMatrix.h:275
VTKData::maxElmPerPart
int maxElmPerPart
Definition: adaptiveData.h:502
adaptivePoint
Definition: adaptiveData.h:78
adaptiveData::_level
int _level
Definition: adaptiveData.h:642
PViewData
Definition: PViewData.h:29
VTKData::vtkLocalCoords
std::vector< PCoords > vtkLocalCoords
Definition: adaptiveData.h:525
adaptiveLine::error
static void error(double AVG, double tol)
Definition: adaptiveData.cpp:177
TYPE_PYR
#define TYPE_PYR
Definition: GmshDefines.h:69
adaptiveVertex::X
double X
Definition: adaptiveData.h:49
adaptiveQuadrangle::numNodes
static int numNodes
Definition: adaptiveData.h:170
adaptivePrism::recurCreate
static void recurCreate(adaptivePrism *p, int maxlevel, int level)
Definition: adaptiveData.cpp:806
PViewData::getNumTetrahedra
virtual int getNumTetrahedra(int step=-1)
Definition: PViewData.h:119
adaptivePrism::recurError
static void recurError(adaptivePrism *p, double AVG, double tol)
Definition: adaptiveData.cpp:880
nodMap::mapping
std::vector< int > mapping
Definition: adaptiveData.h:70
adaptiveHexahedron::V
double V() const
Definition: adaptiveData.h:308
adaptiveData::_points
adaptiveElements< adaptivePoint > * _points
Definition: adaptiveData.h:646
adaptiveData::_pyramids
adaptiveElements< adaptivePyramid > * _pyramids
Definition: adaptiveData.h:653
adaptiveVertex::x
float x
Definition: adaptiveData.h:48
adaptiveLine::recurCreate
static void recurCreate(adaptiveLine *e, int maxlevel, int level)
Definition: adaptiveData.cpp:158
VTKData::vtkFileNodVal
FILE * vtkFileNodVal
Definition: adaptiveData.h:510
PViewDataList::NbSL
int NbSL
Definition: PViewDataList.h:28
adaptiveVertex::add
static adaptiveVertex * add(double x, double y, double z, std::set< adaptiveVertex > &allVertice)
Definition: adaptiveData.cpp:110
TYPE_QUA
#define TYPE_QUA
Definition: GmshDefines.h:67
adaptiveQuadrangle::adaptiveQuadrangle
adaptiveQuadrangle(adaptiveVertex *p1, adaptiveVertex *p2, adaptiveVertex *p3, adaptiveVertex *p4)
Definition: adaptiveData.h:173
PViewDataList::TH
std::vector< double > TH
Definition: PViewDataList.h:39
VTKData::incrementTotElmLev0
void incrementTotElmLev0(int increment)
Definition: adaptiveData.h:571
PViewDataList::SH
std::vector< double > SH
Definition: PViewDataList.h:39
PViewDataList::NbTI
int NbTI
Definition: PViewDataList.h:40
adaptivePrism::allVertices
static std::set< adaptiveVertex > allVertices
Definition: adaptiveData.h:206
adaptiveVertex::valyy
double valyy
Definition: adaptiveData.h:51
adaptiveData::_hexahedra
adaptiveElements< adaptiveHexahedron > * _hexahedra
Definition: adaptiveData.h:651
adaptiveData::upBuildStaticData
void upBuildStaticData(bool newValue)
Definition: adaptiveData.h:682
adaptiveTetrahedron::p
adaptiveVertex * p[4]
Definition: adaptiveData.h:247
PViewDataList::NbSS
int NbSS
Definition: PViewDataList.h:36
adaptiveElements::adaptForVTK
void adaptForVTK(double tol, int numComp, std::vector< PCoords > &coords, std::vector< PValues > &values, double &minVal, double &maxVal)
Definition: adaptiveData.cpp:2240
PViewData::getNumComponents
virtual int getNumComponents(int step, int ent, int ele)
Definition: PViewData.h:152
adaptiveTetrahedron::recurCreate
static void recurCreate(adaptiveTetrahedron *t, int maxlevel, int level)
Definition: adaptiveData.cpp:470
PViewDataList::NbTY
int NbTY
Definition: PViewDataList.h:42
PViewDataList::NbVL
int NbVL
Definition: PViewDataList.h:28
fullMatrix::size1
int size1() const
Definition: fullMatrix.h:274
z
const double z
Definition: GaussQuadratureQuad.cpp:56
adaptiveQuadrangle::create
static void create(int maxlevel)
Definition: adaptiveData.cpp:335
adaptiveTriangle::visible
bool visible
Definition: adaptiveData.h:134
GMSH_PostPlugin
Definition: Plugin.h:83
adaptiveTriangle::numEdges
static int numEdges
Definition: adaptiveData.h:139
VTKData::writeVTKElmData
void writeVTKElmData()
Definition: adaptiveData.cpp:1624
adaptiveTetrahedron::adaptiveTetrahedron
adaptiveTetrahedron(adaptiveVertex *p1, adaptiveVertex *p2, adaptiveVertex *p3, adaptiveVertex *p4)
Definition: adaptiveData.h:254
PViewData::getNumPrisms
virtual int getNumPrisms(int step=-1)
Definition: PViewData.h:121
adaptiveData::_prisms
adaptiveElements< adaptivePrism > * _prisms
Definition: adaptiveData.h:652
VTKData::vtkStep
int vtkStep
Definition: adaptiveData.h:494
adaptiveTetrahedron::e
adaptiveTetrahedron * e[8]
Definition: adaptiveData.h:248
PViewDataList::NbSI
int NbSI
Definition: PViewDataList.h:40
VTKData::vtkNumComp
int vtkNumComp
Definition: adaptiveData.h:496
VTKData::minElmPerPart
int minElmPerPart
Definition: adaptiveData.h:502
picojson::array
value::array array
Definition: picojson.h:194
nodMap::cleanMapping
void cleanMapping()
Definition: adaptiveData.h:73
VTKData::vtkCountFile
int vtkCountFile
Definition: adaptiveData.h:511
adaptiveElements::addInViewForVTK
void addInViewForVTK(int step, PViewData *in, VTKData &myVTKData, bool writeVtk=true, bool buildStaticData=false)
Definition: adaptiveData.cpp:2452
PViewDataList::VS
std::vector< double > VS
Definition: PViewDataList.h:37
adaptiveData::adaptiveData
adaptiveData(PViewData *data, bool outDataInit=true)
Definition: adaptiveData.cpp:1496
TYPE_HEX
#define TYPE_HEX
Definition: GmshDefines.h:71
VTKData::vtkTotNumElmLev0
int vtkTotNumElmLev0
Definition: adaptiveData.h:513
adaptiveData::_quadrangles
adaptiveElements< adaptiveQuadrangle > * _quadrangles
Definition: adaptiveData.h:649
VTKData::vtkCountTotElmLev0
int vtkCountTotElmLev0
Definition: adaptiveData.h:514
TYPE_TET
#define TYPE_TET
Definition: GmshDefines.h:68
adaptivePyramid::adaptivePyramid
adaptivePyramid(adaptiveVertex *p1, adaptiveVertex *p2, adaptiveVertex *p3, adaptiveVertex *p4, adaptiveVertex *p5)
Definition: adaptiveData.h:343
XYZ
Definition: Camera.h:13
VTKData::incrementTotElm
void incrementTotElm(int increment)
Definition: adaptiveData.h:570
nodMap
Definition: adaptiveData.h:68
PViewDataList::NbTQ
int NbTQ
Definition: PViewDataList.h:32
adaptiveData::_inData
PViewData * _inData
Definition: adaptiveData.h:644
VTKData::vtkFileName
std::string vtkFileName
Definition: adaptiveData.h:490
adaptiveLine
Definition: adaptiveData.h:104
adaptivePrism::all
static std::list< adaptivePrism * > all
Definition: adaptiveData.h:205
PViewData::getMin
virtual double getMin(int step=-1, bool onlyVisible=false, int tensorRep=0, int forceNumComponents=0, int componentMap[9]=nullptr)=0
adaptivePyramid::recurCreate
static void recurCreate(adaptivePyramid *h, int maxlevel, int level)
Definition: adaptiveData.cpp:944
adaptiveHexahedron::error
static void error(double AVG, double tol)
Definition: adaptiveData.cpp:719
adaptivePoint::recurError
static void recurError(adaptivePoint *e, double AVG, double tol)
Definition: adaptiveData.cpp:144
VTKData::vtkFileConnect
FILE * vtkFileConnect
Definition: adaptiveData.h:507
VTKData::vtkCountTotVal
int vtkCountTotVal
Definition: adaptiveData.h:519
VTKData::vtkLocalCellType
std::vector< int > vtkLocalCellType
Definition: adaptiveData.h:524
PViewDataList::NbVI
int NbVI
Definition: PViewDataList.h:40
PViewDataList::SI
std::vector< double > SI
Definition: PViewDataList.h:41
PViewData::getNumElements
virtual int getNumElements(int step=-1, int ent=-1)
Definition: PViewData.h:131
computeShapeFunctions
static void computeShapeFunctions(fullMatrix< double > *coeffs, fullMatrix< double > *eexps, double u, double v, double w, fullVector< double > *sf, fullVector< double > *tmp)
Definition: adaptiveData.cpp:66
globalVTKData::vtkGlobalCoords
static std::vector< PCoords > vtkGlobalCoords
Definition: adaptiveData.h:447
adaptivePrism::numNodes
static int numNodes
Definition: adaptiveData.h:207
PViewDataList::NbVT
int NbVT
Definition: PViewDataList.h:30
fullMatrix::mult
void mult(const fullVector< scalar > &x, fullVector< scalar > &y) const
Definition: fullMatrix.h:487
PViewData::getName
virtual std::string getName()
Definition: PViewData.h:70
adaptiveVertex
Definition: adaptiveData.h:46
adaptivePyramid::error
static void error(double AVG, double tol)
Definition: adaptiveData.cpp:1032
VTKData::numPartMinElm
int numPartMinElm
Definition: adaptiveData.h:502
adaptiveData::writeVTK
bool writeVTK
Definition: adaptiveData.h:668
adaptiveElements::init
void init(int level)
Definition: adaptiveData.cpp:1103
adaptiveVertex::valy
double valy
Definition: adaptiveData.h:50
adaptiveData::upWriteVTK
void upWriteVTK(bool newValue)
Definition: adaptiveData.h:683
PViewDataList::NbVQ
int NbVQ
Definition: PViewDataList.h:32
PViewDataList::VI
std::vector< double > VI
Definition: PViewDataList.h:41
adaptiveElements::_interpolGeom
fullMatrix< double > * _interpolGeom
Definition: adaptiveData.h:599
PViewDataList::VT
std::vector< double > VT
Definition: PViewDataList.h:31
adaptiveTetrahedron::recurError
static void recurError(adaptiveTetrahedron *t, double AVG, double tol)
Definition: adaptiveData.cpp:530
adaptiveQuadrangle::numEdges
static int numEdges
Definition: adaptiveData.h:170
VTKData::vtkFile
FILE * vtkFile
Definition: adaptiveData.h:505
adaptivePoint::numEdges
static int numEdges
Definition: adaptiveData.h:85
adaptiveHexahedron::all
static std::list< adaptiveHexahedron * > all
Definition: adaptiveData.h:287
PViewDataList::SQ
std::vector< double > SQ
Definition: PViewDataList.h:33
globalVTKData::vtkGlobalConnectivity
static std::vector< vectInt > vtkGlobalConnectivity
Definition: adaptiveData.h:445
PValues
Definition: adaptiveData.h:385
adaptiveLine::numEdges
static int numEdges
Definition: adaptiveData.h:111
adaptiveElements::addInView
void addInView(double tol, int step, PViewData *in, PViewDataList *out, GMSH_PostPlugin *plug=nullptr)
Definition: adaptiveData.cpp:1356
vectInt
std::vector< int > vectInt
Definition: adaptiveData.h:32
adaptiveQuadrangle::e
adaptiveQuadrangle * e[4]
Definition: adaptiveData.h:167
PViewData::getNumTriangles
virtual int getNumTriangles(int step=-1)
Definition: PViewData.h:116
adaptiveTriangle::e
adaptiveTriangle * e[4]
Definition: adaptiveData.h:136
adaptivePoint::visible
bool visible
Definition: adaptiveData.h:80
PViewDataList::NbST
int NbST
Definition: PViewDataList.h:30
cleanElement
static void cleanElement()
Definition: adaptiveData.cpp:59
VTKData::vtkCountCoord
int vtkCountCoord
Definition: adaptiveData.h:517
PViewDataList::NbSY
int NbSY
Definition: PViewDataList.h:42
VTKData::vtkLocalValues
std::vector< PValues > vtkLocalValues
Definition: adaptiveData.h:526
adaptivePrism
Definition: adaptiveData.h:200
VTKData::incrementTotNod
void incrementTotNod(int increment)
Definition: adaptiveData.h:569
adaptiveVertex::valyx
double valyx
Definition: adaptiveData.h:51
adaptiveLine::numNodes
static int numNodes
Definition: adaptiveData.h:111
adaptiveElements::_coeffsGeom
fullMatrix< double > * _coeffsGeom
Definition: adaptiveData.h:599
adaptiveVertex::Y
double Y
Definition: adaptiveData.h:49