diff --git a/Common/AdaptiveViews.cpp b/Common/AdaptiveViews.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c77a27715ef8c07a6aeb4b1c0bfe887f11f0be04
--- /dev/null
+++ b/Common/AdaptiveViews.cpp
@@ -0,0 +1,1191 @@
+//
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include <stdio.h>
+#include <math.h>
+#include <list>
+#include <set>
+#include "AdaptiveViews.h"
+#include "Plugin.h"
+#include "OS.h"
+
+// A recursive effective implementation
+
+void computeShapeFunctions(Double_Matrix * coeffs, Double_Matrix * eexps,
+ double u, double v, double w, double *sf);
+
+std::set < adapt_point > adapt_point::all_points;
+std::list < adapt_edge * >adapt_edge::all_elems;
+std::list < adapt_triangle * >adapt_triangle::all_elems;
+std::list < adapt_tet * >adapt_tet::all_elems;
+std::list < adapt_quad * >adapt_quad::all_elems;
+std::list < adapt_hex * >adapt_hex::all_elems;
+#define MAX_NB_NOD 8
+int adapt_edge::nbNod = 2;
+int adapt_triangle::nbNod = 3;
+int adapt_tet::nbNod = 4;
+int adapt_quad::nbNod = 4;
+int adapt_hex::nbNod = 8;
+
+adapt_point *adapt_point::New(double x, double y, double z,
+ Double_Matrix * coeffs, Double_Matrix * eexps)
+{
+ adapt_point p;
+ p.x = x;
+ p.y = y;
+ p.z = z;
+ std::set < adapt_point >::iterator it = all_points.find(p);
+ if(it == all_points.end()) {
+ all_points.insert(p);
+ it = all_points.find(p);
+ double *kkk = (double *)(it->shape_functions);
+ computeShapeFunctions(coeffs, eexps, x, y, z, kkk);
+ return (adapt_point *) & (*it);
+ }
+ else
+ return (adapt_point *) & (*it);
+}
+
+void adapt_edge::Create(int maxlevel, Double_Matrix * coeffs,
+ Double_Matrix * eexps)
+{
+ int level = 0;
+ cleanElement < adapt_edge > ();
+ adapt_point *p1 = adapt_point::New(-1, 0, 0, coeffs, eexps);
+ adapt_point *p2 = adapt_point::New(1, 0, 0, coeffs, eexps);
+ adapt_edge *t = new adapt_edge(p1, p2);
+ Recur_Create(t, maxlevel, level, coeffs, eexps);
+}
+
+void adapt_triangle::Create(int maxlevel, Double_Matrix * coeffs,
+ Double_Matrix * eexps)
+{
+ int level = 0;
+ cleanElement < adapt_triangle > ();
+ adapt_point *p1 = adapt_point::New(0, 0, 0, coeffs, eexps);
+ adapt_point *p2 = adapt_point::New(0, 1, 0, coeffs, eexps);
+ adapt_point *p3 = adapt_point::New(1, 0, 0, coeffs, eexps);
+ adapt_triangle *t = new adapt_triangle(p1, p2, p3);
+ Recur_Create(t, maxlevel, level, coeffs, eexps);
+}
+
+void adapt_quad::Create(int maxlevel, Double_Matrix * coeffs,
+ Double_Matrix * eexps)
+{
+ int level = 0;
+ cleanElement < adapt_quad > ();
+ adapt_point *p1 = adapt_point::New(-1, -1, 0, coeffs, eexps);
+ adapt_point *p2 = adapt_point::New(1, -1, 0, coeffs, eexps);
+ adapt_point *p3 = adapt_point::New(1, 1, 0, coeffs, eexps);
+ adapt_point *p4 = adapt_point::New(-1, 1, 0, coeffs, eexps);
+ adapt_quad *q = new adapt_quad(p1, p2, p3, p4);
+ Recur_Create(q, maxlevel, level, coeffs, eexps);
+}
+
+void adapt_tet::Create(int maxlevel, Double_Matrix * coeffs,
+ Double_Matrix * eexps)
+{
+ int level = 0;
+ cleanElement < adapt_tet > ();
+ adapt_point *p1 = adapt_point::New(0, 0, 0, coeffs, eexps);
+ adapt_point *p2 = adapt_point::New(0, 1, 0, coeffs, eexps);
+ adapt_point *p3 = adapt_point::New(1, 0, 0, coeffs, eexps);
+ adapt_point *p4 = adapt_point::New(0, 0, 1, coeffs, eexps);
+ adapt_tet *t = new adapt_tet(p1, p2, p3, p4);
+ Recur_Create(t, maxlevel, level, coeffs, eexps);
+}
+
+void adapt_hex::Create(int maxlevel, Double_Matrix * coeffs,
+ Double_Matrix * eexps)
+{
+ int level = 0;
+ cleanElement < adapt_hex > ();
+ adapt_point *p1 = adapt_point::New(-1, -1, -1, coeffs, eexps);
+ adapt_point *p2 = adapt_point::New(-1, 1, -1, coeffs, eexps);
+ adapt_point *p3 = adapt_point::New(1, 1, -1, coeffs, eexps);
+ adapt_point *p4 = adapt_point::New(1, -1, -1, coeffs, eexps);
+ adapt_point *p11 = adapt_point::New(-1, -1, 1, coeffs, eexps);
+ adapt_point *p21 = adapt_point::New(-1, 1, 1, coeffs, eexps);
+ adapt_point *p31 = adapt_point::New(1, 1, 1, coeffs, eexps);
+ adapt_point *p41 = adapt_point::New(1, -1, 1, coeffs, eexps);
+ adapt_hex *h = new adapt_hex(p1, p2, p3, p4, p11, p21, p31, p41);
+ Recur_Create(h, maxlevel, level, coeffs, eexps);
+}
+
+void adapt_edge::Recur_Create(adapt_edge * e, int maxlevel, int level,
+ Double_Matrix * coeffs, Double_Matrix * eexps)
+{
+ all_elems.push_back(e);
+ if(level++ >= maxlevel)
+ return;
+
+ /*
+ p1 p12 p2
+ */
+
+ adapt_point *p1 = e->p[0];
+ adapt_point *p2 = e->p[1];
+ adapt_point *p12 =
+ adapt_point::New((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5,
+ (p1->z + p2->z) * 0.5, coeffs, eexps);
+ adapt_edge *e1 = new adapt_edge(p1, p12);
+ Recur_Create(e1, maxlevel, level, coeffs, eexps);
+ adapt_edge *e2 = new adapt_edge(p12, p2);
+ Recur_Create(e2, maxlevel, level, coeffs, eexps);
+ e->e[0] = e1;
+ e->e[1] = e2;
+}
+
+void adapt_triangle::Recur_Create(adapt_triangle * t, int maxlevel, int level,
+ Double_Matrix * coeffs,
+ Double_Matrix * eexps)
+{
+ all_elems.push_back(t);
+ if(level++ >= maxlevel)
+ return;
+
+ /*
+ p3
+
+ p13 p23
+
+ p1 p12 p2
+ */
+
+ adapt_point *p1 = t->p[0];
+ adapt_point *p2 = t->p[1];
+ adapt_point *p3 = t->p[2];
+ adapt_point *p12 = adapt_point::New((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5, 0,
+ coeffs, eexps);
+ adapt_point *p13 = adapt_point::New((p1->x + p3->x) * 0.5, (p1->y + p3->y) * 0.5, 0,
+ coeffs, eexps);
+ adapt_point *p23 = adapt_point::New((p3->x + p2->x) * 0.5, (p3->y + p2->y) * 0.5, 0,
+ coeffs, eexps);
+ adapt_triangle *t1 = new adapt_triangle(p1, p12, p13);
+ Recur_Create(t1, maxlevel, level, coeffs, eexps);
+ adapt_triangle *t2 = new adapt_triangle(p2, p23, p12);
+ Recur_Create(t2, maxlevel, level, coeffs, eexps);
+ adapt_triangle *t3 = new adapt_triangle(p3, p13, p23);
+ Recur_Create(t3, maxlevel, level, coeffs, eexps);
+ adapt_triangle *t4 = new adapt_triangle(p12, p23, p13);
+ Recur_Create(t4, maxlevel, level, coeffs, eexps);
+ t->e[0] = t1;
+ t->e[1] = t2;
+ t->e[2] = t3;
+ t->e[3] = t4;
+}
+
+void adapt_quad::Recur_Create(adapt_quad * q, int maxlevel, int level,
+ Double_Matrix * coeffs, Double_Matrix * eexps)
+{
+ all_elems.push_back(q);
+ if(level++ >= maxlevel)
+ return;
+
+ /*
+ p4 p34 p3
+
+ p14 pc p23
+
+ p1 p12 p2
+ */
+
+ adapt_point *p1 = q->p[0];
+ adapt_point *p2 = q->p[1];
+ adapt_point *p3 = q->p[2];
+ adapt_point *p4 = q->p[3];
+ adapt_point *p12 = adapt_point::New((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5, 0,
+ coeffs, eexps);
+ adapt_point *p23 = adapt_point::New((p2->x + p3->x) * 0.5, (p2->y + p3->y) * 0.5, 0,
+ coeffs, eexps);
+ adapt_point *p34 = adapt_point::New((p3->x + p4->x) * 0.5, (p3->y + p4->y) * 0.5, 0,
+ coeffs, eexps);
+ adapt_point *p14 = adapt_point::New((p1->x + p4->x) * 0.5, (p1->y + p4->y) * 0.5, 0,
+ coeffs, eexps);
+ adapt_point *pc = adapt_point::New((p1->x + p2->x + p3->x + p4->x) * 0.25,
+ (p1->y + p2->y + p3->y + p4->y) * 0.25, 0,
+ coeffs, eexps);
+ adapt_quad *q1 = new adapt_quad(p1, p12, pc, p14);
+ Recur_Create(q1, maxlevel, level, coeffs, eexps);
+ adapt_quad *q2 = new adapt_quad(p2, p23, pc, p12);
+ Recur_Create(q2, maxlevel, level, coeffs, eexps);
+ adapt_quad *q3 = new adapt_quad(p3, p34, pc, p23);
+ Recur_Create(q3, maxlevel, level, coeffs, eexps);
+ adapt_quad *q4 = new adapt_quad(p4, p14, pc, p34);
+ Recur_Create(q4, maxlevel, level, coeffs, eexps);
+ q->e[0] = q1;
+ q->e[1] = q2;
+ q->e[2] = q3;
+ q->e[3] = q4;
+}
+
+void adapt_tet::Recur_Create(adapt_tet * t, int maxlevel, int level,
+ Double_Matrix * coeffs, Double_Matrix * eexps)
+{
+ all_elems.push_back(t);
+ if(level++ >= maxlevel)
+ return;
+
+ adapt_point *p0 = t->p[0];
+ adapt_point *p1 = t->p[1];
+ adapt_point *p2 = t->p[2];
+ adapt_point *p3 = t->p[3];
+ adapt_point *pe0 = adapt_point::New((p0->x + p1->x) * 0.5, (p0->y + p1->y) * 0.5,
+ (p0->z + p1->z) * 0.5, coeffs, eexps);
+ adapt_point *pe1 = adapt_point::New((p0->x + p2->x) * 0.5, (p0->y + p2->y) * 0.5,
+ (p0->z + p2->z) * 0.5, coeffs, eexps);
+ adapt_point *pe2 = adapt_point::New((p0->x + p3->x) * 0.5, (p0->y + p3->y) * 0.5,
+ (p0->z + p3->z) * 0.5, coeffs, eexps);
+ adapt_point *pe3 = adapt_point::New((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5,
+ (p1->z + p2->z) * 0.5, coeffs, eexps);
+ adapt_point *pe4 = adapt_point::New((p1->x + p3->x) * 0.5, (p1->y + p3->y) * 0.5,
+ (p1->z + p3->z) * 0.5, coeffs, eexps);
+ adapt_point *pe5 = adapt_point::New((p2->x + p3->x) * 0.5, (p2->y + p3->y) * 0.5,
+ (p2->z + p3->z) * 0.5, coeffs, eexps);
+
+ adapt_tet *t1 = new adapt_tet(p0, pe0, pe2, pe1);
+ Recur_Create(t1, maxlevel, level, coeffs, eexps);
+ adapt_tet *t2 = new adapt_tet(p1, pe0, pe3, pe4);
+ Recur_Create(t2, maxlevel, level, coeffs, eexps);
+ adapt_tet *t3 = new adapt_tet(p2, pe3, pe1, pe5);
+ Recur_Create(t3, maxlevel, level, coeffs, eexps);
+ adapt_tet *t4 = new adapt_tet(p3, pe2, pe4, pe5);
+ Recur_Create(t4, maxlevel, level, coeffs, eexps);
+
+ adapt_tet *t5 = new adapt_tet(pe3, pe5, pe2, pe4);
+ Recur_Create(t5, maxlevel, level, coeffs, eexps);
+ adapt_tet *t6 = new adapt_tet(pe3, pe2, pe0, pe4);
+ Recur_Create(t6, maxlevel, level, coeffs, eexps);
+ adapt_tet *t7 = new adapt_tet(pe2, pe5, pe3, pe1);
+ Recur_Create(t7, maxlevel, level, coeffs, eexps);
+ adapt_tet *t8 = new adapt_tet(pe0, pe2, pe3, pe1);
+ Recur_Create(t8, maxlevel, level, coeffs, eexps);
+
+ t->e[0] = t1;
+ t->e[1] = t2;
+ t->e[2] = t3;
+ t->e[3] = t4;
+ t->e[4] = t5;
+ t->e[5] = t6;
+ t->e[6] = t7;
+ t->e[7] = t8;
+}
+
+void adapt_hex::Recur_Create(adapt_hex * h, int maxlevel, int level,
+ Double_Matrix * coeffs, Double_Matrix * eexps)
+{
+ all_elems.push_back(h);
+ if(level++ >= maxlevel)
+ return;
+
+ adapt_point *p0 = h->p[0];
+ adapt_point *p1 = h->p[1];
+ adapt_point *p2 = h->p[2];
+ adapt_point *p3 = h->p[3];
+ adapt_point *p4 = h->p[4];
+ adapt_point *p5 = h->p[5];
+ adapt_point *p6 = h->p[6];
+ adapt_point *p7 = h->p[7];
+ adapt_point *p01 = adapt_point::New((p0->x + p1->x) * 0.5, (p0->y + p1->y) * 0.5,
+ (p0->z + p1->z) * 0.5, coeffs, eexps);
+ adapt_point *p12 = adapt_point::New((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5,
+ (p1->z + p2->z) * 0.5, coeffs, eexps);
+ adapt_point *p23 = adapt_point::New((p2->x + p3->x) * 0.5, (p2->y + p3->y) * 0.5,
+ (p2->z + p3->z) * 0.5, coeffs, eexps);
+ adapt_point *p03 = adapt_point::New((p3->x + p0->x) * 0.5, (p3->y + p0->y) * 0.5,
+ (p3->z + p0->z) * 0.5, coeffs, eexps);
+ adapt_point *p45 = adapt_point::New((p4->x + p5->x) * 0.5, (p4->y + p5->y) * 0.5,
+ (p4->z + p5->z) * 0.5, coeffs, eexps);
+ adapt_point *p56 = adapt_point::New((p5->x + p6->x) * 0.5, (p5->y + p6->y) * 0.5,
+ (p5->z + p6->z) * 0.5, coeffs, eexps);
+ adapt_point *p67 = adapt_point::New((p6->x + p7->x) * 0.5, (p6->y + p7->y) * 0.5,
+ (p6->z + p7->z) * 0.5, coeffs, eexps);
+ adapt_point *p47 = adapt_point::New((p7->x + p4->x) * 0.5, (p7->y + p4->y) * 0.5,
+ (p7->z + p4->z) * 0.5, coeffs, eexps);
+ adapt_point *p04 = adapt_point::New((p4->x + p0->x) * 0.5, (p4->y + p0->y) * 0.5,
+ (p4->z + p0->z) * 0.5, coeffs, eexps);
+ adapt_point *p15 = adapt_point::New((p5->x + p1->x) * 0.5, (p5->y + p1->y) * 0.5,
+ (p5->z + p1->z) * 0.5, coeffs, eexps);
+ adapt_point *p26 = adapt_point::New((p6->x + p2->x) * 0.5, (p6->y + p2->y) * 0.5,
+ (p6->z + p2->z) * 0.5, coeffs, eexps);
+ adapt_point *p37 = adapt_point::New((p7->x + p3->x) * 0.5, (p7->y + p3->y) * 0.5,
+ (p7->z + p3->z) * 0.5, coeffs, eexps);
+
+ adapt_point *p0145 = adapt_point::New((p45->x + p01->x) * 0.5, (p45->y + p01->y) * 0.5,
+ (p45->z + p01->z) * 0.5, coeffs, eexps);
+ adapt_point *p1256 = adapt_point::New((p12->x + p56->x) * 0.5, (p12->y + p56->y) * 0.5,
+ (p12->z + p56->z) * 0.5, coeffs, eexps);
+ adapt_point *p2367 = adapt_point::New((p23->x + p67->x) * 0.5, (p23->y + p67->y) * 0.5,
+ (p23->z + p67->z) * 0.5, coeffs, eexps);
+ adapt_point *p0347 = adapt_point::New((p03->x + p47->x) * 0.5, (p03->y + p47->y) * 0.5,
+ (p03->z + p47->z) * 0.5, coeffs, eexps);
+ adapt_point *p4756 = adapt_point::New((p47->x + p56->x) * 0.5, (p47->y + p56->y) * 0.5,
+ (p47->z + p56->z) * 0.5, coeffs, eexps);
+ adapt_point *p0312 = adapt_point::New((p03->x + p12->x) * 0.5, (p03->y + p12->y) * 0.5,
+ (p03->z + p12->z) * 0.5, coeffs, eexps);
+
+ adapt_point *pc =
+ adapt_point::New((p0->x + p1->x + p2->x + p3->x + p4->x + p5->x + p6->x + p7->x) * 0.125,
+ (p0->y + p1->y + p2->y + p3->y + p4->y + p5->y + p6->y + p7->y) * 0.125,
+ (p0->z + p1->z + p2->z + p3->z + p4->z + p5->z + p6->z + p7->z) * 0.125,
+ coeffs, eexps);
+
+ adapt_hex *h1 = new adapt_hex(p0, p01, p0312, p03, p04, p0145, pc, p0347); //p0
+ Recur_Create(h1, maxlevel, level, coeffs, eexps);
+ adapt_hex *h2 = new adapt_hex(p01, p0145, p15, p1, p0312, pc, p1256, p12); //p1
+ Recur_Create(h2, maxlevel, level, coeffs, eexps);
+ adapt_hex *h3 = new adapt_hex(p04, p4, p45, p0145, p0347, p47, p4756, pc); //p4
+ Recur_Create(h3, maxlevel, level, coeffs, eexps);
+ adapt_hex *h4 = new adapt_hex(p0145, p45, p5, p15, pc, p4756, p56, p1256); //p5
+ Recur_Create(h4, maxlevel, level, coeffs, eexps);
+ adapt_hex *h5 = new adapt_hex(p0347, p47, p4756, pc, p37, p7, p67, p2367); //p7
+ Recur_Create(h5, maxlevel, level, coeffs, eexps);
+ adapt_hex *h6 = new adapt_hex(pc, p4756, p56, p1256, p2367, p67, p6, p26); //p6
+ Recur_Create(h6, maxlevel, level, coeffs, eexps);
+ adapt_hex *h7 = new adapt_hex(p03, p0347, pc, p0312, p3, p37, p2367, p23); //p3
+ Recur_Create(h7, maxlevel, level, coeffs, eexps);
+ adapt_hex *h8 = new adapt_hex(p0312, pc, p1256, p12, p23, p2367, p26, p2); //p2
+ Recur_Create(h8, maxlevel, level, coeffs, eexps);
+
+ h->e[0] = h1;
+ h->e[1] = h2;
+ h->e[2] = h3;
+ h->e[3] = h4;
+ h->e[4] = h5;
+ h->e[5] = h6;
+ h->e[6] = h7;
+ h->e[7] = h8;
+}
+
+void adapt_edge::Error(double AVG, double tol)
+{
+ adapt_edge *e = *all_elems.begin();
+ Recur_Error(e, AVG, tol);
+}
+
+void adapt_triangle::Error(double AVG, double tol)
+{
+ adapt_triangle *t = *all_elems.begin();
+ Recur_Error(t, AVG, tol);
+}
+
+void adapt_quad::Error(double AVG, double tol)
+{
+ adapt_quad *q = *all_elems.begin();
+ Recur_Error(q, AVG, tol);
+}
+
+void adapt_tet::Error(double AVG, double tol)
+{
+ adapt_tet *t = *all_elems.begin();
+ Recur_Error(t, AVG, tol);
+}
+void adapt_hex::Error(double AVG, double tol)
+{
+ adapt_hex *h = *all_elems.begin();
+ Recur_Error(h, AVG, tol);
+}
+
+void adapt_edge::Recur_Error(adapt_edge * e, double AVG, double tol)
+{
+ if(!e->e[0])
+ e->visible = true;
+ else {
+ double vr;
+ if(!e->e[0]->e[0]) {
+ double v1 = e->e[0]->V();
+ double v2 = e->e[1]->V();
+ vr = (v1 + v2) / 2.;
+ double v = e->V();
+ if(fabs(v - vr) > AVG * tol)
+ //if (fabs(v - vr) > ((fabs(v) + fabs(vr) + AVG * tol) * tol))
+ {
+ e->visible = false;
+ Recur_Error(e->e[0], AVG, tol);
+ Recur_Error(e->e[1], AVG, tol);
+ }
+ else
+ e->visible = true;
+ }
+ else {
+ double v11 = e->e[0]->e[0]->V();
+ double v12 = e->e[0]->e[1]->V();
+
+ double v21 = e->e[1]->e[0]->V();
+ double v22 = e->e[1]->e[1]->V();
+
+ double vr1 = (v11 + v12) / 2.;
+ double vr2 = (v21 + v22) / 2.;
+ vr = (vr1 + vr2) / 2.;
+ if(fabs(e->e[0]->V() - vr1) > AVG * tol ||
+ fabs(e->e[1]->V() - vr2) > AVG * tol ||
+ fabs(e->V() - vr) > AVG * tol) {
+ e->visible = false;
+ Recur_Error(e->e[0], AVG, tol);
+ Recur_Error(e->e[1], AVG, tol);
+ }
+ else
+ e->visible = true;
+ }
+ }
+}
+
+void adapt_triangle::Recur_Error(adapt_triangle * t, double AVG, double tol)
+{
+ if(!t->e[0])
+ t->visible = true;
+ else {
+ double vr;
+ if(!t->e[0]->e[0]) {
+ double v1 = t->e[0]->V();
+ double v2 = t->e[1]->V();
+ double v3 = t->e[2]->V();
+ double v4 = t->e[3]->V();
+ vr = (2 * v1 + 2 * v2 + 2 * v3 + v4) / 7.;
+ double v = t->V();
+ if(fabs(v - vr) > AVG * tol)
+ //if (fabs(v - vr) > ((fabs(v) + fabs(vr) + AVG * tol) * tol))
+ {
+ t->visible = false;
+ Recur_Error(t->e[0], AVG, tol);
+ Recur_Error(t->e[1], AVG, tol);
+ Recur_Error(t->e[2], AVG, tol);
+ Recur_Error(t->e[3], AVG, tol);
+ }
+ else
+ t->visible = true;
+ }
+ else {
+ double v11 = t->e[0]->e[0]->V();
+ double v12 = t->e[0]->e[1]->V();
+ double v13 = t->e[0]->e[2]->V();
+ double v14 = t->e[0]->e[3]->V();
+ double v21 = t->e[1]->e[0]->V();
+ double v22 = t->e[1]->e[1]->V();
+ double v23 = t->e[1]->e[2]->V();
+ double v24 = t->e[1]->e[3]->V();
+ double v31 = t->e[2]->e[0]->V();
+ double v32 = t->e[2]->e[1]->V();
+ double v33 = t->e[2]->e[2]->V();
+ double v34 = t->e[2]->e[3]->V();
+ double v41 = t->e[3]->e[0]->V();
+ double v42 = t->e[3]->e[1]->V();
+ double v43 = t->e[3]->e[2]->V();
+ double v44 = t->e[3]->e[3]->V();
+ double vr1 = (2 * v11 + 2 * v12 + 2 * v13 + v14) / 7.;
+ double vr2 = (2 * v21 + 2 * v22 + 2 * v23 + v24) / 7.;
+ double vr3 = (2 * v31 + 2 * v32 + 2 * v33 + v34) / 7.;
+ double vr4 = (2 * v41 + 2 * v42 + 2 * v43 + v44) / 7.;
+ vr = (2 * vr1 + 2 * vr2 + 2 * vr3 + vr4) / 7.;
+ if(fabs(t->e[0]->V() - vr1) > AVG * tol ||
+ fabs(t->e[1]->V() - vr2) > AVG * tol ||
+ fabs(t->e[2]->V() - vr3) > AVG * tol ||
+ fabs(t->e[3]->V() - vr4) > AVG * tol ||
+ fabs(t->V() - vr) > AVG * tol)
+ //if (fabs(t->e[0]->V() - vr1) > (fabs(t->e[0]->V())+fabs(vr1)+AVG * tol)*tol ||
+ // fabs(t->e[1]->V() - vr2) > (fabs(t->e[1]->V())+fabs(vr2)+AVG * tol)*tol ||
+ // fabs(t->e[2]->V() - vr3) > (fabs(t->e[2]->V())+fabs(vr3)+AVG * tol)*tol ||
+ // fabs(t->e[3]->V() - vr4) > (fabs(t->e[3]->V())+fabs(vr4)+AVG * tol)*tol ||
+ // fabs(t->V() - vr) > (fabs(t->V())+fabs(vr)+AVG * tol) *tol)
+ {
+ t->visible = false;
+ Recur_Error(t->e[0], AVG, tol);
+ Recur_Error(t->e[1], AVG, tol);
+ Recur_Error(t->e[2], AVG, tol);
+ Recur_Error(t->e[3], AVG, tol);
+ }
+ else
+ t->visible = true;
+ }
+ }
+}
+
+void adapt_quad::Recur_Error(adapt_quad * q, double AVG, double tol)
+{
+ if(!q->e[0])
+ q->visible = true;
+ else {
+ double vr;
+ if(!q->e[0]->e[0]) {
+ double v1 = q->e[0]->V();
+ double v2 = q->e[1]->V();
+ double v3 = q->e[2]->V();
+ double v4 = q->e[3]->V();
+ vr = (v1 + v2 + v3 + v4) / 4.;
+ double v = q->V();
+ if(fabs(v - vr) > AVG * tol)
+ //if (fabs(v - vr) > ((fabs(v) + fabs(vr) + AVG * tol) * tol))
+ {
+ q->visible = false;
+ Recur_Error(q->e[0], AVG, tol);
+ Recur_Error(q->e[1], AVG, tol);
+ Recur_Error(q->e[2], AVG, tol);
+ Recur_Error(q->e[3], AVG, tol);
+ }
+ else
+ q->visible = true;
+ }
+ else {
+ double v11 = q->e[0]->e[0]->V();
+ double v12 = q->e[0]->e[1]->V();
+ double v13 = q->e[0]->e[2]->V();
+ double v14 = q->e[0]->e[3]->V();
+ double v21 = q->e[1]->e[0]->V();
+ double v22 = q->e[1]->e[1]->V();
+ double v23 = q->e[1]->e[2]->V();
+ double v24 = q->e[1]->e[3]->V();
+ double v31 = q->e[2]->e[0]->V();
+ double v32 = q->e[2]->e[1]->V();
+ double v33 = q->e[2]->e[2]->V();
+ double v34 = q->e[2]->e[3]->V();
+ double v41 = q->e[3]->e[0]->V();
+ double v42 = q->e[3]->e[1]->V();
+ double v43 = q->e[3]->e[2]->V();
+ double v44 = q->e[3]->e[3]->V();
+ double vr1 = (v11 + v12 + v13 + v14) / 4.;
+ double vr2 = (v21 + v22 + v23 + v24) / 4.;
+ double vr3 = (v31 + v32 + v33 + v34) / 4.;
+ double vr4 = (v41 + v42 + v43 + v44) / 4.;
+ vr = (vr1 + vr2 + vr3 + vr4) / 4.;
+ if(fabs(q->e[0]->V() - vr1) > AVG * tol ||
+ fabs(q->e[1]->V() - vr2) > AVG * tol ||
+ fabs(q->e[2]->V() - vr3) > AVG * tol ||
+ fabs(q->e[3]->V() - vr4) > AVG * tol ||
+ fabs(q->V() - vr) > AVG * tol)
+ //if (fabs(t->e[0]->V() - vr1) > (fabs(t->e[0]->V())+fabs(vr1)+AVG * tol)*tol ||
+ // fabs(t->e[1]->V() - vr2) > (fabs(t->e[1]->V())+fabs(vr2)+AVG * tol)*tol ||
+ // fabs(t->e[2]->V() - vr3) > (fabs(t->e[2]->V())+fabs(vr3)+AVG * tol)*tol ||
+ // fabs(t->e[3]->V() - vr4) > (fabs(t->e[3]->V())+fabs(vr4)+AVG * tol)*tol ||
+ // fabs(t->V() - vr) > (fabs(t->V())+fabs(vr)+AVG * tol) *tol)
+ {
+ q->visible = false;
+ Recur_Error(q->e[0], AVG, tol);
+ Recur_Error(q->e[1], AVG, tol);
+ Recur_Error(q->e[2], AVG, tol);
+ Recur_Error(q->e[3], AVG, tol);
+ }
+ else
+ q->visible = true;
+ }
+ }
+}
+
+void adapt_tet::Recur_Error(adapt_tet * t, double AVG, double tol)
+{
+ if(!t->e[0])
+ t->visible = true;
+ else {
+ const double v1 = t->e[0]->V();
+ const double v2 = t->e[1]->V();
+ const double v3 = t->e[2]->V();
+ const double v4 = t->e[3]->V();
+ const double v5 = t->e[4]->V();
+ const double v6 = t->e[5]->V();
+ const double v7 = t->e[6]->V();
+ const double v8 = t->e[7]->V();
+ const double vr = (v1 + v2 + v3 + v4 + v5 + v6 + v7 + v8) * .125;
+ const double v = t->V();
+ if(!t->e[0]->e[0]) {
+ if(fabs(v - vr) > AVG * tol) {
+ t->visible = false;
+ Recur_Error(t->e[0], AVG, tol);
+ Recur_Error(t->e[1], AVG, tol);
+ Recur_Error(t->e[2], AVG, tol);
+ Recur_Error(t->e[3], AVG, tol);
+ Recur_Error(t->e[4], AVG, tol);
+ Recur_Error(t->e[5], AVG, tol);
+ Recur_Error(t->e[6], AVG, tol);
+ Recur_Error(t->e[7], AVG, tol);
+ }
+ else
+ t->visible = true;
+ }
+ else {
+ double vi[8][8];
+ for(int k = 0; k < 8; k++)
+ for(int l = 0; l < 8; l++)
+ vi[k][l] = t->e[k]->e[l]->V();
+ double vri[8];
+ for(int k = 0; k < 8; k++) {
+ vri[k] = 0.0;
+ for(int l = 0; l < 8; l++) {
+ vri[k] += vi[k][l];
+ }
+ vri[k] /= 8.0;
+ }
+ if(fabs(t->e[0]->V() - vri[0]) > AVG * tol ||
+ fabs(t->e[1]->V() - vri[1]) > AVG * tol ||
+ fabs(t->e[2]->V() - vri[2]) > AVG * tol ||
+ fabs(t->e[3]->V() - vri[3]) > AVG * tol ||
+ fabs(t->e[4]->V() - vri[4]) > AVG * tol ||
+ fabs(t->e[5]->V() - vri[5]) > AVG * tol ||
+ fabs(t->e[6]->V() - vri[6]) > AVG * tol ||
+ fabs(t->e[7]->V() - vri[7]) > AVG * tol ||
+ fabs(v - vr) > AVG * tol) {
+ t->visible = false;
+ Recur_Error(t->e[0], AVG, tol);
+ Recur_Error(t->e[1], AVG, tol);
+ Recur_Error(t->e[2], AVG, tol);
+ Recur_Error(t->e[3], AVG, tol);
+ Recur_Error(t->e[4], AVG, tol);
+ Recur_Error(t->e[5], AVG, tol);
+ Recur_Error(t->e[6], AVG, tol);
+ Recur_Error(t->e[7], AVG, tol);
+ }
+ else
+ t->visible = true;
+ }
+ }
+}
+
+void adapt_hex::Recur_Error(adapt_hex * h, double AVG, double tol)
+{
+ if(!h->e[0])
+ h->visible = true;
+ else {
+ double vr;
+ double v1 = h->e[0]->V();
+ double v2 = h->e[1]->V();
+ double v3 = h->e[2]->V();
+ double v4 = h->e[3]->V();
+ double v5 = h->e[4]->V();
+ double v6 = h->e[5]->V();
+ double v7 = h->e[6]->V();
+ double v8 = h->e[7]->V();
+ vr = (v1 + v2 + v3 + v4 + v5 + v6 + v7 + v8) * .125;
+ double v = h->V();
+ if(fabs(v - vr) > AVG * tol) {
+ h->visible = false;
+ Recur_Error(h->e[0], AVG, tol);
+ Recur_Error(h->e[1], AVG, tol);
+ Recur_Error(h->e[2], AVG, tol);
+ Recur_Error(h->e[3], AVG, tol);
+ Recur_Error(h->e[4], AVG, tol);
+ Recur_Error(h->e[5], AVG, tol);
+ Recur_Error(h->e[6], AVG, tol);
+ Recur_Error(h->e[7], AVG, tol);
+ }
+ else
+ h->visible = true;
+ }
+}
+
+static double t0, t1, t2, t3;
+
+template < class ELEM >
+int Adaptive_Post_View::zoomElement(Post_View * view,
+ int ielem,
+ int level,
+ int levelmax,
+ GMSH_Post_Plugin * plug,
+ List_T * theList, int *counter)
+{
+ const int nbNod = ELEM::nbNod;
+
+ typename std::set < adapt_point >::iterator it = adapt_point::all_points.begin();
+ typename std::set < adapt_point >::iterator ite = adapt_point::all_points.end();
+
+ double c0 = Cpu();
+
+ const int N = _coefs->size1();
+ Double_Vector val(N), res(adapt_point::all_points.size());
+ Double_Vector valx(N), resx(adapt_point::all_points.size());
+ Double_Vector valy(N), resy(adapt_point::all_points.size());
+ Double_Vector valz(N), resz(adapt_point::all_points.size());
+ Double_Matrix xyz(nbNod,3);
+ Double_Matrix XYZ(adapt_point::all_points.size(),3);
+
+ for(int k = 0; k < nbNod; ++k){
+ xyz(k, 0) = (*_STposX)(ielem, k);
+ xyz(k, 1) = (*_STposY)(ielem, k);
+ xyz(k, 2) = (*_STposZ)(ielem, k);
+ }
+
+ for (int k = 0; k < N; ++k){
+ val(k) = (*_STval)(ielem, k);
+ }
+
+ _Interpolate->mult(val,res);
+
+ if(_STvalX){
+ for(int k = 0; k < N; ++k){
+ valx(k) = (*_STvalX)(ielem, k);
+ valy(k) = (*_STvalY)(ielem, k);
+ valz(k) = (*_STvalZ)(ielem, k);
+ }
+ _Interpolate->mult(valx, resx);
+ _Interpolate->mult(valy, resy);
+ _Interpolate->mult(valz, resz);
+ }
+
+ _Geometry->mult(xyz, XYZ);
+
+ double c1 = Cpu();
+
+ int kk = 0;
+ for (; it !=ite; ++it){
+ adapt_point *p = (adapt_point*) &(*it);
+ p->val = res(kk);
+ if (_STvalX){
+ p->valx = resx(kk);
+ p->valy = resy(kk);
+ p->valz = resz(kk);
+ }
+ p->val = res(kk);
+ p->X = XYZ(kk, 0);
+ p->Y = XYZ(kk, 1);
+ p->Z = XYZ(kk, 2);
+ if (minval > p->val) minval = p->val;
+ if (maxval < p->val) maxval = p->val;
+ kk++;
+ }
+ double c2 = Cpu();
+
+ typename std::list < ELEM * >::iterator itt = ELEM::all_elems.begin();
+ typename std::list < ELEM * >::iterator itte = ELEM::all_elems.end();
+
+ for(; itt != itte; itt++) {
+ (*itt)->visible = false;
+ }
+
+ if(!plug || tol != 0.0) {
+ ELEM::Error(maxval - minval, tol);
+ }
+
+ if(plug)
+ plug->assign_specific_visibility();
+
+ double c3 = Cpu();
+ itt = ELEM::all_elems.begin();
+
+ for(; itt != itte; itt++) {
+ if((*itt)->visible && !(*itt)->e[0] && level != levelmax)
+ return 0;
+ }
+
+ itt = ELEM::all_elems.begin();
+ adapt_point **p;
+
+ for(; itt != itte; itt++){
+ if((*itt)->visible){
+ p = (*itt)->p;
+ for(int k=0; k < nbNod; ++k) List_Add(theList, &p[k]->X);
+ for(int k=0; k < nbNod; ++k) List_Add(theList, &p[k]->Y);
+ for(int k=0; k < nbNod; ++k) List_Add(theList, &p[k]->Z);
+ if(_STvalX){
+ for(int k = 0; k < nbNod; ++k){
+ List_Add(theList, &p[k]->valx);
+ List_Add(theList, &p[k]->valy);
+ List_Add(theList, &p[k]->valz);
+ }
+ }
+ else{
+ for (int k = 0; k < nbNod; ++k) List_Add(theList, &p[k]->val);
+ }
+ (*counter)++;
+ }
+ }
+
+ double c4 = Cpu();
+
+ t0 += c1 - c0;
+ t1 += c2 - c1;
+ t2 += c3 - c2;
+ t3 += c4 - c3;
+
+ return 1;
+}
+
+
+/*
+ We first do the adaptive stuff at level 2 and will only process
+ elements that have reached the maximal recursion level
+
+ We compute first the matrix at max recursion level (those should be stored
+ on disk once in the GMSHPLUGINSHOME directory, i'll do that at some point).
+*/
+
+void Adaptive_Post_View::setAdaptiveResolutionLevel(Post_View * view,
+ int level,
+ GMSH_Post_Plugin * plug)
+{
+
+ if(presentTol == tol && presentZoomLevel == level && !plug)
+ return;
+
+ printf
+ ("calling setAdaptive with level = %d and plug = %p tol %g presentTol %g\n",
+ level, (void *)plug, tol, presentTol);
+
+ int *done = new int[_STposX->size1()];
+ for(int i = 0; i < _STposX->size1(); ++i)
+ done[i] = 0;
+ int level_act = (level > 2) ? 2 : level;
+ int nbelm = _STposX->size1();
+
+ int TYP = 0;
+ if(view->NbSL) {
+ TYP = 7;
+ List_Delete(view->SL);
+ view->NbSL = 0;
+ view->SL = List_Create(nbelm * 8, nbelm, sizeof(double));
+ }
+ if (view->NbVT) {
+ TYP = 5;
+ List_Delete(view->VT);
+ view->NbVT = 0;
+ view->VT =List_Create ( nbelm * 36, nbelm , sizeof(double));
+ }
+ if(view->NbST) {
+ TYP = 1;
+ List_Delete(view->ST);
+ view->NbST = 0;
+ view->ST = List_Create(nbelm * 4, nbelm, sizeof(double));
+ }
+ if(view->NbSS) {
+ TYP = 3;
+ List_Delete(view->SS);
+ view->NbSS = 0;
+ view->SS = List_Create(nbelm * 4, nbelm, sizeof(double));
+ }
+ if(view->NbSQ) {
+ TYP = 2;
+ List_Delete(view->SQ);
+ view->NbSQ = 0;
+ view->SQ = List_Create(nbelm * 20, nbelm * 20, sizeof(double));
+ }
+ if(view->NbVQ) {
+ TYP = 6;
+ List_Delete(view->VQ);
+ view->NbVQ = 0;
+ view->VQ = List_Create(nbelm * 60, nbelm * 20, sizeof(double));
+ }
+ if(view->NbSH) {
+ TYP = 4;
+ List_Delete(view->SH);
+ view->NbSH = 0;
+ view->SH = List_Create(nbelm * 4, nbelm, sizeof(double));
+ }
+
+ view->NbTimeStep = 1;
+
+ t0 = t1 = t2 = t3 = 0;
+
+ while (1){
+ if(TYP == 7)
+ setAdaptiveResolutionLevel_TEMPL < adapt_edge > (view, level_act, level, plug,
+ &(view->SL), &(view->NbSL), done);
+ if (TYP == 1)
+ setAdaptiveResolutionLevel_TEMPL <adapt_triangle> (view, level_act, level, plug,
+ &(view->ST), &(view->NbST), done);
+ if (TYP == 5)
+ setAdaptiveResolutionLevel_TEMPL <adapt_triangle> (view, level_act, level, plug,
+ &(view->VT), &(view->NbVT), done);
+ if (TYP == 2)
+ setAdaptiveResolutionLevel_TEMPL <adapt_quad> (view, level_act, level, plug,
+ &(view->SQ), &(view->NbSQ), done);
+ if (TYP == 6)
+ setAdaptiveResolutionLevel_TEMPL <adapt_quad> (view, level_act, level, plug,
+ &(view->VQ), &(view->NbVQ), done);
+ if (TYP == 4)
+ setAdaptiveResolutionLevel_TEMPL <adapt_hex> (view, level_act, level, plug,
+ &(view->SH), &(view->NbSH), done);
+ if (TYP == 3)
+ setAdaptiveResolutionLevel_TEMPL <adapt_tet> (view, level_act, level, plug,
+ &(view->SS), &(view->NbSS), done);
+ int nbDone = 0;
+ for(int i=0; i < _STposX->size1(); ++i) nbDone += done[i];
+ printf("adaptive %d %d %d %d\n", level, level_act, nbDone, _STposX->size1());
+ if (nbDone ==_STposX->size1()) break;
+ if (level_act >= level) break;
+ level_act ++;
+ }
+
+ // recompute min/max, etc.:
+ view->Min = VAL_INF;
+ view->Max = -VAL_INF;
+ EndView(view, 0, view->FileName, view->Name);
+ view->Changed = 1;
+ presentZoomLevel = level;
+ presentTol = tol;
+
+ printf("finished %g %g %g %g\n", t0, t1, t2, t3);
+ delete[]done;
+}
+
+template < class ELEM >
+void Adaptive_Post_View::setAdaptiveResolutionLevel_TEMPL(Post_View * view,
+ int level,
+ int levelmax,
+ GMSH_Post_Plugin *
+ plug,
+ List_T ** myList,
+ int *counter,
+ int *done)
+{
+ const int N = _coefs->size1();
+
+ const int nbelm = _STposX->size1();
+
+ double sf[100];
+ ELEM::Create(level, _coefs, _eexps);
+ std::set < adapt_point >::iterator it = adapt_point::all_points.begin();
+ std::set < adapt_point >::iterator ite = adapt_point::all_points.end();
+
+ if(_Interpolate)
+ delete _Interpolate;
+ if(_Geometry)
+ delete _Geometry;
+ _Interpolate = new Double_Matrix(adapt_point::all_points.size(), N);
+ _Geometry = new Double_Matrix(adapt_point::all_points.size(), ELEM::nbNod);
+
+ int kk = 0;
+ for(; it != ite; ++it) {
+ adapt_point *p = (adapt_point *) & (*it);
+ for(int k = 0; k < N; ++k) {
+ (*_Interpolate) (kk, k) = p->shape_functions[k];
+ }
+ ELEM::GSF(p->x, p->y, p->z, sf);
+ for(int k = 0; k < ELEM::nbNod; k++)
+ (*_Geometry) (kk, k) = sf[k];
+ kk++;
+ }
+
+ for(int i = 0; i < nbelm; ++i) {
+ if(!done[i])
+ done[i] =
+ zoomElement < ELEM > (view, i, level, levelmax, plug, *myList,
+ counter);
+ }
+}
+
+void computeShapeFunctions(Double_Matrix * coeffs, Double_Matrix * eexps,
+ double u, double v, double w, double *sf)
+{
+ static double powsuvw[256];
+ for(int j = 0; j < coeffs->size2(); ++j) {
+ double powu = (*eexps) (j, 0);
+ double powv = (*eexps) (j, 1);
+ double poww = (*eexps) (j, 2);
+ powsuvw[j] = pow(u, powu) * pow(v, powv) * pow(w, poww);
+ }
+
+ for(int i = 0; i < coeffs->size1(); ++i) {
+ sf[i] = 0.0;
+ for(int j = 0; j < coeffs->size2(); ++j) {
+ sf[i] += (*coeffs) (i, j) * powsuvw[j];
+ }
+ }
+}
+
+void Adaptive_Post_View::initWithLowResolution(Post_View * view)
+{
+ List_T *myList;
+ int nbelm;
+ int nbnod;
+
+ int nbComp = 1;
+
+ if (view->NbST){
+ myList = view->ST;
+ nbelm = view->NbST;
+ nbnod = 3;
+ }
+ else if(view->NbSL){
+ myList = view->SL;
+ nbelm = view->NbSL;
+ nbnod = 2;
+ }
+ else if (view->NbVT){
+ myList = view->VT;
+ nbelm = view->NbVT;
+ nbnod = 3;
+ nbComp = 3;
+ }
+ else if (view->NbVQ){
+ myList = view->VQ;
+ nbelm = view->NbVQ;
+ nbnod = 4;
+ nbComp = 3;
+ }
+ else if (view->NbSS){
+ myList = view->SS;
+ nbelm = view->NbSS;
+ nbnod = 4;
+ }
+ else if (view->NbSQ){
+ myList = view->SQ;
+ nbelm = view->NbSQ;
+ nbnod = 4;
+ }
+ else if (view->NbSH){
+ myList = view->SH;
+ nbelm = view->NbSH;
+ nbnod = 8;
+ }
+ else return;
+
+ // if there exists a polynomial representation
+ // of the geometry , then use it
+ if (_coefsGeom)
+ {
+ nbnod = _coefsGeom -> size1 ();
+ }
+
+
+ minval = VAL_INF;
+ maxval = -VAL_INF;
+
+ int nb = List_Nbr(myList) / (nbelm);
+
+ _STposX = new Double_Matrix ( nbelm , nbnod );
+ _STposY = new Double_Matrix ( nbelm , nbnod );
+ _STposZ = new Double_Matrix ( nbelm , nbnod );
+ _STval = new Double_Matrix ( nbelm , (nb-3*nbnod)/nbComp );
+
+ if (nbComp == 3){
+ _STvalX = new Double_Matrix ( nbelm , (nb-3*nbnod)/nbComp );
+ _STvalY = new Double_Matrix ( nbelm , (nb-3*nbnod)/nbComp );
+ _STvalZ = new Double_Matrix ( nbelm , (nb-3*nbnod)/nbComp );
+ }
+ /// Store non interpolated data
+
+ int k=0;
+ for (int i=0;i<List_Nbr(myList);i+=nb){
+ double *x = (double*)List_Pointer_Fast (myList,i);
+ double *y = (double*)List_Pointer_Fast (myList,i+nbnod);
+ double *z = (double*)List_Pointer_Fast (myList,i+2*nbnod);
+ for (int NN=0;NN<nbnod;NN++){
+ (*_STposX) ( k , NN) = x[NN];
+ (*_STposY) ( k , NN) = y[NN];
+ (*_STposZ) ( k , NN) = z[NN];
+ }
+ double *val = (double*)List_Pointer_Fast (myList,i+3*nbnod);
+ if (nbComp == 1){
+ for (int j=0;j<(nb-3*nbnod)/nbComp;j++){
+ (*_STval)(k,j)=val[j];
+ }
+ }
+ else if (nbComp == 3){
+ int size = (nb-3*nbnod)/3;
+ for (int j=0;j<size;j++){
+ int index1 = j;
+ int index2 = j+size;
+ int index3 = j+2*size;
+ // adaptation of the visualization mesh bases on the norm squared of the vector
+ (*_STval)(k,j)=(val[index1]*val[index1] + val[index2]*val[index2] +
+ val[index3]*val[index3]);
+ (*_STvalX)(k,j)=val[index1];
+ (*_STvalY)(k,j)=val[index2];
+ (*_STvalZ)(k,j)=val[index3];
+ }
+ }
+ k++;
+ }
+ setAdaptiveResolutionLevel(view, 0);
+}
+
+Adaptive_Post_View::Adaptive_Post_View(Post_View * view,
+ List_T * _c,
+ List_T * _pol,
+ List_T * _cGeom,
+ List_T * _polGeom)
+ :tol(1.e-3),_coefsGeom(0),_eexpsGeom(0)
+{
+
+ _Interpolate = _Geometry = 0;
+ _coefs = new Double_Matrix ( List_Nbr (_c) , List_Nbr (_c) );
+ _eexps = new Double_Matrix ( List_Nbr (_c) , 3 );
+
+ _STvalX = _STvalY = _STvalZ =0;
+
+ for(int i = 0; i < List_Nbr(_c); ++i) {
+ List_T **line = (List_T **) List_Pointer_Fast(_c, i);
+ List_T **eexp = (List_T **) List_Pointer_Fast(_pol, i);
+
+ double dpowu, dpowv, dpoww;
+
+ List_Read(*eexp, 0, &dpowu);
+ List_Read(*eexp, 1, &dpowv);
+ List_Read(*eexp, 2, &dpoww);
+
+ (*_eexps) (i, 0) = dpowu;
+ (*_eexps) (i, 1) = dpowv;
+ (*_eexps) (i, 2) = dpoww;
+
+ for(int j = 0; j < List_Nbr(*line); ++j) {
+ double val;
+ List_Read(*line, j, &val);
+ (*_coefs) (i, j) = val;
+ }
+ }
+
+ if (_cGeom && _polGeom)
+ {
+ _coefsGeom = new Double_Matrix ( List_Nbr (_cGeom) , List_Nbr (_cGeom) );
+ _eexpsGeom = new Double_Matrix ( List_Nbr (_cGeom) , 3 );
+ for(int i = 0; i < List_Nbr(_cGeom); ++i) {
+ List_T **line = (List_T **) List_Pointer_Fast(_cGeom, i);
+ List_T **eexp = (List_T **) List_Pointer_Fast(_polGeom, i);
+
+ double dpowu, dpowv, dpoww;
+
+ List_Read(*eexp, 0, &dpowu);
+ List_Read(*eexp, 1, &dpowv);
+ List_Read(*eexp, 2, &dpoww);
+
+ (*_eexpsGeom) (i, 0) = dpowu;
+ (*_eexpsGeom) (i, 1) = dpowv;
+ (*_eexpsGeom) (i, 2) = dpoww;
+
+ for(int j = 0; j < List_Nbr(*line); ++j) {
+ double val;
+ List_Read(*line, j, &val);
+ (*_coefsGeom) (i, j) = val;
+ }
+ }
+ }
+ initWithLowResolution(view);
+}
+
+Adaptive_Post_View::~Adaptive_Post_View()
+{
+ delete _coefs;
+ delete _eexps;
+ delete _STposX;
+ delete _STposY;
+ delete _STposZ;
+ delete _STval;
+
+ if (_coefsGeom)
+ delete _coefsGeom;
+ if (_eexpsGeom)
+ delete _eexpsGeom;
+
+ if(_Interpolate)
+ delete _Interpolate;
+ if(_Geometry)
+ delete _Geometry;
+ cleanElement < adapt_edge > ();
+ cleanElement < adapt_triangle > ();
+ cleanElement < adapt_tet > ();
+ cleanElement < adapt_hex > ();
+ cleanElement < adapt_quad > ();
+}
diff --git a/Common/AdaptiveViews.h b/Common/AdaptiveViews.h
new file mode 100644
index 0000000000000000000000000000000000000000..53effd4b2860d3f1f5b777ee002aa04de475c936
--- /dev/null
+++ b/Common/AdaptiveViews.h
@@ -0,0 +1,320 @@
+#ifndef _ADAPTIVE_VIEWS_H_
+#define _ADAPTIVE_VIEWS_H_
+
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include <list>
+#include "List.h"
+#include "GmshMatrix.h"
+
+#define MAX_LEVEL_OF_ZOOM 8
+
+class Post_View;
+class GMSH_Post_Plugin;
+
+// On a triangle, we suppose that there exists an interpolation scheme
+// such that u = \sum_i u_i \phi_i, phi_i being polynomials of order
+// p, i goes from 1...(p+1)(p+2)/2 and phi_i = \sum_j coeffs_{ij}
+// monomials_j and monomials are 1,x,y,x^2,xy,y^2,x^3,x^2y,xy^2,y^3...
+
+class adapt_point
+{
+public :
+ double x,y,z;
+ double X,Y,Z,val,valx,valy,valz;
+ double shape_functions[128];
+ static adapt_point * New (double x, double y, double z,
+ Double_Matrix *coeffs, Double_Matrix *eexps);
+ void print () const
+ {
+ printf("p %g %g\n" ,x,y);
+ }
+ bool operator < (const adapt_point & other) const
+ {
+ if(other.x < x) return true;
+ if(other.x > x) return false;
+ if(other.y < y) return true;
+ if(other.y > y) return false;
+ if(other.z < z) return true;
+ return false;
+ }
+ static std::set<adapt_point> all_points;
+};
+
+class adapt_edge
+{
+public:
+ adapt_edge (adapt_point *p1,adapt_point *p2)
+ : visible(false)
+ {
+ p[0] = p1;
+ p[1] = p2;
+ e[0] = e[1] = 0;
+ }
+ inline double V () const
+ {
+ return (p[0]->val + p[1]->val)/2.;
+ }
+ inline static void GSF (const double u, const double v, double w, double sf[])
+ {
+ sf[0] = (1-u)/2.;
+ sf[1] = (1+u)/2.;
+ }
+ void print ()
+ {
+ printf("p1 %g %g p2 %g %g\n", p[0]->x, p[0]->y, p[1]->x, p[1]->y);
+ }
+ static void Create (int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps) ;
+ static void Recur_Create (adapt_edge *e, int maxlevel, int level,
+ Double_Matrix *coeffs, Double_Matrix *eexps);
+ static void Error (double AVG , double tol );
+ static void Recur_Error (adapt_edge *e, double AVG, double tol);
+ bool visible;
+ adapt_point *p[2];
+ adapt_edge *e[2];
+ static std::list<adapt_edge*> all_elems;
+ static int nbNod;
+};
+
+class adapt_triangle
+{
+public:
+ adapt_triangle (adapt_point *p1, adapt_point *p2, adapt_point *p3)
+ : visible (false)
+ {
+ p[0] = p1;
+ p[1] = p2;
+ p[2] = p3;
+ e[0] = e[1] = e[2] = e[3] = 0;
+ }
+ inline double V () const
+ {
+ return (p[0]->val + p[1]->val + p[2]->val)/3.;
+ }
+ inline static void GSF (const double u, const double v, double w, double sf[])
+ {
+ sf[0] = 1. - u - v;
+ sf[1] = u;
+ sf[2] = v;
+ }
+ void print ()
+ {
+ printf("p1 %g %g p2 %g %g p3 %g %g\n",
+ p[0]->x, p[0]->y, p[1]->x, p[1]->y, p[2]->x, p[2]->y);
+ }
+ static void Create (int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps);
+ static void Recur_Create (adapt_triangle *t, int maxlevel, int level,
+ Double_Matrix *coeffs, Double_Matrix *eexps);
+ static void Error (double AVG, double tol);
+ static void Recur_Error (adapt_triangle *t, double AVG, double tol);
+ bool visible;
+ adapt_point *p[3];
+ adapt_triangle *e[4];
+ static std::list<adapt_triangle*> all_elems;
+ static int nbNod;
+};
+
+class adapt_quad
+{
+public:
+ adapt_quad (adapt_point *p1, adapt_point *p2, adapt_point *p3, adapt_point *p4)
+ : visible (false)
+ {
+ p[0] = p1;
+ p[1] = p2;
+ p[2] = p3;
+ p[3] = p4;
+ e[0] = e[1] = e[2] = e[3] = 0;
+ }
+ inline double V () const
+ {
+ return (p[0]->val + p[1]->val + p[2]->val + p[3]->val)/4.;
+ }
+ inline static void GSF (const double u, const double v, double w, double sf[])
+ {
+ sf[0] = 0.25 * (1. - u) * (1. - v);
+ sf[1] = 0.25 * (1. + u) * (1. - v);
+ sf[2] = 0.25 * (1. + u) * (1. + v);
+ sf[3] = 0.25 * (1. - u) * (1. + v);
+ }
+ void print ()
+ {
+ printf("p1 %g %g p2 %g %g p3 %g %g\n",
+ p[0]->x, p[0]->y, p[1]->x, p[1]->y, p[2]->x, p[2]->y);
+ }
+ static void Create (int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps);
+ static void Recur_Create (adapt_quad *q, int maxlevel, int level,
+ Double_Matrix *coeffs, Double_Matrix *eexps);
+ static void Error (double AVG, double tol);
+ static void Recur_Error (adapt_quad *q, double AVG, double tol);
+ bool visible;
+ adapt_point *p[4];
+ adapt_quad *e[4];
+ static std::list<adapt_quad*> all_elems;
+ static int nbNod;
+};
+
+class adapt_tet
+{
+public:
+ adapt_tet (adapt_point *p1, adapt_point *p2, adapt_point *p3, adapt_point *p4)
+ : visible (false)
+ {
+ p[0] = p1;
+ p[1] = p2;
+ p[2] = p3;
+ p[3] = p4;
+ e[0] = e[1] = e[2] = e[3] = 0;
+ e[4] = e[5] = e[6] = e[7] = 0;
+ }
+ inline static void GSF (const double u, const double v, double w, double sf[])
+ {
+ sf[0] = 1. - u - v - w;
+ sf[1] = u;
+ sf[2] = v;
+ sf[3] = w;
+ }
+ inline double V () const
+ {
+ return (p[0]->val + p[1]->val + p[2]->val + p[3]->val)/4.;
+ }
+ void print ()
+ {
+ printf("p1 %g %g p2 %g %g p3 %g %g\n",
+ p[0]->x, p[0]->y, p[1]->x, p[1]->y, p[2]->x, p[2]->y);
+ }
+ static void Create (int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps);
+ static void Recur_Create (adapt_tet *t, int maxlevel, int level,
+ Double_Matrix *coeffs, Double_Matrix *eexps);
+ static void Error (double AVG, double tol);
+ static void Recur_Error (adapt_tet *t, double AVG, double tol);
+ bool visible;
+ adapt_point *p[4];
+ adapt_tet *e[8];
+ static std::list<adapt_tet*> all_elems;
+ static int nbNod;
+};
+
+class adapt_hex
+{
+public:
+ adapt_hex (adapt_point *p1, adapt_point *p2, adapt_point *p3, adapt_point *p4,
+ adapt_point *p5, adapt_point *p6, adapt_point *p7, adapt_point *p8)
+ : visible (false)
+ {
+ p[0] = p1;
+ p[1] = p2;
+ p[2] = p3;
+ p[3] = p4;
+ p[4] = p5;
+ p[5] = p6;
+ p[6] = p7;
+ p[7] = p8;
+ e[0] = e[1] = e[2] = e[3] = 0;
+ e[4] = e[5] = e[6] = e[7] = 0;
+ }
+ inline static void GSF (const double u, const double v, double w, double sf[])
+ {
+ sf[0] = 0.125 * (1 - u) * (1 - v) * (1 - w);
+ sf[1] = 0.125 * (1 + u) * (1 - v) * (1 - w);
+ sf[2] = 0.125 * (1 + u) * (1 + v) * (1 - w);
+ sf[3] = 0.125 * (1 - u) * (1 + v) * (1 - w);
+ sf[4] = 0.125 * (1 - u) * (1 - v) * (1 + w);
+ sf[5] = 0.125 * (1 + u) * (1 - v) * (1 + w);
+ sf[6] = 0.125 * (1 + u) * (1 + v) * (1 + w);
+ sf[7] = 0.125 * (1 - u) * (1 + v) * (1 + w);
+ }
+ inline double V () const
+ {
+ return (p[0]->val + p[1]->val + p[2]->val+ p[3]->val +
+ p[4]->val + p[5]->val + p[6]->val+ p[7]->val)/8.;
+ }
+ void print ()
+ {
+ printf("p1 %g %g p2 %g %g p3 %g %g\n",
+ p[0]->x, p[0]->y, p[1]->x, p[1]->y, p[2]->x, p[2]->y);
+ }
+ static void Create (int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps);
+ static void Recur_Create (adapt_hex *h, int maxlevel, int level,
+ Double_Matrix *coeffs, Double_Matrix *eexps);
+ static void Error (double AVG, double tol);
+ static void Recur_Error (adapt_hex *h, double AVG, double tol);
+ bool visible;
+ adapt_point *p[8];
+ adapt_hex *e[8];
+ static std::list<adapt_hex*> all_elems;
+ static int nbNod;
+};
+
+class Adaptive_Post_View
+{
+ double tol;
+ double minval, maxval;
+ int presentZoomLevel;
+ double presentTol;
+ Double_Matrix * _eexps;
+ Double_Matrix * _coefs;
+ Double_Matrix * _coefsGeom;
+ Double_Matrix * _eexpsGeom;
+ Double_Matrix * _STposX;
+ Double_Matrix * _STposY;
+ Double_Matrix * _STposZ;
+ Double_Matrix * _STval;
+ // for vectors
+ Double_Matrix * _STvalX;
+ Double_Matrix * _STvalY;
+ Double_Matrix * _STvalZ;
+ Double_Matrix * _Interpolate;
+ Double_Matrix * _Geometry;
+public:
+ Adaptive_Post_View (Post_View *view, List_T *_coeffs, List_T *_eexps, List_T *_coeffsGeom=0, List_T *_eexpsGeom=0);
+ ~Adaptive_Post_View ();
+ int getGlobalResolutionLevel () const { return presentZoomLevel; }
+ void setGlobalResolutionLevel (Post_View * view, int level)
+ {
+ setAdaptiveResolutionLevel(view, level);
+ }
+ template <class ELEM>
+ void setAdaptiveResolutionLevel_TEMPL(Post_View * view, int level, int lemvelmax,
+ GMSH_Post_Plugin *plug, List_T **myList,
+ int *counter, int *done);
+ void setAdaptiveResolutionLevel (Post_View * view , int levelmax,
+ GMSH_Post_Plugin *plug = 0);
+ void initWithLowResolution (Post_View *view);
+ void setTolerance (const double eps) { tol=eps; }
+ double getTolerance () const { return tol; }
+ template <class ELEM>
+ int zoomElement (Post_View * view, int ielem, int level, int levelmax,
+ GMSH_Post_Plugin *plug, List_T *theList, int *counter);
+};
+
+template <class ELEM>
+void cleanElement ()
+{
+ typename std::list<ELEM*>::iterator it = ELEM::all_elems.begin();
+ typename std::list<ELEM*>::iterator ite = ELEM::all_elems.end();
+ for (; it != ite; ++it){
+ delete *it;
+ }
+ ELEM::all_elems.clear();
+ adapt_point::all_points.clear();
+}
+
+#endif
diff --git a/Common/ColorTable.cpp b/Common/ColorTable.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a0f0f9862811f494b52e60c5a99c2b78ce908a64
--- /dev/null
+++ b/Common/ColorTable.cpp
@@ -0,0 +1,464 @@
+// $Id: ColorTable.cpp,v 1.33 2006-11-27 22:22:07 geuzaine Exp $
+//
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// David Colignon
+//
+
+// These colortable routines were inspired by those provided in Vis5d,
+// a program for visualizing five dimensional gridded data sets
+// Copyright (C) 1990 - 1995 Bill Hibbard, Brian Paul, Dave Santek,
+// and Andre Battaiola.
+
+#include "Gmsh.h"
+#include "ColorTable.h"
+#include "Context.h"
+#include "Numeric.h"
+
+extern Context_T CTX;
+
+void ColorTable_InitParam(int number, GmshColorTable * ct)
+{
+ ct->size = 255;
+ for(int i = 0; i < COLORTABLE_NBMAX_PARAM; i++){
+ ct->ipar[i] = 0;
+ ct->dpar[i] = 0.;
+ }
+ ct->ipar[COLORTABLE_MODE] = COLORTABLE_RGB;
+ ct->ipar[COLORTABLE_NUMBER] = number;
+ ct->dpar[COLORTABLE_ALPHA] = 1.0;
+}
+
+static double gray(double s)
+{
+ return s < 0. ? 0. : (s < 1. ? s : 1.);
+}
+
+static double hot_r(double s)
+{
+ return s < 0. ? 0. : (s < 3./8. ? 8./3. * s : 1.);
+}
+
+static double hot_g(double s)
+{
+ return s < 3./8. ? 0. : (s < 6./8. ? 8./3. * (s - 3./8.) : 1.);
+}
+
+static double hot_b(double s)
+{
+ return s < 6./8. ? 0. : (s < 1. ? 8./2. * (s - 6./8.) : 1.);
+}
+
+static double cubic(double a, double b, double c, double d, double x)
+{
+ return a + b * x + c * x * x + d * x * x * x;
+}
+
+void ColorTable_Recompute(GmshColorTable * ct)
+{
+ double s, t, gamma;
+ int r, g, b, a;
+
+ double bias = ct->dpar[COLORTABLE_BIAS];
+ double curvature = ct->dpar[COLORTABLE_CURVATURE];
+ int rotation = ct->ipar[COLORTABLE_ROTATION];
+
+ for(int i = 0; i < ct->size; i++) {
+
+ if(ct->size > 1) {
+ if(i + rotation < 0)
+ s = (double)(i + rotation + ct->size) / (double)(ct->size - 1);
+ else if(i + rotation > ct->size - 1)
+ s = (double)(i + rotation - ct->size) / (double)(ct->size - 1);
+ else
+ s = (double)(i + rotation) / (double)(ct->size - 1);
+ }
+ else
+ s = 0.;
+
+ if(ct->ipar[COLORTABLE_SWAP])
+ s = 1.0 - s;
+
+ switch (ct->ipar[COLORTABLE_NUMBER]) {
+ case 0: // all black
+ r = g = b = 0;
+ break;
+ case 1: // vis5d
+ t = (curvature + 1.4) * (s - (1. + bias) / 2.);
+ r = (int)(128.0 + 127.0 * atan(7.0 * t) / 1.57);
+ g = (int)(128.0 + 127.0 * (2 * exp(-7 * t * t) - 1));
+ b = (int)(128.0 + 127.0 * atan(-7.0 * t) / 1.57);
+ break;
+ case 2: // matlab "jet"
+ {
+ double ii = (double)(s-bias)*128.;
+ if(ii < 0) ii = 0;
+ if(ii > 128) ii = 128;
+ double rr =
+ ii <= 46 ? 0. :
+ ii >= 111 ? -0.03125*(ii - 111) + 1. :
+ ii >= 78 ? 1. :
+ 0.03125*(ii - 46);
+ double gg =
+ ii <= 14 || ii >= 111 ? 0. :
+ ii >= 79 ? -0.03125*(ii - 111) :
+ ii <= 46 ? 0.03125*(ii - 14) :
+ 1.;
+ double bb =
+ ii >= 79 ? 0. :
+ ii >= 47 ? -0.03125*(ii - 79) :
+ ii <= 14 ? 0.03125*(ii - 14) + 1.:
+ 1.;
+ r = (int)(rr*255.);
+ g = (int)(gg*255.);
+ b = (int)(bb*255.);
+ }
+ break;
+ case 3: // lucie, samcef (?)
+ if(s - bias <= 0.) {
+ r = 0;
+ g = 0;
+ b = 255;
+ }
+ else if(s - bias <= 0.40) {
+ r = 0;
+ g = (int)((s - bias) * 637.5);
+ b = (int)(255. - (s - bias) * 637.5);
+ }
+ else if(s - bias <= 0.60) {
+ r = (int)(1275. * (s - bias - 0.4));
+ g = 255;
+ b = 0;
+ }
+ else if(s - bias <= 1.) {
+ r = 255;
+ g = (int)(255. - 637.5 * (s - bias - 0.6));
+ b = 0;
+ }
+ else {
+ r = 255;
+ g = 0;
+ b = 0;
+ }
+ break;
+ case 4: // rainbow
+ if(s - bias <= 0.) {
+ r = 0;
+ g = 0;
+ b = 255;
+ }
+ else if(s - bias <= 0.25 + curvature) {
+ curvature = (curvature == -0.25) ? -0.26 : curvature;
+ r = 0;
+ g = (int)((s - bias) * (255. / (0.25 + curvature)));
+ b = 255;
+ }
+ else if(s - bias <= 0.50) {
+ curvature = (curvature == 0.25) ? 0.26 : curvature;
+ r = 0;
+ g = 255;
+ b = (int)(255. - (255. / (0.25 - curvature)) * (s - bias - 0.25 - curvature));
+ }
+ else if(s - bias <= 0.75 - curvature) {
+ curvature = (curvature == 0.25) ? 0.26 : curvature;
+ r = (int)((s - bias - 0.5) * (255. / (0.25 - curvature)));
+ g = 255;
+ b = 0;
+ }
+ else if(s - bias <= 1.) {
+ curvature = (curvature == -0.25) ? -0.26 : curvature;
+ r = 255;
+ g = (int)(255. - (255. / (0.25 + curvature)) * (s - bias - 0.75 + curvature));
+ b = 0;
+ }
+ else {
+ r = 255;
+ g = 0;
+ b = 0;
+ }
+ break;
+ case 5: // emc2000 (rainbow with black and white)
+ if(s - bias <= 0.) {
+ r = 0;
+ g = 0;
+ b = 0;
+ }
+ else if(s - bias <= 0.2) {
+ r = (int)(57 * (1 - 100 * ((s - bias) - 0.1) * ((s - bias) - 0.1)));
+ g = 0;
+ b = (int)((s - bias) * (255. / 0.2));
+ }
+ else if(s - bias <= 0.3624) {
+ r = 0;
+ g = (int)((s - bias - 0.2) * (255. / 0.1624));
+ b = 255;
+ }
+ else if(s - bias <= 0.50) {
+ r = 0;
+ g = 255;
+ b = (int)(255. - (255. / 0.1376) * (s - bias - 0.3624));
+ }
+ else if(s - bias <= 0.6376) {
+ r = (int)((s - bias - 0.5) * (255. / 0.1376));
+ g = 255;
+ b = 0;
+ }
+ else if(s - bias <= 0.8) {
+ r = 255;
+ g = (int)(255. - (255. / 0.1624) * (s - bias - 0.6376));
+ b = 0;
+ }
+ else if(s - bias <= 1.0) {
+ r = 255;
+ g = (int)((255. / 0.2) * (s - bias - 0.8));
+ b = (int)(-3187.66 * (s - bias) * (s - bias) + 7012.76 * (s - bias) - 3570.61);
+ }
+ else {
+ r = 255;
+ g = 255;
+ b = 255;
+ }
+ break;
+ case 6: // darkblue->red->yellow->white
+ r = (int)(255. * cubic(-0.0506169, 2.81633, -1.87033, 0.0524573, s-bias));
+ g = (int)(255. * cubic(0.0485868, -1.26109, 6.3074, -4.12498, s-bias));
+ b = (int)(255. * cubic(0.364662, 1.50814, -7.36756, 6.51847, s-bias));
+ break;
+ case 7: // matlab "hot"
+ r = (int)(255. * hot_r(s-bias));
+ g = (int)(255. * hot_g(s-bias));
+ b = (int)(255. * hot_b(s-bias));
+ break;
+ case 8: // matlab "pink"
+ r = (int)(255. * sqrt((2.*gray(s-bias) + hot_r(s-bias))/3.));
+ g = (int)(255. * sqrt((2.*gray(s-bias) + hot_g(s-bias))/3.));
+ b = (int)(255. * sqrt((2.*gray(s-bias) + hot_b(s-bias))/3.));
+ break;
+ case 9: // grayscale
+ if(s - bias <= 0.) {
+ r = g = b = 0;
+ }
+ else if(s - bias <= 1.) {
+ r = g = b = (int)(255 * (1. - curvature) * (s - bias));
+ }
+ else {
+ r = g = b = (int)(255 * (1. - curvature));
+ }
+ break;
+ case 10: // all white
+ r = g = b = 255;
+ break;
+ case 11: // matlab "hsv"
+ {
+ double H = 6. * s + 1.e-10, R, G, B;
+ HSV_to_RGB(H, 1., 1., &R, &G, &B);
+ r = (int)(255 * R);
+ g = (int)(255 * G);
+ b = (int)(255 * B);
+ }
+ break;
+ case 12: // spectrum (truncated hsv)
+ {
+ double H = 5. * s + 1.e-10, R, G, B;
+ HSV_to_RGB(H, 1., 1., &R, &G, &B);
+ r = (int)(255 * R);
+ g = (int)(255 * G);
+ b = (int)(255 * B);
+ }
+ break;
+ case 13: // matlab "bone"
+ r = (int)(255. * (7.*gray(s-bias) + hot_b(s-bias))/8.);
+ g = (int)(255. * (7.*gray(s-bias) + hot_g(s-bias))/8.);
+ b = (int)(255. * (7.*gray(s-bias) + hot_r(s-bias))/8.);
+ break;
+ case 14: // matlab "spring"
+ r = (int)(255. * 1.);
+ g = (int)(255. * gray(s-bias));
+ b = (int)(255. * (1. - gray(s-bias)));
+ break;
+ case 15: // matlab "summer"
+ r = (int)(255. * gray(s-bias));
+ g = (int)(255. * (0.5+gray(s-bias)/2.));
+ b = (int)(255. * 0.4);
+ break;
+ case 16: // matlab "autumn"
+ r = (int)(255. * 1.);
+ g = (int)(255. * gray(s-bias));
+ b = (int)(255. * 0.);
+ break;
+ case 17: // matlab "winter"
+ r = (int)(255. * 0.);
+ g = (int)(255. * gray(s-bias));
+ b = (int)(255. * (0.5+(1.-gray(s-bias))/2.));
+ break;
+ case 18: // matlab "cool"
+ r = (int)(255. * gray(s-bias));
+ g = (int)(255. * (1.-gray(s-bias)));
+ b = (int)(255. * 1.);
+ break;
+ case 19: // matlab "copper"
+ r = (int)(255. * DMIN(1., gray(s-bias) * 1.25));
+ g = (int)(255. * DMIN(1., gray(s-bias) * 0.7812));
+ b = (int)(255. * DMIN(1., gray(s-bias) * 0.4975));
+ break;
+ default:
+ r = g = b = 0;
+ break;
+ }
+
+ double aa = 1.0;
+ if(ct->dpar[COLORTABLE_ALPHAPOW])
+ aa = pow(s ? s : 1.e-10, ct->dpar[COLORTABLE_ALPHAPOW]);
+ a = (int)(255. * aa * ct->dpar[COLORTABLE_ALPHA]);
+
+ if(ct->dpar[COLORTABLE_BETA]) {
+ if(ct->dpar[COLORTABLE_BETA] > 0.0)
+ gamma = 1. - ct->dpar[COLORTABLE_BETA];
+ else
+ gamma = 1. / (1.001 + ct->dpar[COLORTABLE_BETA]); // beta is thresholded to [-1,1]
+ r = (int)(255. * pow((double)r / 255., gamma));
+ g = (int)(255. * pow((double)g / 255., gamma));
+ b = (int)(255. * pow((double)b / 255., gamma));
+ }
+
+ if(ct->ipar[COLORTABLE_INVERT]) {
+ r = 255 - r;
+ g = 255 - g;
+ b = 255 - b;
+ }
+
+ // clamp to [0,255]
+ r = r < 0 ? 0 : (r > 255 ? 255 : r);
+ g = g < 0 ? 0 : (g > 255 ? 255 : g);
+ b = b < 0 ? 0 : (b > 255 ? 255 : b);
+ a = a < 0 ? 0 : (a > 255 ? 255 : a);
+
+ ct->table[i] = CTX.PACK_COLOR(r, g, b, a);
+ }
+
+}
+
+static GmshColorTable clip;
+
+void ColorTable_Copy(GmshColorTable * ct)
+{
+ clip.size = ct->size;
+ memcpy(clip.table, ct->table, ct->size * sizeof(unsigned int));
+ memcpy(clip.ipar, ct->ipar, COLORTABLE_NBMAX_PARAM * sizeof(int));
+ memcpy(clip.dpar, ct->dpar, COLORTABLE_NBMAX_PARAM * sizeof(double));
+}
+
+void ColorTable_Paste(GmshColorTable * ct)
+{
+ ct->size = clip.size;
+ memcpy(ct->table, clip.table, clip.size * sizeof(unsigned int));
+ memcpy(ct->ipar, clip.ipar, COLORTABLE_NBMAX_PARAM * sizeof(int));
+ memcpy(ct->dpar, clip.dpar, COLORTABLE_NBMAX_PARAM * sizeof(double));
+}
+
+int ColorTable_Diff(GmshColorTable * ct1, GmshColorTable * ct2)
+{
+ if(ct1->size != ct2->size) return 1;
+ for(int i = 0; i < ct1->size; i++)
+ if(ct1->table[i] != ct2->table[i]) return 1;
+ return 0;
+}
+
+void ColorTable_Print(GmshColorTable * ct, FILE * fp)
+{
+ int i, r, g, b, a;
+ char tmp1[1024], tmp2[1024];
+
+ strcpy(tmp1, "");
+ for(i = 0; i < ct->size; i++) {
+ r = CTX.UNPACK_RED(ct->table[i]);
+ g = CTX.UNPACK_GREEN(ct->table[i]);
+ b = CTX.UNPACK_BLUE(ct->table[i]);
+ a = CTX.UNPACK_ALPHA(ct->table[i]);
+ if(i && !(i % 4)) {
+ if(fp)
+ fprintf(fp, "%s\n", tmp1);
+ else
+ Msg(DIRECT, tmp1);
+ strcpy(tmp1, "");
+ }
+ sprintf(tmp2, "{%d, %d, %d, %d}", r, g, b, a);
+ strcat(tmp1, tmp2);
+ if(i != ct->size - 1)
+ strcat(tmp1, ", ");
+ }
+ if(fp)
+ fprintf(fp, "%s\n", tmp1);
+ else
+ Msg(DIRECT, tmp1);
+}
+
+int ColorTable_IsAlpha(GmshColorTable * ct)
+{
+ int i, a;
+ for(i = 0; i < ct->size; i++) {
+ a = CTX.UNPACK_ALPHA(ct->table[i]);
+ if(a < 255)
+ return 1;
+ }
+ return 0;
+}
+
+// HSV/RBG conversion routines
+
+void HSV_to_RGB(double H, double S, double V,
+ double *R, double *G, double *B)
+{
+ if(S < 5.0e-6) {
+ *R = *G = *B = V;
+ }
+ else {
+ int i = (int)H;
+ double f = H - (float)i;
+ double p1 = V*(1.0-S);
+ double p2 = V*(1.0-S*f);
+ double p3 = V*(1.0-S*(1.0-f));
+ switch(i){
+ case 0: *R = V; *G = p3; *B = p1; break;
+ case 1: *R = p2; *G = V; *B = p1; break;
+ case 2: *R = p1; *G = V; *B = p3; break;
+ case 3: *R = p1; *G = p2; *B = V; break;
+ case 4: *R = p3; *G = p1; *B = V; break;
+ case 5: *R = V; *G = p1; *B = p2; break;
+ }
+ }
+}
+
+void RGB_to_HSV(double R, double G, double B,
+ double *H, double *S, double *V)
+{
+ double maxv = R > G ? R : G; if(B > maxv) maxv = B;
+ *V = maxv;
+ if(maxv>0){
+ double minv = R < G ? R : G; if(B < minv) minv = B;
+ *S = 1.0 - double(minv)/maxv;
+ if(maxv > minv){
+ if(maxv == R){ *H = (G-B)/double(maxv-minv); if (*H<0) *H += 6.0; }
+ else if(maxv == G) *H = 2.0+(B-R)/double(maxv-minv);
+ else *H = 4.0+(R-G)/double(maxv-minv);
+ }
+ }
+}
diff --git a/Common/ColorTable.h b/Common/ColorTable.h
new file mode 100644
index 0000000000000000000000000000000000000000..11bde849ec16469597f95aaaadeff0657c139dd5
--- /dev/null
+++ b/Common/ColorTable.h
@@ -0,0 +1,67 @@
+#ifndef _COLORTABLE_H_
+#define _COLORTABLE_H_
+
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#define COLORTABLE_NBMAX_PARAM 10
+#define COLORTABLE_NBMAX_COLOR 1024
+
+typedef struct{
+ unsigned int table[COLORTABLE_NBMAX_COLOR];
+ int size; // must be >= 2
+ int ipar[COLORTABLE_NBMAX_PARAM];
+ double dpar[COLORTABLE_NBMAX_PARAM];
+}GmshColorTable;
+
+// COLORTABLE_MODE
+
+#define COLORTABLE_RGB 1
+#define COLORTABLE_HSV 2
+
+// integrer parameters indices
+
+#define COLORTABLE_NUMBER 0 // predefined curve index
+#define COLORTABLE_INVERT 1 // invert (rbg<->255-rgb)
+#define COLORTABLE_SWAP 2 // swap (min<->max)
+#define COLORTABLE_ROTATION 3 // rotation
+#define COLORTABLE_MODE 4 // mode (rgb, hsv)
+
+// double parameters indices
+
+#define COLORTABLE_CURVATURE 0 // curvature
+#define COLORTABLE_BIAS 1 // offset
+#define COLORTABLE_ALPHA 2 // alpha channel value
+#define COLORTABLE_BETA 3 // beta coeff for brighten
+#define COLORTABLE_ALPHAPOW 4 // alpha channel power value
+
+void ColorTable_InitParam(int number, GmshColorTable *ct);
+void ColorTable_Recompute(GmshColorTable *ct);
+void ColorTable_Copy(GmshColorTable *ct);
+void ColorTable_Paste(GmshColorTable *ct);
+void ColorTable_Print(GmshColorTable *ct, FILE *fp) ;
+int ColorTable_IsAlpha(GmshColorTable *ct) ;
+int ColorTable_Diff(GmshColorTable *ct1, GmshColorTable *ct2);
+
+void RGB_to_HSV(double R, double G, double B,
+ double *H, double *S, double *V);
+void HSV_to_RGB(double H, double S, double V,
+ double *R, double *G, double *B);
+
+#endif
diff --git a/Common/OctreePost.cpp b/Common/OctreePost.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c8f6b61e6ba1bf592317ed22590d924430245215
--- /dev/null
+++ b/Common/OctreePost.cpp
@@ -0,0 +1,405 @@
+// $Id: OctreePost.cpp,v 1.2 2006-11-27 22:22:07 geuzaine Exp $
+//
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "Octree.h"
+#include "OctreePost.h"
+#include "List.h"
+#include "Views.h"
+#include "Numeric.h"
+#include "Message.h"
+#include "ShapeFunctions.h"
+
+static void minmax(int n, double *X, double *Y, double *Z,
+ double *min, double *max)
+{
+ min[0] = X[0];
+ min[1] = Y[0];
+ min[2] = Z[0];
+ max[0] = X[0];
+ max[1] = Y[0];
+ max[2] = Z[0];
+
+ for(int i = 1; i < n; i++) {
+ min[0] = (X[i] < min[0]) ? X[i] : min[0];
+ min[1] = (Y[i] < min[1]) ? Y[i] : min[1];
+ min[2] = (Z[i] < min[2]) ? Z[i] : min[2];
+ max[0] = (X[i] > max[0]) ? X[i] : max[0];
+ max[1] = (Y[i] > max[1]) ? Y[i] : max[1];
+ max[2] = (Z[i] > max[2]) ? Z[i] : max[2];
+ }
+}
+
+static void centroid(int n, double *X, double *Y, double *Z, double *c)
+{
+ const double oc = 1./(double)n;
+ c[0] = X[0];
+ c[1] = Y[0];
+ c[2] = Z[0];
+ for(int i = 1; i < n; i++) {
+ c[0] += X[i];
+ c[1] += Y[i];
+ c[2] += Z[i];
+ }
+ c[0] *= oc;
+ c[1] *= oc;
+ c[2] *= oc;
+}
+
+void linBB(void *a, double *min, double *max)
+{
+ double *X = (double*) a, *Y = &X[2], *Z = &X[4];
+ minmax(2, X, Y, Z, min, max);
+}
+
+void triBB(void *a, double *min, double *max)
+{
+ double *X = (double*) a, *Y = &X[3], *Z = &X[6];
+ minmax(3, X, Y, Z, min, max);
+}
+
+void quaBB(void *a, double *min, double *max)
+{
+ double *X = (double*) a, *Y = &X[4], *Z = &X[8];
+ minmax(4, X, Y, Z, min, max);
+}
+
+void tetBB(void *a, double *min, double *max)
+{
+ double *X = (double*) a, *Y = &X[4], *Z = &X[8];
+ minmax(4, X, Y, Z, min, max);
+}
+
+void hexBB(void *a, double *min, double *max)
+{
+ double *X = (double*) a, *Y = &X[8], *Z = &X[16];
+ minmax(8, X, Y, Z, min, max);
+}
+
+void priBB(void *a, double *min, double *max)
+{
+ double *X = (double*) a, *Y = &X[6], *Z = &X[12];
+ minmax(6, X, Y, Z, min, max);
+}
+
+void pyrBB(void *a, double *min, double *max)
+{
+ double *X = (double*) a, *Y = &X[5], *Z = &X[10];
+ minmax(5, X, Y, Z, min, max);
+}
+
+int linInEle(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[2], *Z = &X[4], uvw[3];
+ line lin(X, Y, Z);
+ lin.xyz2uvw(x, uvw);
+ return lin.isInside(uvw[0], uvw[1], uvw[2]);
+}
+
+int triInEle(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[3], *Z = &X[6], uvw[3];
+ triangle tri(X, Y, Z);
+ tri.xyz2uvw(x, uvw);
+ return tri.isInside(uvw[0], uvw[1], uvw[2]);
+}
+
+int quaInEle(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[4], *Z = &X[8], uvw[3];
+ quadrangle qua(X, Y, Z);
+ qua.xyz2uvw(x, uvw);
+ return qua.isInside(uvw[0], uvw[1], uvw[2]);
+}
+
+int tetInEle(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[4], *Z = &X[8], uvw[3];
+ tetrahedron tet(X, Y, Z);
+ tet.xyz2uvw(x, uvw);
+ return tet.isInside(uvw[0], uvw[1], uvw[2]);
+}
+
+int hexInEle(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[8], *Z = &X[16], uvw[3];
+ hexahedron hex(X, Y, Z);
+ hex.xyz2uvw(x, uvw);
+ return hex.isInside(uvw[0], uvw[1], uvw[2]);
+}
+
+int priInEle(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[6], *Z = &X[12], uvw[3];
+ prism pri(X, Y, Z);
+ pri.xyz2uvw(x, uvw);
+ return pri.isInside(uvw[0], uvw[1], uvw[2]);
+}
+
+int pyrInEle(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[5], *Z = &X[10], uvw[3];
+ pyramid pyr(X, Y, Z);
+ pyr.xyz2uvw(x, uvw);
+ return pyr.isInside(uvw[0], uvw[1], uvw[2]);
+}
+
+void linCentroid(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[2], *Z = &X[4];
+ centroid(2, X, Y, Z, x);
+}
+
+void triCentroid(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[3], *Z = &X[6];
+ centroid(3, X, Y, Z, x);
+}
+
+void quaCentroid(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[4], *Z = &X[8];
+ centroid(4, X, Y, Z, x);
+}
+
+void tetCentroid(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[4], *Z = &X[8];
+ centroid(4, X, Y, Z, x);
+}
+
+void hexCentroid(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[8], *Z = &X[16];
+ centroid(8, X, Y, Z, x);
+}
+
+void priCentroid(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[6], *Z = &X[12];
+ centroid(6, X, Y, Z, x);
+}
+
+void pyrCentroid(void *a, double *x)
+{
+ double *X = (double*) a, *Y = &X[5], *Z = &X[10];
+ centroid(5, X, Y, Z, x);
+}
+
+static void addListOfStuff(Octree *o, List_T *l, int nbelm)
+{
+ if(!l) return;
+
+ for(int i = 0; i < List_Nbr(l); i += nbelm){
+ double * X = (double *)List_Pointer_Fast(l, i);
+ Octree_Insert(X, o);
+ }
+}
+
+OctreePost::~OctreePost()
+{
+ Octree_Delete(ST); Octree_Delete(VT); Octree_Delete(TT);
+ Octree_Delete(SQ); Octree_Delete(VQ); Octree_Delete(TQ);
+ Octree_Delete(SS); Octree_Delete(VS); Octree_Delete(TS);
+ Octree_Delete(SH); Octree_Delete(VH); Octree_Delete(TH);
+ Octree_Delete(SI); Octree_Delete(VI); Octree_Delete(TI);
+ Octree_Delete(SY); Octree_Delete(VY); Octree_Delete(TY);
+}
+
+OctreePost::OctreePost(Post_View *v)
+ : theView(v)
+{
+ double min [3] = {v->BBox[0], v->BBox[2], v->BBox[4]};
+
+ double size[3] = {v->BBox[1] - v->BBox[0],
+ v->BBox[3] - v->BBox[2],
+ v->BBox[5] - v->BBox[4]};
+
+ const int maxElePerBucket = 100; // memory vs. speed trade-off
+
+ SL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
+ addListOfStuff(SL, v->SL, 6 + 2 * v->NbTimeStep);
+ Octree_Arrange(SL);
+ VL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
+ addListOfStuff(VL, v->VL, 6 + 6 * v->NbTimeStep);
+ Octree_Arrange(VL);
+ TL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
+ addListOfStuff(TL, v->TL, 6 + 18 * v->NbTimeStep);
+ Octree_Arrange(TL);
+
+ ST = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
+ addListOfStuff(ST, v->ST, 9 + 3 * v->NbTimeStep);
+ Octree_Arrange(ST);
+ VT = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
+ addListOfStuff(VT, v->VT, 9 + 9 * v->NbTimeStep);
+ Octree_Arrange(VT);
+ TT = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
+ addListOfStuff(TT, v->TT, 9 + 27 * v->NbTimeStep);
+ Octree_Arrange(TT);
+
+ SQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
+ addListOfStuff(SQ, v->SQ, 12 + 4 * v->NbTimeStep);
+ Octree_Arrange(SQ);
+ VQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
+ addListOfStuff(VQ, v->VQ, 12 + 12 * v->NbTimeStep);
+ Octree_Arrange(VQ);
+ TQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
+ addListOfStuff(TQ, v->TQ, 12 + 36 * v->NbTimeStep);
+ Octree_Arrange(TQ);
+
+ SS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
+ addListOfStuff(SS, v->SS, 12 + 4 * v->NbTimeStep);
+ Octree_Arrange(SS);
+ VS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
+ addListOfStuff(VS, v->VS, 12 + 12 * v->NbTimeStep);
+ Octree_Arrange(VS);
+ TS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
+ addListOfStuff(TS, v->TS, 12 + 36 * v->NbTimeStep);
+ Octree_Arrange(TS);
+
+ SH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
+ addListOfStuff(SH, v->SH, 24 + 8 * v->NbTimeStep);
+ Octree_Arrange(SH);
+ VH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
+ addListOfStuff(VH, v->VH, 24 + 24 * v->NbTimeStep);
+ Octree_Arrange(VH);
+ TH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
+ addListOfStuff(TH, v->TH, 24 + 72 * v->NbTimeStep);
+ Octree_Arrange(TH);
+
+ SI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
+ addListOfStuff(SI, v->SI, 18 + 6 * v->NbTimeStep);
+ Octree_Arrange(SI);
+ VI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
+ addListOfStuff(VI, v->VI, 18 + 18 * v->NbTimeStep);
+ Octree_Arrange(VI);
+ TI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
+ addListOfStuff(TI, v->TI, 18 + 54 * v->NbTimeStep);
+ Octree_Arrange(TI);
+
+ SY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
+ addListOfStuff(SY, v->SY, 15 + 5 * v->NbTimeStep);
+ Octree_Arrange(SY);
+ VY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
+ addListOfStuff(VY, v->VY, 15 + 15 * v->NbTimeStep);
+ Octree_Arrange(VY);
+ TY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
+ addListOfStuff(TY, v->TY, 15 + 45 * v->NbTimeStep);
+ Octree_Arrange(TY);
+}
+
+bool OctreePost::getValue(void *in, int dim, int nbNod, int nbComp,
+ double P[3], int timestep, double *values, double *size_elem)
+{
+ if(!in) return false;
+
+ double *X = (double*)in, *Y = &X[nbNod], *Z = &X[2*nbNod], *V = &X[3*nbNod], U[3];
+
+ elementFactory factory;
+ element *e = factory.create(nbNod, dim, X, Y, Z);
+ if(!e) return false;
+
+ e->xyz2uvw(P, U);
+ if(timestep < 0){
+ for(int i = 0; i < theView->NbTimeStep; i++)
+ for(int j = 0; j < nbComp; j++)
+ values[nbComp * i + j] = e->interpolate(&V[nbNod * nbComp * i + j],
+ U[0], U[1], U[2], nbComp);
+ }
+ else{
+ for(int j = 0; j < nbComp; j++)
+ values[j] = e->interpolate(&V[nbNod * nbComp * timestep + j],
+ U[0], U[1], U[2], nbComp);
+ }
+
+ if(size_elem)
+ *size_elem = e->maxEdgeLength();
+
+ delete e;
+ return true;
+}
+
+bool OctreePost::searchScalar(double x, double y, double z, double *values,
+ int timestep, double *size_elem)
+{
+ double P[3] = {x, y, z};
+
+ if(timestep < 0)
+ for(int i = 0; i < theView->NbTimeStep; i++)
+ values[i] = 0.0;
+ else
+ values[0] = 0.0;
+
+ if(getValue(Octree_Search(P, SS), 3, 4, 1, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, SH), 3, 8, 1, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, SI), 3, 6, 1, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, SY), 3, 5, 1, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, ST), 2, 3, 1, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, SQ), 2, 4, 1, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, SL), 1, 2, 1, P, timestep, values, size_elem)) return true;
+
+ return false;
+}
+
+bool OctreePost::searchVector(double x, double y, double z, double *values,
+ int timestep, double *size_elem)
+{
+ double P[3] = {x, y, z};
+
+ if(timestep < 0)
+ for(int i = 0; i < 3 * theView->NbTimeStep; i++)
+ values[i] = 0.0;
+ else
+ for(int i = 0; i < 3; i++)
+ values[i] = 0.0;
+
+ if(getValue(Octree_Search(P, VS), 3, 4, 3, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, VH), 3, 8, 3, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, VI), 3, 6, 3, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, VY), 3, 5, 3, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, VT), 2, 3, 3, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, VQ), 2, 4, 3, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, VL), 1, 2, 3, P, timestep, values, size_elem)) return true;
+
+ return false;
+}
+
+bool OctreePost::searchTensor(double x, double y, double z, double *values,
+ int timestep, double *size_elem)
+{
+ double P[3] = {x, y, z};
+
+ if(timestep < 0)
+ for(int i = 0; i < 9 * theView->NbTimeStep; i++)
+ values[i] = 0.0;
+ else
+ for(int i = 0; i < 9; i++)
+ values[i] = 0.0;
+
+ if(getValue(Octree_Search(P, TS), 3, 4, 9, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, TH), 3, 8, 9, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, TI), 3, 6, 9, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, TY), 3, 5, 9, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, TT), 2, 3, 9, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, TQ), 2, 4, 9, P, timestep, values, size_elem)) return true;
+ if(getValue(Octree_Search(P, TL), 1, 2, 9, P, timestep, values, size_elem)) return true;
+
+ return false;
+}
diff --git a/Common/OctreePost.h b/Common/OctreePost.h
new file mode 100644
index 0000000000000000000000000000000000000000..e7960d984f14026a2d5e86f7377d7b5c37bf866c
--- /dev/null
+++ b/Common/OctreePost.h
@@ -0,0 +1,55 @@
+#ifndef _OCTREE_POST_H_
+#define _OCTREE_POST_H_
+
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "Octree.h"
+
+class Post_View;
+
+class OctreePost
+{
+ Octree *SL, *VL, *TL;
+ Octree *ST, *VT, *TT;
+ Octree *SQ, *VQ, *TQ;
+ Octree *SS, *VS, *TS;
+ Octree *SH, *VH, *TH;
+ Octree *SI, *VI, *TI;
+ Octree *SY, *VY, *TY;
+ Post_View *theView;
+ bool getValue(void *in, int dim, int nbNod, int nbComp,
+ double P[3], int timestep, double *values, double *size_elem);
+ public :
+ OctreePost(Post_View *);
+ ~OctreePost();
+ // search for the value of the View at point x, y, z. Values are
+ // interpolated using standard first order shape functions in the
+ // post element. If several time steps are present, they are all
+ // interpolated unless time step is set to a different value than
+ // -1.
+ bool searchScalar(double x, double y, double z, double *values,
+ int timestep = -1, double *size_elem = 0);
+ bool searchVector(double x, double y, double z, double *values,
+ int timestep = -1, double *size_elem = 0);
+ bool searchTensor(double x, double y, double z, double *values,
+ int timestep = -1, double *size_elem = 0);
+};
+
+#endif
diff --git a/Common/Views.cpp b/Common/Views.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7e77b7f8f1bd1405e86ac78a12c839962b4975d9
--- /dev/null
+++ b/Common/Views.cpp
@@ -0,0 +1,1489 @@
+// $Id: Views.cpp,v 1.195 2007-04-21 19:39:59 geuzaine Exp $
+//
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Nicolas Tardieu
+
+#include <set>
+#include "Gmsh.h"
+#include "Numeric.h"
+#include "Views.h"
+#include "Context.h"
+#include "Options.h"
+#include "ColorTable.h"
+#include "SmoothData.h"
+#include "BackgroundMesh.h"
+
+#if defined(HAVE_MATH_EVAL)
+#include "matheval.h"
+#endif
+
+extern Context_T CTX;
+
+#if defined(HAVE_FLTK)
+void UpdateViewsInGUI();
+#endif
+
+// Static reference view
+
+Post_View *Post_ViewReference = NULL;
+
+// FIXME: the whole View interface should be rewritten in C++...
+
+int fcmpPostViewNum(const void *v1, const void *v2)
+{
+ return ((*(Post_View **) v1)->Num - (*(Post_View **) v2)->Num);
+}
+
+int fcmpPostViewAliasOf(const void *v1, const void *v2)
+{
+ return ((*(Post_View **) v1)->AliasOf - (*(Post_View **) v2)->AliasOf);
+}
+
+int fcmpPostViewName(const void *v1, const void *v2)
+{
+ return strcmp((*(Post_View **) v1)->Name, (*(Post_View **) v2)->Name);
+}
+
+int fcmpPostViewVisibility(const void *v1, const void *v2)
+{
+ return ((*(Post_View **) v2)->Visible - (*(Post_View **) v1)->Visible);
+}
+
+Post_View *BeginView(int allocate)
+{
+ Post_View *v = (Post_View *)Malloc(sizeof(Post_View));
+ static int UniqueNum = 0;
+
+ if(!CTX.post.list)
+ CTX.post.list = List_Create(100, 100, sizeof(Post_View *));
+
+ // Important notes:
+ // - each view *must* have a unique number
+ // - the view list is assumned to be sorted with increasing nums
+
+ if(!CTX.post.force_num) {
+ v->Num = ++UniqueNum;
+ List_Add(CTX.post.list, &v);
+ }
+ else {
+ v->Num = CTX.post.force_num;
+ List_Replace(CTX.post.list, &v, fcmpPostViewNum);
+ // FIXME: need to check here if the old view is used as a field in
+ // a background mesh (and if it is, remove that field). Until we
+ // do this, let's just invalidate all the lc fields
+ BGMReset();
+ }
+
+ int i = List_ISearch(CTX.post.list, &v, fcmpPostViewNum);
+ List_Read(CTX.post.list, i, &v);
+
+ v->Index = i;
+ v->Dirty = 1;
+ v->TriVertexArray = NULL;
+ v->LinVertexArray = NULL;
+ v->NbTimeStep = 1;
+ v->TimeStepMin = NULL;
+ v->TimeStepMax = NULL;
+ v->NbSP = v->NbVP = v->NbTP = 0;
+ v->NbSL = v->NbVL = v->NbTL = v->NbSL2 = v->NbVL2 = v->NbTL2 = 0;
+ v->NbST = v->NbVT = v->NbTT = v->NbST2 = v->NbVT2 = v->NbTT2 = 0;
+ v->NbSQ = v->NbVQ = v->NbTQ = v->NbSQ2 = v->NbVQ2 = v->NbTQ2 = 0;
+ v->NbSS = v->NbVS = v->NbTS = v->NbSS2 = v->NbVS2 = v->NbTS2 = 0;
+ v->NbSH = v->NbVH = v->NbTH = v->NbSH2 = v->NbVH2 = v->NbTH2 = 0;
+ v->NbSI = v->NbVI = v->NbTI = v->NbSI2 = v->NbVI2 = v->NbTI2 = 0;
+ v->NbSY = v->NbVY = v->NbTY = v->NbSY2 = v->NbVY2 = v->NbTY2 = 0;
+ v->NbT2 = v->NbT3 = 0;
+
+ if(allocate) {
+ v->DataSize = sizeof(double);
+ v->Grains = new std::map < int , List_T * > ;
+ v->DisplayListsOfGrains= new std::map < int , int > ;
+
+#define LCD List_Create(1, 1000, sizeof(double))
+ v->Time = LCD;
+ v->SP = LCD; v->VP = LCD; v->TP = LCD;
+ v->SL = LCD; v->VL = LCD; v->TL = LCD; v->SL2 = LCD; v->VL2 = LCD; v->TL2 = LCD;
+ v->ST = LCD; v->VT = LCD; v->TT = LCD; v->ST2 = LCD; v->VT2 = LCD; v->TT2 = LCD;
+ v->SQ = LCD; v->VQ = LCD; v->TQ = LCD; v->SQ2 = LCD; v->VQ2 = LCD; v->TQ2 = LCD;
+ v->SS = LCD; v->VS = LCD; v->TS = LCD; v->SS2 = LCD; v->VS2 = LCD; v->TS2 = LCD;
+ v->SH = LCD; v->VH = LCD; v->TH = LCD; v->SH2 = LCD; v->VH2 = LCD; v->TH2 = LCD;
+ v->SI = LCD; v->VI = LCD; v->TI = LCD; v->SI2 = LCD; v->VI2 = LCD; v->TI2 = LCD;
+ v->SY = LCD; v->VY = LCD; v->TY = LCD; v->SY2 = LCD; v->VY2 = LCD; v->TY2 = LCD;
+#undef LCD
+
+ v->T2D = List_Create(1, 100, sizeof(double));
+ v->T2C = List_Create(1, 100, sizeof(char));
+ v->T3D = List_Create(1, 100, sizeof(double));
+ v->T3C = List_Create(1, 100, sizeof(char));
+ }
+ else {
+ v->Time = NULL;
+ v->Grains = 0;
+ v->DisplayListsOfGrains = 0;
+ v->SP = v->VP = v->TP = NULL;
+ v->SL = v->VL = v->TL = v->SL2 = v->VL2 = v->TL2 = NULL;
+ v->ST = v->VT = v->TT = v->ST2 = v->VT2 = v->TT2 = NULL;
+ v->SQ = v->VQ = v->TQ = v->SQ2 = v->VQ2 = v->TQ2 = NULL;
+ v->SS = v->VS = v->TS = v->SS2 = v->VS2 = v->TS2 = NULL;
+ v->SH = v->VH = v->TH = v->SH2 = v->VH2 = v->TH2 = NULL;
+ v->SI = v->VI = v->TI = v->SI2 = v->VI2 = v->TI2 = NULL;
+ v->SY = v->VY = v->TY = v->SY2 = v->VY2 = v->TY2 = NULL;
+ v->T2D = v->T2C = NULL;
+ v->T3D = v->T3C = NULL;
+ }
+
+ // Copy all options from the reference view initialized in InitOptions()
+ CopyViewOptions(Post_ViewReference, v);
+
+ v->Changed = 1;
+ v->Links = 0;
+ v->AliasOf = 0;
+ v->ScalarOnly = 1;
+ v->TextOnly = 1;
+ v->normals = new smooth_normals(v->AngleSmoothNormals);
+ v->Min = VAL_INF;
+ v->Max = -VAL_INF;
+ v->adaptive = 0;
+ for(i = 0; i < 3; i++) {
+ v->BBox[2 * i] = v->TmpBBox[2 * i] = VAL_INF;
+ v->BBox[2 * i + 1] = v->TmpBBox[2 * i + 1] = -VAL_INF;
+ }
+ for(i = 0; i < 3; i++)
+ v->GenRaise_f[i] = NULL;
+
+ return v;
+}
+
+void Post_View::reset_normals()
+{
+ if(normals)
+ delete normals;
+ normals = new smooth_normals(AngleSmoothNormals);
+}
+
+double ComputeVonMises(double *V)
+{
+ static const double THIRD = 1.e0 / 3.e0;
+ double tr = (V[0] + V[4] + V[8]) * THIRD;
+ double v11 = V[0] - tr;
+ double v12 = V[1];
+ double v13 = V[2];
+ double v21 = V[3];
+ double v22 = V[4] - tr;
+ double v23 = V[5];
+ double v31 = V[6];
+ double v32 = V[7];
+ double v33 = V[8] - tr;
+ return sqrt(1.5 * (v11 * v11 + v12 * v12 + v13 * v13 +
+ v21 * v21 + v22 * v22 + v23 * v23 +
+ v31 * v31 + v32 * v32 + v33 * v33));
+}
+
+void Stat_Element(Post_View *v, int type, int nbnod, int N,
+ double *X, double *Y, double *Z, double *V)
+{
+ int i;
+ double l0;
+
+ switch (type) {
+
+ case 0: // scalar
+ if(v->Min == VAL_INF || v->Max == -VAL_INF) {
+ v->NbTimeStep = N / nbnod;
+ if(v->TimeStepMin) Free(v->TimeStepMin);
+ if(v->TimeStepMax) Free(v->TimeStepMax);
+ v->TimeStepMin = (double*)Malloc(v->NbTimeStep * sizeof(double));
+ v->TimeStepMax = (double*)Malloc(v->NbTimeStep * sizeof(double));
+ for(i = 0; i < v->NbTimeStep; i++){
+ v->TimeStepMin[i] = VAL_INF;
+ v->TimeStepMax[i] = -VAL_INF;
+ }
+ }
+ else if(N / nbnod < v->NbTimeStep){
+ // if some elts have less steps, reduce the total number!
+ v->NbTimeStep = N / nbnod;
+ }
+
+ for(i = 0; i < N; i++) {
+ l0 = V[i];
+ if(l0 < v->Min)
+ v->Min = l0;
+ if(l0 > v->Max)
+ v->Max = l0;
+ int ts = i / nbnod;
+ if(ts < v->NbTimeStep){ // security
+ if(l0 < v->TimeStepMin[ts])
+ v->TimeStepMin[ts] = l0;
+ if(l0 > v->TimeStepMax[ts])
+ v->TimeStepMax[ts] = l0;
+ }
+ }
+ break;
+
+ case 1: // vector
+ if(v->Min == VAL_INF || v->Max == -VAL_INF) {
+ v->NbTimeStep = N / (3 * nbnod);
+ if(v->TimeStepMin) Free(v->TimeStepMin);
+ if(v->TimeStepMax) Free(v->TimeStepMax);
+ v->TimeStepMin = (double*)Malloc(v->NbTimeStep * sizeof(double));
+ v->TimeStepMax = (double*)Malloc(v->NbTimeStep * sizeof(double));
+ for(i = 0; i < v->NbTimeStep; i++){
+ v->TimeStepMin[i] = VAL_INF;
+ v->TimeStepMax[i] = -VAL_INF;
+ }
+ }
+ else if(N / (3 * nbnod) < v->NbTimeStep){
+ v->NbTimeStep = N / (3 * nbnod);
+ }
+
+ for(i = 0; i < N; i += 3) {
+ l0 = sqrt(DSQR(V[i]) + DSQR(V[i + 1]) + DSQR(V[i + 2]));
+ if(l0 < v->Min)
+ v->Min = l0;
+ if(l0 > v->Max)
+ v->Max = l0;
+ int ts = i / (3 * nbnod);
+ if(ts < v->NbTimeStep){ // security
+ if(l0 < v->TimeStepMin[ts])
+ v->TimeStepMin[ts] = l0;
+ if(l0 > v->TimeStepMax[ts])
+ v->TimeStepMax[ts] = l0;
+ }
+ }
+ v->ScalarOnly = 0;
+ break;
+
+ case 2: // tensor
+ if(v->Min == VAL_INF || v->Max == -VAL_INF) {
+ v->NbTimeStep = N / (9 * nbnod);
+ if(v->TimeStepMin) Free(v->TimeStepMin);
+ if(v->TimeStepMax) Free(v->TimeStepMax);
+ v->TimeStepMin = (double*)Malloc(v->NbTimeStep * sizeof(double));
+ v->TimeStepMax = (double*)Malloc(v->NbTimeStep * sizeof(double));
+ for(i = 0; i < v->NbTimeStep; i++){
+ v->TimeStepMin[i] = VAL_INF;
+ v->TimeStepMax[i] = -VAL_INF;
+ }
+ }
+ else if(N / (9 * nbnod) < v->NbTimeStep){
+ v->NbTimeStep = N / (9 * nbnod);
+ }
+
+ for(i = 0; i < N; i += 9) {
+ // by lack of any current better solution, tensors are displayed
+ // as their Von Mises invariant (J2 invariant)
+ l0 = ComputeVonMises(V+i);
+ if(l0 < v->Min)
+ v->Min = l0;
+ if(l0 > v->Max)
+ v->Max = l0;
+ int ts = i / (9 * nbnod);
+ if(ts < v->NbTimeStep){ // security
+ if(l0 < v->TimeStepMin[ts])
+ v->TimeStepMin[ts] = l0;
+ if(l0 > v->TimeStepMax[ts])
+ v->TimeStepMax[ts] = l0;
+ }
+ }
+ v->ScalarOnly = 0;
+ break;
+
+ }
+
+ for(i = 0; i < nbnod; i++) {
+ if(X[i] < v->BBox[0]) v->BBox[0] = X[i];
+ if(X[i] > v->BBox[1]) v->BBox[1] = X[i];
+ if(Y[i] < v->BBox[2]) v->BBox[2] = Y[i];
+ if(Y[i] > v->BBox[3]) v->BBox[3] = Y[i];
+ if(Z[i] < v->BBox[4]) v->BBox[4] = Z[i];
+ if(Z[i] > v->BBox[5]) v->BBox[5] = Z[i];
+ }
+
+ v->TextOnly = 0;
+}
+
+void Stat_List(Post_View * v, List_T * listelm, int type, int nbelm,
+ int nbnod)
+{
+ int i, nb;
+ if(nbelm) {
+ nb = List_Nbr(listelm) / nbelm;
+ for(i = 0; i < List_Nbr(listelm); i += nb)
+ Stat_Element(v, type, nbnod, nb - 3 * nbnod,
+ (double *)List_Pointer_Fast(listelm, i),
+ (double *)List_Pointer_Fast(listelm, i + 1 * nbnod),
+ (double *)List_Pointer_Fast(listelm, i + 2 * nbnod),
+ (double *)List_Pointer_Fast(listelm, i + 3 * nbnod));
+ }
+}
+
+void Stat_Text(Post_View * v, List_T *D, List_T *C, int nb)
+{
+ for(int i = 0; i < List_Nbr(D); i += nb){
+ double beg, end;
+ List_Read(D, i+nb-1, &beg);
+ if(i > List_Nbr(D) - 2*nb)
+ end = (double)List_Nbr(C);
+ else
+ List_Read(D, i+nb+nb-1, &end);
+ char *c = (char*)List_Pointer(C, (int)beg);
+ int nbtime = 0;
+ for(int j = 0; j < (int)(end-beg); j++)
+ if(c[j] == '\0') nbtime++;
+ if(nbtime > v->NbTimeStep)
+ v->NbTimeStep = nbtime;
+ }
+
+ if(nb == 5){
+ for(int i = 0; i < List_Nbr(D); i += nb){
+ double x, y, z;
+ List_Read(D, i, &x);
+ List_Read(D, i+1, &y);
+ List_Read(D, i+2, &z);
+ if(x < v->BBox[0]) v->BBox[0] = x;
+ if(x > v->BBox[1]) v->BBox[1] = x;
+ if(y < v->BBox[2]) v->BBox[2] = y;
+ if(y > v->BBox[3]) v->BBox[3] = y;
+ if(z < v->BBox[4]) v->BBox[4] = z;
+ if(z > v->BBox[5]) v->BBox[5] = z;
+ }
+ }
+}
+
+void EndView(Post_View * v, int add_in_gui, char *file_name, char *name)
+{
+ int i;
+ double d;
+
+ // Stat text strings first, to get the max value of NbTimeStep for
+ // srings-only (strings are designed to degrade gracefully when some
+ // have fewer time steps than others). If there are any elements in
+ // the view, this value will be replaced by the minimum number of
+ // time steps common to all elements.
+ Stat_Text(v, v->T2D, v->T2C, 4);
+ Stat_Text(v, v->T3D, v->T3C, 5);
+
+ // convert all curved (geometrically 2nd order) elements into linear
+ // elements *AND* free all the data associated with the curved
+ // elements
+ v->splitCurvedElements();
+
+ // Points
+ Stat_List(v, v->SP, 0, v->NbSP, 1);
+ Stat_List(v, v->VP, 1, v->NbVP, 1);
+ Stat_List(v, v->TP, 2, v->NbTP, 1);
+
+ // Lines
+ Stat_List(v, v->SL, 0, v->NbSL, 2);
+ Stat_List(v, v->VL, 1, v->NbVL, 2);
+ Stat_List(v, v->TL, 2, v->NbTL, 2);
+
+ // Triangles
+ Stat_List(v, v->ST, 0, v->NbST, 3);
+ Stat_List(v, v->VT, 1, v->NbVT, 3);
+ Stat_List(v, v->TT, 2, v->NbTT, 3);
+
+ // Quadrangles
+ Stat_List(v, v->SQ, 0, v->NbSQ, 4);
+ Stat_List(v, v->VQ, 1, v->NbVQ, 4);
+ Stat_List(v, v->TQ, 2, v->NbTQ, 4);
+
+ // Tetrahedra
+ Stat_List(v, v->SS, 0, v->NbSS, 4);
+ Stat_List(v, v->VS, 1, v->NbVS, 4);
+ Stat_List(v, v->TS, 2, v->NbTS, 4);
+
+ // Hexahedra
+ Stat_List(v, v->SH, 0, v->NbSH, 8);
+ Stat_List(v, v->VH, 1, v->NbVH, 8);
+ Stat_List(v, v->TH, 2, v->NbTH, 8);
+
+ // Prisms
+ Stat_List(v, v->SI, 0, v->NbSI, 6);
+ Stat_List(v, v->VI, 1, v->NbVI, 6);
+ Stat_List(v, v->TI, 2, v->NbTI, 6);
+
+ // Pyramids
+ Stat_List(v, v->SY, 0, v->NbSY, 5);
+ Stat_List(v, v->VY, 1, v->NbVY, 5);
+ Stat_List(v, v->TY, 2, v->NbTY, 5);
+
+ // Dummy time values if none (or too few) provided (e.g. using old
+ // parsed format)
+ if(v->Time && List_Nbr(v->Time) < v->NbTimeStep) {
+ for(i = List_Nbr(v->Time); i < v->NbTimeStep; i++) {
+ d = (double)i;
+ List_Add(v->Time, &d);
+ }
+ }
+
+ opt_view_name(v->Index, GMSH_SET | GMSH_GUI, name);
+ opt_view_filename(v->Index, GMSH_SET | GMSH_GUI, file_name);
+ opt_view_nb_timestep(v->Index, GMSH_GUI, 0);
+ opt_view_timestep(v->Index, GMSH_SET | GMSH_GUI, v->TimeStep);
+
+ if(CTX.post.smooth)
+ v->smooth();
+
+#if defined(HAVE_FLTK)
+ if(add_in_gui)
+ UpdateViewsInGUI();
+#endif
+
+ v->Dirty = 0; //the view is complete, we may draw it
+
+ Msg(DEBUG, "Added View[%d]", v->Index);
+}
+
+void AliasView(int index, int withoptions)
+{
+ if(index < 0 || index >= List_Nbr(CTX.post.list)) {
+ return;
+ }
+
+ Post_View v, *pv, **ppv;
+
+ Post_View *v1 = *(Post_View **) List_Pointer(CTX.post.list, index);
+
+ Post_View *v2 = BeginView(0);
+ EndView(v2, 0, v1->FileName, v1->Name);
+
+ if(!v1->AliasOf) {
+ v2->AliasOf = v1->Num;
+ v1->Links++;
+ }
+ else {
+ v.Num = v1->AliasOf;
+ pv = &v;
+ if(!(ppv = (Post_View **) List_PQuery(CTX.post.list, &pv, fcmpPostViewNum))) {
+ v2->AliasOf = v1->Num;
+ v1->Links++;
+ }
+ else {
+ v2->AliasOf = (*ppv)->Num;
+ (*ppv)->Links++;
+ }
+ }
+
+ // When we create an alias we just point to a reference view: we
+ // DON'T allocate a new data set!
+ v2->Time = v1->Time;
+ v2->TimeStepMin = v1->TimeStepMin;
+ v2->TimeStepMax = v1->TimeStepMax;
+
+ v2->NbSP = v1->NbSP; v2->SP = v1->SP;
+ v2->NbVP = v1->NbVP; v2->VP = v1->VP;
+ v2->NbTP = v1->NbTP; v2->TP = v1->TP;
+
+ v2->NbSL = v1->NbSL; v2->SL = v1->SL;
+ v2->NbVL = v1->NbVL; v2->VL = v1->VL;
+ v2->NbTL = v1->NbTL; v2->TL = v1->TL;
+
+ v2->NbST = v1->NbST; v2->ST = v1->ST;
+ v2->NbVT = v1->NbVT; v2->VT = v1->VT;
+ v2->NbTT = v1->NbTT; v2->TT = v1->TT;
+
+ v2->NbSQ = v1->NbSQ; v2->SQ = v1->SQ;
+ v2->NbVQ = v1->NbVQ; v2->VQ = v1->VQ;
+ v2->NbTQ = v1->NbTQ; v2->TQ = v1->TQ;
+
+ v2->NbSS = v1->NbSS; v2->SS = v1->SS;
+ v2->NbVS = v1->NbVS; v2->VS = v1->VS;
+ v2->NbTS = v1->NbTS; v2->TS = v1->TS;
+
+ v2->NbSH = v1->NbSH; v2->SH = v1->SH;
+ v2->NbVH = v1->NbVH; v2->VH = v1->VH;
+ v2->NbTH = v1->NbTH; v2->TH = v1->TH;
+
+ v2->NbSI = v1->NbSI; v2->SI = v1->SI;
+ v2->NbVI = v1->NbVI; v2->VI = v1->VI;
+ v2->NbTI = v1->NbTI; v2->TI = v1->TI;
+
+ v2->NbSY = v1->NbSY; v2->SY = v1->SY;
+ v2->NbVY = v1->NbVY; v2->VY = v1->VY;
+ v2->NbTY = v1->NbTY; v2->TY = v1->TY;
+
+ v2->NbT2 = v1->NbT2; v2->T2D = v1->T2D; v2->T2C = v1->T2C;
+ v2->NbT3 = v1->NbT3; v2->T3D = v1->T3D; v2->T3C = v1->T3C;
+
+ v2->DataSize = v1->DataSize;
+ v2->ScalarOnly = v1->ScalarOnly;
+ v2->TextOnly = v1->TextOnly;
+ v2->Min = v1->Min;
+ v2->Max = v1->Max;
+ v2->NbTimeStep = v1->NbTimeStep;
+ for(int i=0 ; i<6 ; i++)
+ v2->BBox[i] = v1->BBox[i];
+
+ if(withoptions)
+ CopyViewOptions(v1, v2);
+
+#if defined(HAVE_FLTK)
+ UpdateViewsInGUI();
+#endif
+}
+
+bool RemoveViewByIndex(int index)
+{
+ if(index < 0 || index >= List_Nbr(CTX.post.list)) {
+ return false;
+ }
+
+ Post_View *v = *(Post_View **) List_Pointer(CTX.post.list, index);
+ FreeView(v);
+ List_PSuppress(CTX.post.list, index);
+
+ // recalculate the indices
+ for(int i = 0; i < List_Nbr(CTX.post.list); i++){
+ v = *(Post_View **) List_Pointer(CTX.post.list, i);
+ v->Index = i;
+ }
+
+ Msg(DEBUG, "Removed View[%d] (%d views left)", index, List_Nbr(CTX.post.list));
+ return true;
+}
+
+bool RemoveViewByNumber(int num)
+{
+ Post_View v, *pv;
+
+ v.Num = num;
+ pv = &v;
+ int i = List_ISearch(CTX.post.list, &pv, fcmpPostViewNum);
+
+ return RemoveViewByIndex(i);
+}
+
+void FreeView(Post_View * v)
+{
+ Post_View vv, *pvv, **ppvv;
+ int i, numdup, free = 1;
+
+ if(v->AliasOf) {
+ vv.Num = v->AliasOf;
+ pvv = &vv;
+ Msg(DEBUG, "This view is a duplicata");
+ if(!(ppvv = (Post_View **) List_PQuery(CTX.post.list, &pvv, fcmpPostViewNum))) {
+ Msg(DEBUG, " -the original view is gone");
+ numdup = 0;
+ for(i = 0; i < List_Nbr(CTX.post.list); i++)
+ numdup += ((*(Post_View **) List_Pointer(CTX.post.list, i))->AliasOf
+ == v->AliasOf);
+ if(numdup == 1) {
+ Msg(DEBUG, " -there are no other duplicata, so I can free");
+ free = 1;
+ }
+ else {
+ Msg(DEBUG, " -there are still duplicata, so I cannot free");
+ free = 0;
+ }
+ }
+ else {
+ (*ppvv)->Links--;
+ free = 0;
+ Msg(DEBUG, " -the original still exists, so I cannot free");
+ }
+ }
+
+ if(free && !v->Links) {
+ Msg(DEBUG, "Freeing View!");
+ List_Delete(v->Time);
+ Free(v->TimeStepMin);
+ Free(v->TimeStepMax);
+ List_Delete(v->SP); List_Delete(v->VP); List_Delete(v->TP);
+ List_Delete(v->SL); List_Delete(v->VL); List_Delete(v->TL);
+ List_Delete(v->ST); List_Delete(v->VT); List_Delete(v->TT);
+ List_Delete(v->SQ); List_Delete(v->VQ); List_Delete(v->TQ);
+ List_Delete(v->SS); List_Delete(v->VS); List_Delete(v->TS);
+ List_Delete(v->SH); List_Delete(v->VH); List_Delete(v->TH);
+ List_Delete(v->SI); List_Delete(v->VI); List_Delete(v->TI);
+ List_Delete(v->SY); List_Delete(v->VY); List_Delete(v->TY);
+ List_Delete(v->T2D); List_Delete(v->T2C);
+ List_Delete(v->T3D); List_Delete(v->T3C);
+ delete v->Grains;
+ delete v->DisplayListsOfGrains;
+ // Note: all the second order elements have already been freed in xxxx
+ if(v->normals) delete v->normals;
+ if(v->TriVertexArray) delete v->TriVertexArray;
+ if(v->LinVertexArray) delete v->LinVertexArray;
+ if(v->adaptive) delete v->adaptive;
+ FreeGeneralizedRaise(v);
+ Free(v);
+ }
+}
+
+void CopyViewOptions(Post_View * src, Post_View * dest)
+{
+ dest->Type = src->Type;
+ dest->AutoPosition = src->AutoPosition;
+ dest->AxesAutoPosition = src->AxesAutoPosition;
+ dest->Position[0] = src->Position[0];
+ dest->Position[1] = src->Position[1];
+ dest->Size[0] = src->Size[0];
+ dest->Size[1] = src->Size[1];
+ for(int i = 0; i < 6; i++) dest->AxesPosition[i] = src->AxesPosition[i];
+ for(int i = 0; i < 3; i++) dest->AxesTics[i] = src->AxesTics[i];
+ strcpy(dest->Format, src->Format);
+ for(int i = 0; i < 3; i++) strcpy(dest->AxesFormat[i], src->AxesFormat[i]);
+ for(int i = 0; i < 3; i++) strcpy(dest->AxesLabel[i], src->AxesLabel[i]);
+ dest->CustomMin = src->CustomMin;
+ dest->CustomMax = src->CustomMax;
+ dest->Offset[0] = src->Offset[0];
+ dest->Offset[1] = src->Offset[1];
+ dest->Offset[2] = src->Offset[2];
+ for(int i = 0; i < 3; i++)
+ for(int j = 0; j < 3; j++)
+ dest->Transform[i][j] = src->Transform[i][j];
+ dest->Raise[0] = src->Raise[0];
+ dest->Raise[1] = src->Raise[1];
+ dest->Raise[2] = src->Raise[2];
+ dest->ArrowSize = src->ArrowSize;
+ dest->ArrowSizeProportional = src->ArrowSizeProportional;
+ dest->ArrowRelHeadRadius = src->ArrowRelHeadRadius;
+ dest->ArrowRelStemLength = src->ArrowRelStemLength;
+ dest->ArrowRelStemRadius = src->ArrowRelStemRadius;
+ dest->Normals = src->Normals;
+ dest->Tangents = src->Tangents;
+ dest->DisplacementFactor = src->DisplacementFactor;
+ dest->Explode = src->Explode;
+ dest->Visible = src->Visible;
+ dest->IntervalsType = src->IntervalsType;
+ dest->SaturateValues = src->SaturateValues;
+ dest->Boundary = src->Boundary;
+ dest->NbIso = src->NbIso;
+ dest->Light = src->Light;
+ dest->LightTwoSide = src->LightTwoSide;
+ dest->LightLines = src->LightLines;
+ dest->SmoothNormals = src->SmoothNormals;
+ dest->AngleSmoothNormals = src->AngleSmoothNormals;
+ dest->ShowElement = src->ShowElement;
+ dest->ShowTime = src->ShowTime;
+ dest->ShowScale = src->ShowScale;
+ dest->DrawPoints = src->DrawPoints;
+ dest->DrawLines = src->DrawLines;
+ dest->DrawTriangles = src->DrawTriangles;
+ dest->DrawQuadrangles = src->DrawQuadrangles;
+ dest->DrawTetrahedra = src->DrawTetrahedra;
+ dest->DrawHexahedra = src->DrawHexahedra;
+ dest->DrawPrisms = src->DrawPrisms;
+ dest->DrawPyramids = src->DrawPyramids;
+ dest->DrawScalars = src->DrawScalars;
+ dest->DrawVectors = src->DrawVectors;
+ dest->DrawTensors = src->DrawTensors;
+ dest->DrawStrings = src->DrawStrings;
+ dest->ScaleType = src->ScaleType;
+ dest->RangeType = src->RangeType;
+ dest->VectorType = src->VectorType;
+ dest->GlyphLocation = src->GlyphLocation;
+ dest->TensorType = src->TensorType;
+ dest->TimeStep = src->TimeStep;
+ dest->PointSize = src->PointSize;
+ dest->LineWidth = src->LineWidth;
+ dest->PointType = src->PointType;
+ dest->LineType = src->LineType;
+ dest->Axes = src->Axes;
+ dest->ExternalViewIndex = src->ExternalViewIndex;
+ dest->ViewIndexForGenRaise = src->ViewIndexForGenRaise;
+ dest->UseGenRaise = src->UseGenRaise;
+ dest->FakeTransparency = src->FakeTransparency;
+ dest->GenRaiseFactor = src->GenRaiseFactor;
+ strcpy(dest->GenRaiseX, src->GenRaiseX);
+ strcpy(dest->GenRaiseY, src->GenRaiseY);
+ strcpy(dest->GenRaiseZ, src->GenRaiseZ);
+ dest->UseStipple = src->UseStipple;
+ for(int i = 0; i < 10; i++){
+ dest->Stipple[i][0] = src->Stipple[i][0];
+ dest->Stipple[i][1] = src->Stipple[i][1];
+ strcpy(dest->StippleString[i], src->StippleString[i]);
+ }
+ dest->color.point = src->color.point;
+ dest->color.line = src->color.line;
+ dest->color.triangle = src->color.triangle;
+ dest->color.quadrangle = src->color.quadrangle;
+ dest->color.tetrahedron = src->color.tetrahedron;
+ dest->color.hexahedron = src->color.hexahedron;
+ dest->color.prism = src->color.prism;
+ dest->color.pyramid = src->color.pyramid;
+ dest->color.tangents = src->color.tangents;
+ dest->color.normals = src->color.normals;
+ dest->color.text2d = src->color.text2d;
+ dest->color.text3d = src->color.text3d;
+ dest->color.axes = src->color.axes;
+ ColorTable_Copy(&src->CT);
+ ColorTable_Paste(&dest->CT);
+}
+
+GmshColorTable *Get_ColorTable(int num)
+{
+ Post_View *v;
+
+ if(!CTX.post.list){
+ v = Post_ViewReference;
+ if(!v) return NULL;
+ }
+ else{
+ Post_View **vv = (Post_View **) List_Pointer_Test(CTX.post.list, num);
+ if(!vv) return NULL;
+ v = *vv;
+ }
+ v->Changed = 1; // let's assume that if we get the ct we change it
+ return &v->CT;
+}
+
+void Print_ColorTable(int num, int diff, char *prefix, FILE * file)
+{
+ char tmp[1024];
+ Post_View *v;
+
+ if(!CTX.post.list){
+ v = Post_ViewReference;
+ if(!v) return;
+ }
+ else{
+ Post_View **vv = (Post_View **) List_Pointer_Test(CTX.post.list, num);
+ if(!vv) return;
+ v = *vv;
+ }
+
+ if(diff){
+ // compare the current colormap with a vanilla colormap having the
+ // parameters
+ GmshColorTable ref;
+ ColorTable_InitParam(v->CT.ipar[COLORTABLE_NUMBER], &ref);
+ for(int i = 0; i < COLORTABLE_NBMAX_PARAM; i++){
+ ref.ipar[i] = v->CT.ipar[i];
+ ref.dpar[i] = v->CT.dpar[i];
+ }
+ ColorTable_Recompute(&ref);
+ if(!ColorTable_Diff(&ref, &v->CT))
+ return;
+ }
+
+ sprintf(tmp, "%s = {", prefix);
+ if(file)
+ fprintf(file, "%s\n", tmp);
+ else
+ Msg(DIRECT, tmp);
+ ColorTable_Print(&v->CT, file);
+ sprintf(tmp, "};");
+ if(file)
+ fprintf(file, "%s\n", tmp);
+ else
+ Msg(DIRECT, tmp);
+}
+
+Post_View *Create2DGraph(char *xname, char *yname,
+ int nbdata, double *x, double *y)
+{
+ Post_View *v = BeginView(1);
+ for(int i = 0; i < nbdata; i++) {
+ double d;
+ if(x){
+ List_Add(v->SP, &x[i]);
+ }
+ else{
+ d = nbdata > 1 ? (double)i/(double)(nbdata - 1) : 0;
+ List_Add(v->SP, &d);
+ }
+ d = 0.;
+ List_Add(v->SP, &d);
+ List_Add(v->SP, &d);
+ List_Add(v->SP, &y[i]);
+ v->NbSP++;
+ }
+ char filename[1024];
+ sprintf(filename, "%s.pos", yname);
+ EndView(v, 1, filename, yname);
+ v->Type = DRAW_POST_2D_SPACE;
+ v->Axes = 2;
+ strcpy(v->AxesLabel[0], xname);
+ return v;
+}
+
+// Transform curved elements into linear ones and delete all the
+// curved element data. This is a temporary solution, until we can use
+// an Adaptive_Post_View on curved elements, too.
+
+void splitCurvedElement(List_T **in, int *nbin, List_T *out, int *nbout,
+ int nodin, int nodout, int nbcomp, int nbsplit, int split[][8],
+ int remove=1)
+{
+ if(*nbin){
+ int nb = List_Nbr(*in) / *nbin;
+ int nbstep = (nb - 3 * nodin) / (nodin * nbcomp); // we don't know this yet for the view
+ for(int i = 0; i < List_Nbr(*in); i += nb) {
+ double *coord = (double *)List_Pointer_Fast(*in, i);
+ double *val = (double *)List_Pointer_Fast(*in, i + 3 * nodin);
+ for(int j = 0; j < nbsplit; j++){
+ for(int k = 0; k < nodout; k++)
+ List_Add(out, &coord[split[j][k]]);
+ for(int k = 0; k < nodout; k++)
+ List_Add(out, &coord[nodin + split[j][k]]);
+ for(int k = 0; k < nodout; k++)
+ List_Add(out, &coord[2 * nodin + split[j][k]]);
+ for(int ts = 0; ts < nbstep; ts++){
+ for(int k = 0; k < nodout; k++){
+ for(int l = 0; l < nbcomp; l++){
+ List_Add(out, &val[nodin * nbcomp * ts + nbcomp * split[j][k] + l]);
+ }
+ }
+ }
+ (*nbout)++;
+ }
+ }
+ }
+
+ if(remove){
+ *nbin = 0;
+ List_Delete(*in);
+ *in = NULL;
+ }
+}
+
+void Post_View::splitCurvedElements()
+{
+ // we could keep track of the starting index in SL, VL, ..., so that
+ // we could draw the boundaries correctly
+
+ int lin[2][8] = { // 2-split
+ {0,2}, {2,1}
+ };
+ splitCurvedElement(&SL2, &NbSL2, SL, &NbSL, 3,2, 1, 2, lin);
+ splitCurvedElement(&VL2, &NbVL2, VL, &NbVL, 3,2, 3, 2, lin);
+ splitCurvedElement(&TL2, &NbTL2, TL, &NbTL, 3,2, 9, 2, lin);
+
+ int tri[4][8] = { // 4-split
+ {0,3,5}, {1,4,3}, {2,5,4}, {3,4,5}
+ };
+ splitCurvedElement(&ST2, &NbST2, ST, &NbST, 6,3, 1, 4, tri);
+ splitCurvedElement(&VT2, &NbVT2, VT, &NbVT, 6,3, 3, 4, tri);
+ splitCurvedElement(&TT2, &NbTT2, TT, &NbTT, 6,3, 9, 4, tri);
+
+ int qua[4][8] = { // 4-split
+ {0,4,8,7}, {1,5,8,4}, {2,6,8,5}, {3,7,8,6}
+ };
+ splitCurvedElement(&SQ2, &NbSQ2, SQ, &NbSQ, 9,4, 1, 4, qua);
+ splitCurvedElement(&VQ2, &NbVQ2, VQ, &NbVQ, 9,4, 3, 4, qua);
+ splitCurvedElement(&TQ2, &NbTQ2, TQ, &NbTQ, 9,4, 9, 4, qua);
+
+ int tet[8][8] = { // 8-split
+ {0,4,6,7}, {1,5,4,9}, {2,6,5,8}, {3,9,7,8},
+ {4,6,7,8}, {4,6,5,8}, {4,5,9,8}, {4,7,9,8}
+ };
+ splitCurvedElement(&SS2, &NbSS2, SS, &NbSS, 10,4, 1, 8, tet);
+ splitCurvedElement(&VS2, &NbVS2, VS, &NbVS, 10,4, 3, 8, tet);
+ splitCurvedElement(&TS2, &NbTS2, TS, &NbTS, 10,4, 9, 8, tet);
+
+ int hex[8][8] = { // 8-split
+ {0,8,20,9, 10,21,26,22}, {8,1,11,20, 21,12,23,26},
+ {9,20,13,3, 22,26,24,15}, {20,11,2,13, 26,23,14,24},
+ {10,21,26,22, 4,16,25,17}, {21,12,23,26, 16,5,18,25},
+ {22,26,24,15, 17,25,19,7}, {26,23,14,24, 25,18,6,19}
+ };
+ splitCurvedElement(&SH2, &NbSH2, SH, &NbSH, 27,8, 1, 8, hex);
+ splitCurvedElement(&VH2, &NbVH2, VH, &NbVH, 27,8, 3, 8, hex);
+ splitCurvedElement(&TH2, &NbTH2, TH, &NbTH, 27,8, 9, 8, hex);
+
+ int pri[8][8] = { // 8-split
+ {0,6,7, 8,15,16}, {1,9,6, 10,17,15}, {2,7,9, 11,16,17}, {6,9,7, 15,17,16},
+ {8,15,16, 3,12,13}, {10,17,15, 4,14,12}, {11,16,17, 5,13,14}, {15,17,16, 12,14,13}
+ };
+ splitCurvedElement(&SI2, &NbSI2, SI, &NbSI, 18,6, 1, 8, pri);
+ splitCurvedElement(&VI2, &NbVI2, VI, &NbVI, 18,6, 3, 8, pri);
+ splitCurvedElement(&TI2, &NbTI2, TI, &NbTI, 18,6, 9, 8, pri);
+
+ int pyr[6][8] = { // 6 pyramids
+ {0,5,13,6, 7}, {5,1,8,13, 9}, {6,13,10,3, 12}, {13,8,2,10, 11},
+ {7,9,11,12, 4}, {7,12,11,9, 13}
+ };
+ splitCurvedElement(&SY2, &NbSY2, SY, &NbSY, 14,5, 1, 6, pyr, 0); // don't remove yet
+ splitCurvedElement(&VY2, &NbVY2, VY, &NbVY, 14,5, 3, 6, pyr, 0);
+ splitCurvedElement(&TY2, &NbTY2, TY, &NbTY, 14,5, 9, 6, pyr, 0);
+
+ int pyr2[4][8] = { // + 4 tets to fill the holes
+ {6,12,7,13}, {7,9,5,13}, {9,11,8,13}, {12,10,11,13}
+ };
+ splitCurvedElement(&SY2, &NbSY2, SS, &NbSS, 14,4, 1, 4, pyr2);
+ splitCurvedElement(&VY2, &NbVY2, VS, &NbVS, 14,4, 3, 4, pyr2);
+ splitCurvedElement(&TY2, &NbTY2, TS, &NbTS, 14,4, 9, 4, pyr2);
+}
+
+// Smoothing
+
+void generate_connectivities(List_T *list, int nbList,
+ int nbTimeStep, int nbVert, smooth_data &data)
+{
+ if(!nbList) return;
+
+ double *vals = new double[nbTimeStep];
+ int nb = List_Nbr(list) / nbList;
+ for(int i = 0; i < List_Nbr(list); i += nb) {
+ double *x = (double *)List_Pointer_Fast(list, i);
+ double *y = (double *)List_Pointer_Fast(list, i + nbVert);
+ double *z = (double *)List_Pointer_Fast(list, i + 2 * nbVert);
+ double *v = (double *)List_Pointer_Fast(list, i + 3 * nbVert);
+ for(int j = 0; j < nbVert; j++) {
+ for(int k = 0; k < nbTimeStep; k++)
+ vals[k] = v[j + k * nbVert];
+ data.add(x[j], y[j], z[j], nbTimeStep, vals);
+ }
+ }
+ delete [] vals;
+}
+
+void smooth_list(List_T *list, int nbList,
+ double *min, double *max, double *tsmin, double *tsmax,
+ int nbTimeStep, int nbVert, smooth_data &data)
+{
+ if(!nbList)
+ return;
+
+ double *vals = new double[nbTimeStep];
+ int nb = List_Nbr(list)/nbList;
+ for(int i = 0; i < List_Nbr(list); i += nb) {
+ double *x = (double *)List_Pointer_Fast(list, i);
+ double *y = (double *)List_Pointer_Fast(list, i + nbVert);
+ double *z = (double *)List_Pointer_Fast(list, i + 2 * nbVert);
+ double *v = (double *)List_Pointer_Fast(list, i + 3 * nbVert);
+ for(int j = 0; j < nbVert; j++) {
+ if(data.get(x[j], y[j], z[j], nbTimeStep, vals)){
+ for(int k = 0; k < nbTimeStep; k++) {
+ double dd = vals[k];
+ v[j + k * nbVert] = dd;
+ if(dd < *min) *min = dd;
+ if(dd > *max) *max = dd;
+ if(dd < tsmin[k]) tsmin[k] = dd;
+ if(dd > tsmax[k]) tsmax[k] = dd;
+ }
+ }
+ }
+ }
+ delete [] vals;
+}
+
+void Post_View::smooth()
+{
+ double old_eps = xyzv::eps;
+ xyzv::eps = CTX.lc * 1.e-8;
+
+ if(NbSL || NbST || NbSQ || NbSS || NbSH || NbSI || NbSY) {
+ Msg(INFO, "Smoothing scalar primitives in View[%d]", Index);
+ Min = VAL_INF;
+ Max = -VAL_INF;
+ for(int k = 0; k < NbTimeStep; k++) {
+ TimeStepMin[k] = VAL_INF;
+ TimeStepMax[k] = -VAL_INF;
+ }
+ smooth_data data;
+ generate_connectivities(SL, NbSL, NbTimeStep, 2, data);
+ generate_connectivities(ST, NbST, NbTimeStep, 3, data);
+ generate_connectivities(SQ, NbSQ, NbTimeStep, 4, data);
+ generate_connectivities(SS, NbSS, NbTimeStep, 4, data);
+ generate_connectivities(SH, NbSH, NbTimeStep, 8, data);
+ generate_connectivities(SI, NbSI, NbTimeStep, 6, data);
+ generate_connectivities(SY, NbSY, NbTimeStep, 5, data);
+ smooth_list(SL, NbSL, &Min, &Max, TimeStepMin, TimeStepMax, NbTimeStep, 2, data);
+ smooth_list(ST, NbST, &Min, &Max, TimeStepMin, TimeStepMax, NbTimeStep, 3, data);
+ smooth_list(SQ, NbSQ, &Min, &Max, TimeStepMin, TimeStepMax, NbTimeStep, 4, data);
+ smooth_list(SS, NbSS, &Min, &Max, TimeStepMin, TimeStepMax, NbTimeStep, 4, data);
+ smooth_list(SH, NbSH, &Min, &Max, TimeStepMin, TimeStepMax, NbTimeStep, 8, data);
+ smooth_list(SI, NbSI, &Min, &Max, TimeStepMin, TimeStepMax, NbTimeStep, 6, data);
+ smooth_list(SY, NbSY, &Min, &Max, TimeStepMin, TimeStepMax, NbTimeStep, 5, data);
+ Changed = 1;
+ }
+
+ xyzv::eps = old_eps;
+}
+
+// Combine views (merge elements or merge time steps)
+
+struct nameidx{
+ char name[256];
+ List_T *indices;
+};
+
+static int fcmp_name(const void *a, const void *b){
+ char *name1 = ((struct nameidx*)a)->name;
+ char *name2 = ((struct nameidx*)b)->name;
+ return strcmp(name1, name2);
+}
+
+static void combine(List_T * a, List_T * b)
+{
+ if(!a || !b)
+ return;
+ for(int i = 0; i < List_Nbr(a); i++) {
+ List_Add(b, List_Pointer(a, i));
+ }
+}
+
+static void combine_strings(Post_View *a, Post_View *b)
+{
+ for(int i = 0; i < List_Nbr(a->T2D); i+=4){
+ List_Add(b->T2D, List_Pointer(a->T2D, i));
+ List_Add(b->T2D, List_Pointer(a->T2D, i+1));
+ List_Add(b->T2D, List_Pointer(a->T2D, i+2));
+ double d = List_Nbr(b->T2C);
+ List_Add(b->T2D, &d);
+ double beg, end;
+ List_Read(a->T2D, i+3, &beg);
+ if(i > List_Nbr(a->T2D) - 8)
+ end = (double)List_Nbr(a->T2C);
+ else
+ List_Read(a->T2D, i+3+4, &end);
+ char *c = (char*)List_Pointer(a->T2C, (int)beg);
+ for(int j = 0; j < (int)(end-beg); j++)
+ List_Add(b->T2C, &c[j]);
+ b->NbT2++;
+ }
+ for(int i = 0; i < List_Nbr(a->T3D); i+=5){
+ List_Add(b->T3D, List_Pointer(a->T3D, i));
+ List_Add(b->T3D, List_Pointer(a->T3D, i+1));
+ List_Add(b->T3D, List_Pointer(a->T3D, i+2));
+ List_Add(b->T3D, List_Pointer(a->T3D, i+3));
+ double d = List_Nbr(b->T3C);
+ List_Add(b->T3D, &d);
+ double beg, end;
+ List_Read(a->T3D, i+4, &beg);
+ if(i > List_Nbr(a->T3D) - 10)
+ end = (double)List_Nbr(a->T3C);
+ else
+ List_Read(a->T3D, i+4+5, &end);
+ char *c = (char*)List_Pointer(a->T3C, (int)beg);
+ for(int j = 0; j < (int)(end-beg); j++)
+ List_Add(b->T3C, &c[j]);
+ b->NbT3++;
+ }
+}
+
+static void combine_space(struct nameidx *id, List_T *to_remove)
+{
+ int index;
+
+ // sanity check
+ int nbt = 0;
+ for(int i = 0; i < List_Nbr(id->indices); i++) {
+ List_Read(id->indices, i, &index);
+ Post_View *v = *(Post_View **)List_Pointer(CTX.post.list, index);
+ if(!i){
+ nbt = v->NbTimeStep;
+ }
+ else{
+ if(v->NbTimeStep != nbt){
+ Msg(GERROR, "Cannot combine views having different number of time steps");
+ return;
+ }
+ }
+ }
+
+ Post_View *vm = BeginView(1);
+ for(int i = 0; i < List_Nbr(id->indices); i++) {
+ List_Read(id->indices, i, &index);
+ Post_View *v = *(Post_View **)List_Pointer(CTX.post.list, index);
+ List_Insert(to_remove, &v->Num, fcmp_int);
+ combine(v->SP,vm->SP); vm->NbSP += v->NbSP;
+ combine(v->VP,vm->VP); vm->NbVP += v->NbVP;
+ combine(v->TP,vm->TP); vm->NbTP += v->NbTP;
+ combine(v->SL,vm->SL); vm->NbSL += v->NbSL;
+ combine(v->VL,vm->VL); vm->NbVL += v->NbVL;
+ combine(v->TL,vm->TL); vm->NbTL += v->NbTL;
+ combine(v->ST,vm->ST); vm->NbST += v->NbST;
+ combine(v->VT,vm->VT); vm->NbVT += v->NbVT;
+ combine(v->TT,vm->TT); vm->NbTT += v->NbTT;
+ combine(v->SQ,vm->SQ); vm->NbSQ += v->NbSQ;
+ combine(v->VQ,vm->VQ); vm->NbVQ += v->NbVQ;
+ combine(v->TQ,vm->TQ); vm->NbTQ += v->NbTQ;
+ combine(v->SS,vm->SS); vm->NbSS += v->NbSS;
+ combine(v->VS,vm->VS); vm->NbVS += v->NbVS;
+ combine(v->TS,vm->TS); vm->NbTS += v->NbTS;
+ combine(v->SH,vm->SH); vm->NbSH += v->NbSH;
+ combine(v->VH,vm->VH); vm->NbVH += v->NbVH;
+ combine(v->TH,vm->TH); vm->NbTH += v->NbTH;
+ combine(v->SI,vm->SI); vm->NbSI += v->NbSI;
+ combine(v->VI,vm->VI); vm->NbVI += v->NbVI;
+ combine(v->TI,vm->TI); vm->NbTI += v->NbTI;
+ combine(v->SY,vm->SY); vm->NbSY += v->NbSY;
+ combine(v->VY,vm->VY); vm->NbVY += v->NbVY;
+ combine(v->TY,vm->TY); vm->NbTY += v->NbTY;
+ combine_strings(v,vm);
+ }
+
+#if 0
+ // debug strings:
+ for(int i=0; i<List_Nbr(vm->T2D); i++)
+ printf("%g ", *(double*)List_Pointer(vm->T2D, i));
+ printf("\n");
+ for(int i=0; i<List_Nbr(vm->T2C); i++)
+ printf("%c ", *(char*)List_Pointer(vm->T2C, i));
+ printf("\n");
+#endif
+
+ // finalize
+ char name[256], filename[256], tmp[256];
+ if(!strcmp(id->name, "__all__"))
+ strcpy(tmp, "all");
+ else if(!strcmp(id->name, "__vis__"))
+ strcpy(tmp, "visible");
+ else
+ strcpy(tmp, id->name);
+ sprintf(name, "%s_Combine", tmp);
+ sprintf(filename, "%s_Combine.pos", tmp);
+ EndView(vm, 0, filename, name);
+}
+
+static void combine_time(struct nameidx *id, List_T *to_remove)
+{
+ int index, *nbe=0, *nbe2=0, nbn, nbn2, nbc, nbc2;
+ List_T *list=0, *list2=0;
+
+ if(List_Nbr(id->indices) < 2){
+ return; // nothing to do
+ }
+
+ Post_View *vm = BeginView(1);
+
+ // use the first view as the reference
+ List_Read(id->indices, 0, &index);
+ Post_View *v = *(Post_View **)List_Pointer(CTX.post.list, index);
+ for(int i = 0; i < VIEW_NB_ELEMENT_TYPES; i++){
+ vm->get_raw_data(i, &list, &nbe, &nbc, &nbn);
+ v->get_raw_data(i, &list2, &nbe2, &nbc2, &nbn2);
+ *nbe = *nbe2;
+ }
+ vm->NbT2 = v->NbT2;
+ vm->NbT3 = v->NbT3;
+
+ // merge values for all element types
+ for(int i = 0; i < VIEW_NB_ELEMENT_TYPES; i++){
+ vm->get_raw_data(i, &list, &nbe, &nbc, &nbn);
+ for(int j = 0; j < *nbe; j++){
+ for(int k = 0; k < List_Nbr(id->indices); k++){
+ List_Read(id->indices, k, &index);
+ v = *(Post_View **)List_Pointer(CTX.post.list, index);
+ v->get_raw_data(i, &list2, &nbe2, &nbc2, &nbn2);
+ if(*nbe && *nbe == *nbe2){
+ List_Insert(to_remove, &v->Num, fcmp_int);
+ int nb2 = List_Nbr(list2) / *nbe2;
+ if(!k){
+ // copy coordinates of elm j (we are always here as
+ // expected, since the ref view is the first one!)
+ for(int l = 0; l < 3*nbn2; l++){
+ List_Add(list, List_Pointer(list2, j*nb2+l));
+ }
+ }
+ // copy values of elm j
+ for(int l = 0; l < v->NbTimeStep*nbc2*nbn2; l++){
+ List_Add(list, List_Pointer(list2, j*nb2+3*nbn2+l));
+ }
+ }
+ }
+ }
+ }
+
+ // and a bit of spaghetti code for you now:
+
+ // merge 2d strings
+ for(int j = 0; j < vm->NbT2; j++){
+ for(int k = 0; k < List_Nbr(id->indices); k++){
+ List_Read(id->indices, k, &index);
+ v = *(Post_View **)List_Pointer(CTX.post.list, index);
+ if(vm->NbT2 == v->NbT2){
+ List_Insert(to_remove, &v->Num, fcmp_int);
+ if(!k){
+ // copy coordinates
+ List_Add(vm->T2D, List_Pointer(v->T2D, j*4));
+ List_Add(vm->T2D, List_Pointer(v->T2D, j*4+1));
+ List_Add(vm->T2D, List_Pointer(v->T2D, j*4+2));
+ // index
+ double d = List_Nbr(vm->T2C);
+ List_Add(vm->T2D, &d);
+ }
+ // copy char values
+ double beg, end;
+ List_Read(v->T2D, j*4+3, &beg);
+ if(j == vm->NbT2 - 1)
+ end = (double)List_Nbr(v->T2C);
+ else
+ List_Read(v->T2D, j*4+4+3, &end);
+ char *c = (char*)List_Pointer(v->T2C, (int)beg);
+ for(int l = 0; l < (int)(end-beg); l++)
+ List_Add(vm->T2C, &c[l]);
+ }
+ }
+ }
+
+ // merge 3d strings
+ for(int j = 0; j < vm->NbT3; j++){
+ for(int k = 0; k < List_Nbr(id->indices); k++){
+ List_Read(id->indices, k, &index);
+ v = *(Post_View **)List_Pointer(CTX.post.list, index);
+ if(vm->NbT3 == v->NbT3){
+ List_Insert(to_remove, &v->Num, fcmp_int);
+ if(!k){
+ // copy coordinates
+ List_Add(vm->T3D, List_Pointer(v->T3D, j*5));
+ List_Add(vm->T3D, List_Pointer(v->T3D, j*5+1));
+ List_Add(vm->T3D, List_Pointer(v->T3D, j*5+2));
+ List_Add(vm->T3D, List_Pointer(v->T3D, j*5+3));
+ // index
+ double d = List_Nbr(vm->T3C);
+ List_Add(vm->T3D, &d);
+ }
+ // copy char values
+ double beg, end;
+ List_Read(v->T3D, j*5+4, &beg);
+ if(j == vm->NbT3 - 1)
+ end = (double)List_Nbr(v->T3C);
+ else
+ List_Read(v->T3D, j*5+5+4, &end);
+ char *c = (char*)List_Pointer(v->T3C, (int)beg);
+ for(int l = 0; l < (int)(end-beg); l++)
+ List_Add(vm->T3C, &c[l]);
+ }
+ }
+ }
+
+ // create the time data
+ for(int i = 0; i < List_Nbr(id->indices); i++){
+ List_Read(id->indices, i, &index);
+ v = *(Post_View **)List_Pointer(CTX.post.list, index);
+ for(int j = 0; j < List_Nbr(v->Time); j++){
+ List_Add(vm->Time, List_Pointer(v->Time, j));
+ }
+ }
+
+ // if all the time values are the same, it probably means that the
+ // original views didn't have any time data: let's put some indices,
+ // then
+ double time0 = 0.0;
+ if(List_Nbr(vm->Time)) List_Read(vm->Time, 0, &time0);
+ int nbtime = List_Nbr(vm->Time), ok = 0;
+ for(int i = 1; i < nbtime; i++){
+ if(time0 != *(double*)List_Pointer(vm->Time, i)){
+ ok = 1;
+ break;
+ }
+ }
+ if(!ok){
+ List_Reset(vm->Time);
+ for(int i = 0; i < nbtime; i++){
+ double time = i;
+ List_Add(vm->Time, &time);
+ }
+ }
+
+ // finalize
+ char name[256], filename[256], tmp[256];
+ if(!strcmp(id->name, "__all__"))
+ strcpy(tmp, "all");
+ else if(!strcmp(id->name, "__vis__"))
+ strcpy(tmp, "visible");
+ else
+ strcpy(tmp, id->name);
+ sprintf(name, "%s_Combine", tmp);
+ sprintf(filename, "%s_Combine.pos", tmp);
+ EndView(vm, 0, filename, name);
+}
+
+void CombineViews(int time, int how, int remove)
+{
+ // time==0: combine the elements
+ // time==1: combine the timesteps
+
+ // how==0: try to combine all visible views
+ // how==1: try to combine all views
+ // how==2: try to combine all views having identical names
+
+ List_T *ids = List_Create(10, 10, sizeof(nameidx));
+ List_T *to_remove = List_Create(10, 10, sizeof(int));
+ struct nameidx *pid;
+
+ for(int i = 0; i < List_Nbr(CTX.post.list); i++) {
+ Post_View *v = *(Post_View **) List_Pointer(CTX.post.list, i);
+ if(how || v->Visible) {
+ nameidx id;
+ // this might potentially lead to unwanted results if there are
+ // views named "__all__" or "__vis__", but what the heck...
+ if(how == 2)
+ strcpy(id.name, v->Name);
+ else if(how == 1)
+ strcpy(id.name, "__all__");
+ else
+ strcpy(id.name, "__vis__");
+ if((pid = (struct nameidx*)List_PQuery(ids, &id, fcmp_name))){
+ List_Add(pid->indices, &i);
+ }
+ else{
+ id.indices = List_Create(10, 10, sizeof(int));
+ List_Add(id.indices, &i);
+ List_Add(ids, &id);
+ }
+ }
+ }
+
+ for(int i = 0; i < List_Nbr(ids); i++){
+ pid = (struct nameidx*)List_Pointer(ids, i);
+ if(time)
+ combine_time(pid, to_remove);
+ else
+ combine_space(pid, to_remove);
+ List_Delete(pid->indices);
+ }
+ List_Delete(ids);
+
+ // remove original view?
+ if(remove){
+ for(int i = 0; i < List_Nbr(to_remove); i++){
+ int num;
+ List_Read(to_remove, i, &num);
+ RemoveViewByNumber(num);
+ }
+ }
+ List_Delete(to_remove);
+
+#if defined(HAVE_FLTK)
+ UpdateViewsInGUI();
+#endif
+}
+
+// generic access functions
+
+int Post_View::empty(){
+ if(NbSP || NbVP || NbTP ||
+ NbSL || NbVL || NbTL ||
+ NbST || NbVT || NbTT ||
+ NbSQ || NbVQ || NbTQ ||
+ NbSS || NbVS || NbTS ||
+ NbSH || NbVH || NbTH ||
+ NbSI || NbVI || NbTI ||
+ NbSY || NbVY || NbTY ||
+ NbT2 || NbT3)
+ return 0;
+ else
+ return 1;
+}
+
+void Post_View::get_raw_data(int type, List_T **list, int **nbe, int *nbc, int *nbn){
+ switch(type){
+ case 0 : *list = SP; *nbe = &NbSP; *nbc = 1; *nbn = 1; break;
+ case 1 : *list = VP; *nbe = &NbVP; *nbc = 3; *nbn = 1; break;
+ case 2 : *list = TP; *nbe = &NbTP; *nbc = 9; *nbn = 1; break;
+ case 3 : *list = SL; *nbe = &NbSL; *nbc = 1; *nbn = 2; break;
+ case 4 : *list = VL; *nbe = &NbVL; *nbc = 3; *nbn = 2; break;
+ case 5 : *list = TL; *nbe = &NbTL; *nbc = 9; *nbn = 2; break;
+ case 6 : *list = ST; *nbe = &NbST; *nbc = 1; *nbn = 3; break;
+ case 7 : *list = VT; *nbe = &NbVT; *nbc = 3; *nbn = 3; break;
+ case 8 : *list = TT; *nbe = &NbTT; *nbc = 9; *nbn = 3; break;
+ case 9 : *list = SQ; *nbe = &NbSQ; *nbc = 1; *nbn = 4; break;
+ case 10: *list = VQ; *nbe = &NbVQ; *nbc = 3; *nbn = 4; break;
+ case 11: *list = TQ; *nbe = &NbTQ; *nbc = 9; *nbn = 4; break;
+ case 12: *list = SS; *nbe = &NbSS; *nbc = 1; *nbn = 4; break;
+ case 13: *list = VS; *nbe = &NbVS; *nbc = 3; *nbn = 4; break;
+ case 14: *list = TS; *nbe = &NbTS; *nbc = 9; *nbn = 4; break;
+ case 15: *list = SH; *nbe = &NbSH; *nbc = 1; *nbn = 8; break;
+ case 16: *list = VH; *nbe = &NbVH; *nbc = 3; *nbn = 8; break;
+ case 17: *list = TH; *nbe = &NbTH; *nbc = 9; *nbn = 8; break;
+ case 18: *list = SI; *nbe = &NbSI; *nbc = 1; *nbn = 6; break;
+ case 19: *list = VI; *nbe = &NbVI; *nbc = 3; *nbn = 6; break;
+ case 20: *list = TI; *nbe = &NbTI; *nbc = 9; *nbn = 6; break;
+ case 21: *list = SY; *nbe = &NbSY; *nbc = 1; *nbn = 5; break;
+ case 22: *list = VY; *nbe = &NbVY; *nbc = 3; *nbn = 5; break;
+ case 23: *list = TY; *nbe = &NbTY; *nbc = 9; *nbn = 5; break;
+ default: Msg(GERROR, "Wrong type in Post_View::get_raw_data"); break;
+ }
+}
+
+// Generalized raise
+
+void InitGeneralizedRaise(Post_View *v)
+{
+ FreeGeneralizedRaise(v);
+
+ char *expr[3] = { v->GenRaiseX, v->GenRaiseY, v->GenRaiseZ };
+#if defined(HAVE_MATH_EVAL)
+ for(int i = 0; i < 3; i++) {
+ if(strlen(expr[i])) {
+ if(!(v->GenRaise_f[i] = evaluator_create(expr[i])))
+ Msg(GERROR, "Invalid expression '%s'", expr[i]);
+ }
+ }
+#else
+ for(int i = 0; i < 3; i++) {
+ if(!strcmp(expr[i], "v0")) v->GenRaise_f[i] = (void*)0;
+ else if(!strcmp(expr[i], "v1")) v->GenRaise_f[i] = (void*)1;
+ else if(!strcmp(expr[i], "v2")) v->GenRaise_f[i] = (void*)2;
+ else if(!strcmp(expr[i], "v3")) v->GenRaise_f[i] = (void*)3;
+ else if(!strcmp(expr[i], "v4")) v->GenRaise_f[i] = (void*)4;
+ else if(!strcmp(expr[i], "v5")) v->GenRaise_f[i] = (void*)5;
+ else if(!strcmp(expr[i], "v6")) v->GenRaise_f[i] = (void*)6;
+ else if(!strcmp(expr[i], "v7")) v->GenRaise_f[i] = (void*)7;
+ else if(!strcmp(expr[i], "v8")) v->GenRaise_f[i] = (void*)8;
+ else if(strlen(expr[i])) {
+ Msg(GERROR, "Invalid expression '%s'", expr[i]);
+ return;
+ }
+ }
+#endif
+}
+
+void FreeGeneralizedRaise(Post_View *v)
+{
+ for(int i = 0; i < 3; i++){
+#if defined(HAVE_MATH_EVAL)
+ if(v->GenRaise_f[i])
+ evaluator_destroy(v->GenRaise_f[i]);
+ v->GenRaise_f[i] = NULL;
+#else
+ v->GenRaise_f[i] = (void*)-1;
+#endif
+ }
+}
+
+void ApplyGeneralizedRaise(Post_View * v, int numNodes, int numComp, double *vals,
+ double *x, double *y, double *z)
+{
+ double *coords[3] = { x, y, z };
+
+ for(int k = 0; k < numNodes; k++) {
+ double d[9] = {0., 0., 0., 0., 0., 0., 0., 0., 0.};
+ for(int l = 0; l < numComp; l++)
+ d[l] = vals[numComp * k + l];
+#if defined(HAVE_MATH_EVAL)
+ char *names[] = { "v0", "v1", "v2", "v3", "v4", "v5", "v6", "v7", "v8" ,
+ "x", "y", "z" };
+ double values[] = { d[0], d[1], d[2], d[3], d[4], d[5], d[6], d[7], d[8],
+ x[k], y[k], z[k] };
+ for(int i = 0; i < 3; i++) {
+ if(v->GenRaise_f[i])
+ coords[i][k] +=
+ evaluator_evaluate(v->GenRaise_f[i], sizeof(names) /
+ sizeof(names[0]), names, values) * v->GenRaiseFactor;
+ }
+#else
+ for(int i = 0; i < 3; i++){
+ int comp = (int)v->GenRaise_f[i];
+ if(comp >= 0)
+ coords[i][k] += d[comp] * v->GenRaiseFactor;
+ }
+#endif
+ }
+}
diff --git a/Common/Views.h b/Common/Views.h
new file mode 100644
index 0000000000000000000000000000000000000000..3ba4b108997e317dbcd5857126ad7df4d5125bbd
--- /dev/null
+++ b/Common/Views.h
@@ -0,0 +1,231 @@
+#ifndef _VIEWS_H_
+#define _VIEWS_H_
+
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include <map>
+#include "ColorTable.h"
+#include "List.h"
+#include "VertexArray.h"
+#include "SmoothData.h"
+#include "AdaptiveViews.h"
+
+#define VIEW_NB_ELEMENT_TYPES (8*3)
+#define VAL_INF 1.e200
+
+#define POST_POINT 0
+#define POST_LINE 1
+#define POST_TRIANGLE 2
+#define POST_QUADRANGLE 3
+#define POST_TETRAHEDRON 4
+#define POST_HEXAHEDRON 5
+#define POST_PRISM 6
+#define POST_PYRAMID 7
+
+class Post_View{
+ public :
+ // intrinsic to a view
+ int Num, Index, Changed, AliasOf, Links, Dirty;
+ char FileName[256], Name[256];
+
+ // the data
+ int DataSize; // size(double) or sizeof(float)
+ int NbTimeStep, ScalarOnly, TextOnly;
+ double Min, Max, BBox[6], *TimeStepMin, *TimeStepMax;
+ List_T *Time;
+ int NbSP, NbVP, NbTP;
+ List_T *SP, *VP, *TP; // points
+ int NbSL, NbVL, NbTL, NbSL2, NbVL2, NbTL2;
+ List_T *SL, *VL, *TL, *SL2, *VL2, *TL2; // lines
+ int NbST, NbVT, NbTT, NbST2, NbVT2, NbTT2;
+ List_T *ST, *VT, *TT, *ST2, *VT2, *TT2; // triangles
+ int NbSQ, NbVQ, NbTQ, NbSQ2, NbVQ2, NbTQ2;
+ List_T *SQ, *VQ, *TQ, *SQ2, *VQ2, *TQ2; // quadrangles
+ int NbSS, NbVS, NbTS, NbSS2, NbVS2, NbTS2;
+ List_T *SS, *VS, *TS, *SS2, *VS2, *TS2; // tetrahedra
+ int NbSH, NbVH, NbTH, NbSH2, NbVH2, NbTH2;
+ List_T *SH, *VH, *TH, *SH2, *VH2, *TH2; // hexahedra
+ int NbSI, NbVI, NbTI, NbSI2, NbVI2, NbTI2;
+ List_T *SI, *VI, *TI, *SI2, *VI2, *TI2; // prisms
+ int NbSY, NbVY, NbTY, NbSY2, NbVY2, NbTY2;
+ List_T *SY, *VY, *TY, *SY2, *VY2, *TY2; // pyramids
+ int NbT2, NbT3;
+ List_T *T2D, *T2C, *T3D, *T3C; // 2D and 3D text strings
+ std::map < int , List_T * > *Grains; // For LMGC90, grains shapes
+ std::map < int , int > *DisplayListsOfGrains; // For LMGC90, grains shapes
+
+
+ // vertex arrays to draw triangles and lines efficiently
+ VertexArray *TriVertexArray, *LinVertexArray;
+
+ // options
+ int Type;
+ int Position[2], Size[2], AutoPosition;
+ char Format[256];
+ int Axes, AxesAutoPosition, AxesTics[3];
+ char AxesFormat[3][256], AxesLabel[3][256];
+ double AxesPosition[6];
+ double CustomMin, CustomMax;
+ double Offset[3], Raise[3], Transform[3][3], DisplacementFactor, Explode;
+ double ArrowSize, ArrowRelHeadRadius, ArrowRelStemRadius, ArrowRelStemLength;
+ double Normals, Tangents;
+ int Visible, IntervalsType, NbIso, ArrowSizeProportional;
+ int Light, LightTwoSide, LightLines, SmoothNormals;
+ double AngleSmoothNormals;
+ int SaturateValues, FakeTransparency;
+ int ShowElement, ShowTime, ShowScale;
+ int ScaleType, RangeType;
+ int VectorType, TensorType, GlyphLocation;
+ int TimeStep;
+ int DrawStrings;
+ int DrawPoints, DrawLines, DrawTriangles, DrawQuadrangles;
+ int DrawTetrahedra, DrawHexahedra, DrawPrisms, DrawPyramids;
+ int DrawScalars, DrawVectors, DrawTensors;
+ int Boundary, PointType, LineType;
+ double PointSize, LineWidth;
+ GmshColorTable CT;
+ int UseStipple, Stipple[10][2];
+ char StippleString[10][32];
+ int ExternalViewIndex, ViewIndexForGenRaise;
+ int UseGenRaise;
+ double GenRaiseFactor;
+ char GenRaiseX[256], GenRaiseY[256], GenRaiseZ[256];
+ void *GenRaise_f[3];
+
+ // color options
+ struct{
+ unsigned int point, line, triangle, quadrangle;
+ unsigned int tetrahedron, hexahedron, prism, pyramid;
+ unsigned int tangents, normals;
+ unsigned int text2d, text3d, axes;
+ } color;
+
+ // dynamic
+ double (*GVFI)(double min, double max, int nb, int index);
+ int (*GIFV)(double min, double max, int nb, double value);
+ int ExternalElementIndex;
+ double ExternalMin, ExternalMax;
+ double TmpBBox[6]; // dynamically computed during drawing
+
+ // smooth the view
+ void smooth();
+ // smooth normals
+ smooth_normals *normals;
+ void reset_normals();
+
+ // some generic access functions
+ int empty();
+ void get_raw_data(int type, List_T **list, int **nbe, int *nbc, int *nbn);
+
+ // transforms curved elements into linear ones
+ void splitCurvedElements();
+
+ // If the view is high order, coeffs are interpreated as
+ // coefficients of a high order interpolation. So, a pre-pro
+ // is done at the end of the view that sets the view to the
+ // minimal resolution. Then, we can interactively modify the
+ // resolution.
+ Adaptive_Post_View *adaptive;
+ void setGlobalResolutionLevel(int level)
+ {
+ if (adaptive)
+ adaptive->setGlobalResolutionLevel(this, level);
+ }
+ void setAdaptiveResolutionLevel (int level, GMSH_Post_Plugin *plug = 0)
+ {
+ if (adaptive)
+ adaptive->setAdaptiveResolutionLevel(this, level, plug);
+ }
+
+};
+
+// Type
+#define DRAW_POST_3D 1
+#define DRAW_POST_2D_SPACE 2
+#define DRAW_POST_2D_TIME 3
+
+// IntervalsType
+#define DRAW_POST_ISO 1
+#define DRAW_POST_CONTINUOUS 2
+#define DRAW_POST_DISCRETE 3
+#define DRAW_POST_NUMERIC 4
+
+// VectorType
+#define DRAW_POST_SEGMENT 1
+#define DRAW_POST_ARROW 2
+#define DRAW_POST_PYRAMID 3
+#define DRAW_POST_ARROW3D 4
+#define DRAW_POST_DISPLACEMENT 5
+
+// GlyphLocation
+#define DRAW_POST_LOCATE_COG 1
+#define DRAW_POST_LOCATE_VERTEX 2
+
+// TensorType
+#define DRAW_POST_VONMISES 1
+#define DRAW_POST_LMGC90 2
+#define DRAW_POST_LMGC90_TYPE 3
+#define DRAW_POST_LMGC90_COORD 4
+#define DRAW_POST_LMGC90_PRES 5
+#define DRAW_POST_LMGC90_SN 6
+#define DRAW_POST_LMGC90_DEPX 7
+#define DRAW_POST_LMGC90_DEPY 8
+#define DRAW_POST_LMGC90_DEPZ 9
+#define DRAW_POST_LMGC90_DEPAV 10
+#define DRAW_POST_LMGC90_DEPNORM 11
+// RangeType
+#define DRAW_POST_RANGE_DEFAULT 1
+#define DRAW_POST_RANGE_CUSTOM 2
+#define DRAW_POST_RANGE_PER_STEP 3
+
+// ScaleType
+#define DRAW_POST_LINEAR 1
+#define DRAW_POST_LOGARITHMIC 2
+#define DRAW_POST_DOUBLELOGARITHMIC 3 // for vorticity e.g.
+
+// Public functions
+
+int fcmpPostViewNum(const void *v1, const void *v2);
+int fcmpPostViewAliasOf(const void *v1, const void *v2);
+
+Post_View *BeginView(int alloc);
+void EndView(Post_View *v, int AddInUI, char *FileName, char *Name);
+void AliasView(int num, int withoptions);
+void FreeView(Post_View *v);
+bool RemoveViewByIndex(int index);
+bool RemoveViewByNumber(int num);
+int ReadView(char *filename);
+void WriteView(Post_View *v, char *filename, int format, int append);
+void CopyViewOptions(Post_View *src, Post_View *dest);
+void CombineViews(int time, int how, int remove);
+
+Post_View *Create2DGraph(char *xname, char *yname, int nbdata, double *x, double *y);
+
+GmshColorTable *Get_ColorTable(int num);
+void Print_ColorTable(int num, int diff, char *prefix, FILE *file);
+
+double ComputeVonMises(double* val);
+
+void InitGeneralizedRaise(Post_View *v);
+void FreeGeneralizedRaise(Post_View *v);
+void ApplyGeneralizedRaise(Post_View * v, int numNodes, int numComp, double *vals,
+ double *x, double *y, double *z);
+
+#endif
diff --git a/Common/ViewsIO.cpp b/Common/ViewsIO.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..860e6e7ae3cf27ee0187e0a638e7c221d817364a
--- /dev/null
+++ b/Common/ViewsIO.cpp
@@ -0,0 +1,839 @@
+// $Id: ViewsIO.cpp,v 1.9 2006-11-27 22:22:08 geuzaine Exp $
+//
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include <set>
+#include "Gmsh.h"
+#include "Numeric.h"
+#include "Views.h"
+#include "Context.h"
+
+extern Context_T CTX;
+
+#if defined(HAVE_FLTK)
+void UpdateViewsInGUI();
+#endif
+
+// Read view from file
+
+int ReadView(char *filename)
+{
+ FILE *fp = fopen(filename, "rb");
+ if(!fp){
+ Msg(GERROR, "Unable to open file '%s'", filename);
+ return 0;
+ }
+
+ char str[256], name[256];
+ int i, nb, format, size, testone, swap, t2l, t3l;
+ double version;
+ Post_View *v;
+
+ while(1) {
+
+ do {
+ if(!fgets(str, 256, fp) || feof(fp))
+ break;
+ } while(str[0] != '$');
+
+ if(feof(fp))
+ break;
+
+ if(!strncmp(&str[1], "PostFormat", 10)) {
+ if(!fscanf(fp, "%lf %d %d\n", &version, &format, &size)){
+ Msg(GERROR, "Read error");
+ return 0;
+ }
+ if(version < 1.0) {
+ Msg(GERROR, "This post-processing file is too old (version %g < 1.0)",
+ version);
+ return 0;
+ }
+ if(size == sizeof(double))
+ Msg(DEBUG, "Data is in double precision format (size==%d)", size);
+ else if(size == sizeof(float))
+ Msg(DEBUG, "Data is in single precision format (size==%d)", size);
+ else {
+ Msg(GERROR, "Unknown data size (%d) in post-processing file", size);
+ return 0;
+ }
+ if(format == 0)
+ format = LIST_FORMAT_ASCII;
+ else if(format == 1)
+ format = LIST_FORMAT_BINARY;
+ else {
+ Msg(GERROR, "Unknown format for view");
+ return 0;
+ }
+ }
+
+ if(!strncmp(&str[1], "View", 4)) {
+ v = BeginView(0);
+ if(version <= 1.0) {
+ Msg(DEBUG, "Detected post-processing view format <= 1.0");
+ if(!fscanf(fp, "%s %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
+ name, &v->NbTimeStep,
+ &v->NbSP, &v->NbVP, &v->NbTP, &v->NbSL, &v->NbVL, &v->NbTL,
+ &v->NbST, &v->NbVT, &v->NbTT, &v->NbSS, &v->NbVS, &v->NbTS)){
+ Msg(GERROR, "Read error");
+ return 0;
+ }
+ v->NbT2 = t2l = v->NbT3 = t3l = 0;
+ }
+ else if(version == 1.1) {
+ Msg(DEBUG, "Detected post-processing view format 1.1");
+ if(!fscanf(fp,
+ "%s %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
+ name, &v->NbTimeStep, &v->NbSP, &v->NbVP, &v->NbTP, &v->NbSL,
+ &v->NbVL, &v->NbTL, &v->NbST, &v->NbVT, &v->NbTT, &v->NbSS,
+ &v->NbVS, &v->NbTS, &v->NbT2, &t2l, &v->NbT3, &t3l)){
+ Msg(GERROR, "Read error");
+ return 0;
+ }
+ }
+ else if(version == 1.2 || version == 1.3) {
+ Msg(DEBUG, "Detected post-processing view format %g", version);
+ if(!fscanf(fp, "%s %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d "
+ "%d %d %d %d %d %d %d %d %d %d %d %d %d\n",
+ name, &v->NbTimeStep,
+ &v->NbSP, &v->NbVP, &v->NbTP, &v->NbSL, &v->NbVL, &v->NbTL,
+ &v->NbST, &v->NbVT, &v->NbTT, &v->NbSQ, &v->NbVQ, &v->NbTQ,
+ &v->NbSS, &v->NbVS, &v->NbTS, &v->NbSH, &v->NbVH, &v->NbTH,
+ &v->NbSI, &v->NbVI, &v->NbTI, &v->NbSY, &v->NbVY, &v->NbTY,
+ &v->NbT2, &t2l, &v->NbT3, &t3l)){
+ Msg(GERROR, "Read error");
+ return 0;
+ }
+ }
+ else if(version == 1.4) {
+ Msg(DEBUG, "Detected post-processing view format 1.4");
+ if(!fscanf(fp, "%s %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d "
+ "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d "
+ "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
+ name, &v->NbTimeStep,
+ &v->NbSP, &v->NbVP, &v->NbTP, &v->NbSL, &v->NbVL, &v->NbTL,
+ &v->NbST, &v->NbVT, &v->NbTT, &v->NbSQ, &v->NbVQ, &v->NbTQ,
+ &v->NbSS, &v->NbVS, &v->NbTS, &v->NbSH, &v->NbVH, &v->NbTH,
+ &v->NbSI, &v->NbVI, &v->NbTI, &v->NbSY, &v->NbVY, &v->NbTY,
+ &v->NbSL2, &v->NbVL2, &v->NbTL2, &v->NbST2, &v->NbVT2, &v->NbTT2,
+ &v->NbSQ2, &v->NbVQ2, &v->NbTQ2, &v->NbSS2, &v->NbVS2, &v->NbTS2,
+ &v->NbSH2, &v->NbVH2, &v->NbTH2, &v->NbSI2, &v->NbVI2, &v->NbTI2,
+ &v->NbSY2, &v->NbVY2, &v->NbTY2, &v->NbT2, &t2l, &v->NbT3, &t3l)){
+ Msg(GERROR, "Read error");
+ return 0;
+ }
+ }
+ else {
+ Msg(GERROR, "Unknown post-processing file format (version %g)",
+ version);
+ return 0;
+ }
+
+ for(i = 0; i < (int)strlen(name); i++)
+ if(name[i] == '^')
+ name[i] = ' ';
+
+ swap = 0;
+ if(format == LIST_FORMAT_BINARY) {
+ if(!fread(&testone, sizeof(int), 1, fp)){
+ Msg(GERROR, "Read error");
+ return 0;
+ }
+ if(testone != 1) {
+ Msg(INFO, "Swapping bytes from binary file");
+ swap = 1;
+ }
+ }
+
+ v->DataSize = size;
+
+ // Time values
+ v->Time = List_CreateFromFile(v->NbTimeStep, 100, size, fp, format, swap);
+
+ // Note: if nb==0, we still allocates the lists (so that they
+ // are ready to be filled later, e.g. in plugins)
+
+#define LCD List_CreateFromFile(nb, 1000, size, fp, format, swap)
+
+ // Points
+ nb = v->NbSP ? v->NbSP * (v->NbTimeStep * 1 + 3) : 0; v->SP = LCD;
+ nb = v->NbVP ? v->NbVP * (v->NbTimeStep * 3 + 3) : 0; v->VP = LCD;
+ nb = v->NbTP ? v->NbTP * (v->NbTimeStep * 9 + 3) : 0; v->TP = LCD;
+
+ // Lines
+ nb = v->NbSL ? v->NbSL * (v->NbTimeStep * 2 * 1 + 6) : 0; v->SL = LCD;
+ nb = v->NbVL ? v->NbVL * (v->NbTimeStep * 2 * 3 + 6) : 0; v->VL = LCD;
+ nb = v->NbTL ? v->NbTL * (v->NbTimeStep * 2 * 9 + 6) : 0; v->TL = LCD;
+
+ // Triangles
+ nb = v->NbST ? v->NbST * (v->NbTimeStep * 3 * 1 + 9) : 0; v->ST = LCD;
+ nb = v->NbVT ? v->NbVT * (v->NbTimeStep * 3 * 3 + 9) : 0; v->VT = LCD;
+ nb = v->NbTT ? v->NbTT * (v->NbTimeStep * 3 * 9 + 9) : 0; v->TT = LCD;
+
+ // Quadrangles
+ nb = v->NbSQ ? v->NbSQ * (v->NbTimeStep * 4 * 1 + 12) : 0; v->SQ = LCD;
+ nb = v->NbVQ ? v->NbVQ * (v->NbTimeStep * 4 * 3 + 12) : 0; v->VQ = LCD;
+ nb = v->NbTQ ? v->NbTQ * (v->NbTimeStep * 4 * 9 + 12) : 0; v->TQ = LCD;
+
+ // Tetrahedra
+ nb = v->NbSS ? v->NbSS * (v->NbTimeStep * 4 * 1 + 12) : 0; v->SS = LCD;
+ nb = v->NbVS ? v->NbVS * (v->NbTimeStep * 4 * 3 + 12) : 0; v->VS = LCD;
+ nb = v->NbTS ? v->NbTS * (v->NbTimeStep * 4 * 9 + 12) : 0; v->TS = LCD;
+
+ // Hexahedra
+ nb = v->NbSH ? v->NbSH * (v->NbTimeStep * 8 * 1 + 24) : 0; v->SH = LCD;
+ nb = v->NbVH ? v->NbVH * (v->NbTimeStep * 8 * 3 + 24) : 0; v->VH = LCD;
+ nb = v->NbTH ? v->NbTH * (v->NbTimeStep * 8 * 9 + 24) : 0; v->TH = LCD;
+
+ // Prisms
+ nb = v->NbSI ? v->NbSI * (v->NbTimeStep * 6 * 1 + 18) : 0; v->SI = LCD;
+ nb = v->NbVI ? v->NbVI * (v->NbTimeStep * 6 * 3 + 18) : 0; v->VI = LCD;
+ nb = v->NbTI ? v->NbTI * (v->NbTimeStep * 6 * 9 + 18) : 0; v->TI = LCD;
+
+ // Pyramids
+ nb = v->NbSY ? v->NbSY * (v->NbTimeStep * 5 * 1 + 15) : 0; v->SY = LCD;
+ nb = v->NbVY ? v->NbVY * (v->NbTimeStep * 5 * 3 + 15) : 0; v->VY = LCD;
+ nb = v->NbTY ? v->NbTY * (v->NbTimeStep * 5 * 9 + 15) : 0; v->TY = LCD;
+
+ // 2nd order Lines
+ nb = v->NbSL2 ? v->NbSL2 * (v->NbTimeStep * 3 * 1 + 9) : 0; v->SL2 = LCD;
+ nb = v->NbVL2 ? v->NbVL2 * (v->NbTimeStep * 3 * 3 + 9) : 0; v->VL2 = LCD;
+ nb = v->NbTL2 ? v->NbTL2 * (v->NbTimeStep * 3 * 9 + 9) : 0; v->TL2 = LCD;
+
+ // 2nd order Triangles
+ nb = v->NbST2 ? v->NbST2 * (v->NbTimeStep * 6 * 1 + 18) : 0; v->ST2 = LCD;
+ nb = v->NbVT2 ? v->NbVT2 * (v->NbTimeStep * 6 * 3 + 18) : 0; v->VT2 = LCD;
+ nb = v->NbTT2 ? v->NbTT2 * (v->NbTimeStep * 6 * 9 + 18) : 0; v->TT2 = LCD;
+
+ // 2nd order Quadrangles
+ nb = v->NbSQ2 ? v->NbSQ2 * (v->NbTimeStep * 9 * 1 + 27) : 0; v->SQ2 = LCD;
+ nb = v->NbVQ2 ? v->NbVQ2 * (v->NbTimeStep * 9 * 3 + 27) : 0; v->VQ2 = LCD;
+ nb = v->NbTQ2 ? v->NbTQ2 * (v->NbTimeStep * 9 * 9 + 27) : 0; v->TQ2 = LCD;
+
+ // 2nd order Tetrahedra
+ nb = v->NbSS2 ? v->NbSS2 * (v->NbTimeStep * 10 * 1 + 30) : 0; v->SS2 = LCD;
+ nb = v->NbVS2 ? v->NbVS2 * (v->NbTimeStep * 10 * 3 + 30) : 0; v->VS2 = LCD;
+ nb = v->NbTS2 ? v->NbTS2 * (v->NbTimeStep * 10 * 9 + 30) : 0; v->TS2 = LCD;
+
+ // 2nd order Hexahedra
+ nb = v->NbSH2 ? v->NbSH2 * (v->NbTimeStep * 27 * 1 + 81) : 0; v->SH2 = LCD;
+ nb = v->NbVH2 ? v->NbVH2 * (v->NbTimeStep * 27 * 3 + 81) : 0; v->VH2 = LCD;
+ nb = v->NbTH2 ? v->NbTH2 * (v->NbTimeStep * 27 * 9 + 81) : 0; v->TH2 = LCD;
+
+ // 2nd order Prisms
+ nb = v->NbSI2 ? v->NbSI2 * (v->NbTimeStep * 18 * 1 + 54) : 0; v->SI2 = LCD;
+ nb = v->NbVI2 ? v->NbVI2 * (v->NbTimeStep * 18 * 3 + 54) : 0; v->VI2 = LCD;
+ nb = v->NbTI2 ? v->NbTI2 * (v->NbTimeStep * 18 * 9 + 54) : 0; v->TI2 = LCD;
+
+ // 2nd order Pyramids
+ nb = v->NbSY2 ? v->NbSY2 * (v->NbTimeStep * 14 * 1 + 42) : 0; v->SY2 = LCD;
+ nb = v->NbVY2 ? v->NbVY2 * (v->NbTimeStep * 14 * 3 + 42) : 0; v->VY2 = LCD;
+ nb = v->NbTY2 ? v->NbTY2 * (v->NbTimeStep * 14 * 9 + 42) : 0; v->TY2 = LCD;
+
+#undef LCD
+
+ // 2D strings
+ nb = v->NbT2 ? v->NbT2 * 4 : 0;
+ v->T2D = List_CreateFromFile(nb, 100, size, fp, format, swap);
+ if(version <= 1.2)
+ v->T2C = List_CreateFromFileOld(t2l, 100, sizeof(char), fp, format, swap);
+ else
+ v->T2C = List_CreateFromFile(t2l, 100, sizeof(char), fp, format, swap);
+
+ // 3D strings
+ nb = v->NbT3 ? v->NbT3 * 5 : 0;
+ v->T3D = List_CreateFromFile(nb, 100, size, fp, format, swap);
+ if(version <= 1.2)
+ v->T3C = List_CreateFromFileOld(t3l, 100, sizeof(char), fp, format, swap);
+ else
+ v->T3C = List_CreateFromFile(t3l, 100, sizeof(char), fp, format, swap);
+
+ Msg(DEBUG,
+ "Read View '%s' (%d TimeSteps): "
+ "SP(%d/%d) VP(%d/%d) TP(%d/%d) "
+ "SL(%d/%d) VL(%d/%d) TL(%d/%d) "
+ "ST(%d/%d) VT(%d/%d) TT(%d/%d) "
+ "SQ(%d/%d) VQ(%d/%d) TQ(%d/%d) "
+ "SS(%d/%d) VS(%d/%d) TS(%d/%d) "
+ "SH(%d/%d) VH(%d/%d) TH(%d/%d) "
+ "SI(%d/%d) VI(%d/%d) TI(%d/%d) "
+ "SY(%d/%d) VY(%d/%d) TY(%d/%d) "
+ "SL2(%d/%d) VL2(%d/%d) TL2(%d/%d) "
+ "ST2(%d/%d) VT2(%d/%d) TT2(%d/%d) "
+ "SQ2(%d/%d) VQ2(%d/%d) TQ2(%d/%d) "
+ "SS2(%d/%d) VS2(%d/%d) TS2(%d/%d) "
+ "SH2(%d/%d) VH2(%d/%d) TH2(%d/%d) "
+ "SI2(%d/%d) VI2(%d/%d) TI2(%d/%d) "
+ "SY2(%d/%d) VY2(%d/%d) TY2(%d/%d) "
+ "T2(%d/%d/%d) T3(%d/%d/%d) ",
+ name, v->NbTimeStep,
+ v->NbSP, List_Nbr(v->SP), v->NbVP, List_Nbr(v->VP), v->NbTP, List_Nbr(v->TP),
+ v->NbSL, List_Nbr(v->SL), v->NbVL, List_Nbr(v->VL), v->NbTL, List_Nbr(v->TL),
+ v->NbST, List_Nbr(v->ST), v->NbVT, List_Nbr(v->VT), v->NbTT, List_Nbr(v->TT),
+ v->NbSQ, List_Nbr(v->SQ), v->NbVQ, List_Nbr(v->VQ), v->NbTQ, List_Nbr(v->TQ),
+ v->NbSS, List_Nbr(v->SS), v->NbVS, List_Nbr(v->VS), v->NbTS, List_Nbr(v->TS),
+ v->NbSH, List_Nbr(v->SH), v->NbVH, List_Nbr(v->VH), v->NbTH, List_Nbr(v->TH),
+ v->NbSI, List_Nbr(v->SI), v->NbVI, List_Nbr(v->VI), v->NbTI, List_Nbr(v->TI),
+ v->NbSY, List_Nbr(v->SY), v->NbVY, List_Nbr(v->VY), v->NbTY, List_Nbr(v->TY),
+ v->NbSL2, List_Nbr(v->SL2), v->NbVL2, List_Nbr(v->VL2), v->NbTL2, List_Nbr(v->TL2),
+ v->NbST2, List_Nbr(v->ST2), v->NbVT2, List_Nbr(v->VT2), v->NbTT2, List_Nbr(v->TT2),
+ v->NbSQ2, List_Nbr(v->SQ2), v->NbVQ2, List_Nbr(v->VQ2), v->NbTQ2, List_Nbr(v->TQ2),
+ v->NbSS2, List_Nbr(v->SS2), v->NbVS2, List_Nbr(v->VS2), v->NbTS2, List_Nbr(v->TS2),
+ v->NbSH2, List_Nbr(v->SH2), v->NbVH2, List_Nbr(v->VH2), v->NbTH2, List_Nbr(v->TH2),
+ v->NbSI2, List_Nbr(v->SI2), v->NbVI2, List_Nbr(v->VI2), v->NbTI2, List_Nbr(v->TI2),
+ v->NbSY2, List_Nbr(v->SY2), v->NbVY2, List_Nbr(v->VY2), v->NbTY2, List_Nbr(v->TY2),
+ v->NbT2, List_Nbr(v->T2D), List_Nbr(v->T2C),
+ v->NbT3, List_Nbr(v->T3D), List_Nbr(v->T3C));
+
+ // don't update the ui after each view, but only at the end
+ EndView(v, 0, filename, name);
+ }
+
+ do {
+ if(!fgets(str, 256, fp) || feof(fp))
+ Msg(GERROR, "Prematured end of file");
+ } while(str[0] != '$');
+
+ }
+
+#if defined(HAVE_FLTK)
+ UpdateViewsInGUI();
+#endif
+
+ return 1;
+}
+
+// Write view to file in Parsed, ASCII or Binary format
+
+static void write_parsed_time(List_T *list, FILE *fp)
+{
+ if(List_Nbr(list) > 1) {
+ fprintf(fp, "TIME{");
+ for(int i = 0; i < List_Nbr(list); i ++) {
+ if(i) fprintf(fp, ",");
+ fprintf(fp, "%.16g", *(double *)List_Pointer(list, i));
+ }
+ fprintf(fp, "};\n");
+ }
+}
+
+static void write_parsed_elements(char *str, int nbnod, int nb, List_T *list, FILE *fp)
+{
+ if(nb) {
+ int n = List_Nbr(list) / nb;
+ for(int i = 0; i < List_Nbr(list); i += n) {
+ double *x = (double *)List_Pointer(list, i);
+ double *y = (double *)List_Pointer(list, i + nbnod);
+ double *z = (double *)List_Pointer(list, i + 2 * nbnod);
+ fprintf(fp, "%s(", str);
+ for(int j = 0; j < nbnod; j++) {
+ if(j) fprintf(fp, ",");
+ fprintf(fp, "%.16g,%.16g,%.16g", x[j], y[j], z[j]);
+ }
+ fprintf(fp, "){");
+ for(int j = 3 * nbnod; j < n; j++) {
+ if(j - 3 * nbnod) fprintf(fp, ",");
+ fprintf(fp, "%.16g", *(double *)List_Pointer(list, i + j));
+ }
+ fprintf(fp, "};\n");
+ }
+ }
+}
+
+static void write_parsed_strings(int nbc, int nb, List_T *TD, List_T *TC, FILE *fp)
+{
+ if(!nb || (nbc != 4 && nbc != 5)) return;
+ for(int j = 0; j < List_Nbr(TD); j += nbc){
+ double x, y, z, style, start, end;
+ List_Read(TD, j, &x);
+ List_Read(TD, j+1, &y);
+ if(nbc == 5)
+ List_Read(TD, j+2, &z);
+ List_Read(TD, j+nbc-2, &style);
+ if(nbc == 4)
+ fprintf(fp, "T2(%g,%g,%g){", x, y, style);
+ else
+ fprintf(fp, "T3(%g,%g,%g,%g){", x, y, z, style);
+ List_Read(TD, j+nbc-1, &start);
+ if(j+nbc*2-1 < List_Nbr(TD))
+ List_Read(TD, j+nbc*2-1, &end);
+ else
+ end = List_Nbr(TC);
+ int l = 0;
+ while(l < end-start){
+ char *str = (char*)List_Pointer(TC, (int)start + l);
+ if(l) fprintf(fp, ",");
+ fprintf(fp, "\"%s\"", str);
+ l += strlen(str)+1;
+ }
+ fprintf(fp, "};\n");
+ }
+}
+
+void WriteViewPOS(Post_View *v, FILE *file, int binary=0, int parsed=1, int append=0)
+{
+ char name[256];
+ int f, One = 1;
+
+ if(!parsed && !append){
+ fprintf(file, "$PostFormat /* Gmsh 1.3, %s */\n",
+ binary ? "binary" : "ascii");
+ fprintf(file, "1.3 %d %d\n", binary, (int)sizeof(double));
+ fprintf(file, "$EndPostFormat\n");
+ }
+
+ strcpy(name, v->Name);
+ for(int i = 0; i < (int)strlen(name); i++)
+ if(name[i] == ' ') name[i] = '^';
+
+ if(!parsed){
+ fprintf(file, "$View /* %s */\n", v->Name);
+ fprintf(file, "%s ", name);
+ fprintf(file, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d "
+ "%d %d %d %d %d %d %d %d %d %d %d %d\n",
+ List_Nbr(v->Time),
+ v->NbSP, v->NbVP, v->NbTP, v->NbSL, v->NbVL, v->NbTL,
+ v->NbST, v->NbVT, v->NbTT, v->NbSQ, v->NbVQ, v->NbTQ,
+ v->NbSS, v->NbVS, v->NbTS, v->NbSH, v->NbVH, v->NbTH,
+ v->NbSI, v->NbVI, v->NbTI, v->NbSY, v->NbVY, v->NbTY,
+ v->NbT2, List_Nbr(v->T2C), v->NbT3, List_Nbr(v->T3C));
+ if(binary) {
+ f = LIST_FORMAT_BINARY;
+ if(!fwrite(&One, sizeof(int), 1, file)){
+ Msg(GERROR, "Write error");
+ return;
+ }
+ }
+ else
+ f = LIST_FORMAT_ASCII;
+ List_WriteToFile(v->Time, file, f);
+ List_WriteToFile(v->SP, file, f);
+ List_WriteToFile(v->VP, file, f);
+ List_WriteToFile(v->TP, file, f);
+ List_WriteToFile(v->SL, file, f);
+ List_WriteToFile(v->VL, file, f);
+ List_WriteToFile(v->TL, file, f);
+ List_WriteToFile(v->ST, file, f);
+ List_WriteToFile(v->VT, file, f);
+ List_WriteToFile(v->TT, file, f);
+ List_WriteToFile(v->SQ, file, f);
+ List_WriteToFile(v->VQ, file, f);
+ List_WriteToFile(v->TQ, file, f);
+ List_WriteToFile(v->SS, file, f);
+ List_WriteToFile(v->VS, file, f);
+ List_WriteToFile(v->TS, file, f);
+ List_WriteToFile(v->SH, file, f);
+ List_WriteToFile(v->VH, file, f);
+ List_WriteToFile(v->TH, file, f);
+ List_WriteToFile(v->SI, file, f);
+ List_WriteToFile(v->VI, file, f);
+ List_WriteToFile(v->TI, file, f);
+ List_WriteToFile(v->SY, file, f);
+ List_WriteToFile(v->VY, file, f);
+ List_WriteToFile(v->TY, file, f);
+ List_WriteToFile(v->T2D, file, f);
+ List_WriteToFile(v->T2C, file, f);
+ List_WriteToFile(v->T3D, file, f);
+ List_WriteToFile(v->T3C, file, f);
+ fprintf(file, "\n");
+ fprintf(file, "$EndView\n");
+ }
+ else{
+ fprintf(file, "View \"%s\" {\n", v->Name);
+ write_parsed_time(v->Time, file);
+ write_parsed_elements("SP", 1, v->NbSP, v->SP, file);
+ write_parsed_elements("VP", 1, v->NbVP, v->VP, file);
+ write_parsed_elements("TP", 1, v->NbTP, v->TP, file);
+ write_parsed_elements("SL", 2, v->NbSL, v->SL, file);
+ write_parsed_elements("VL", 2, v->NbVL, v->VL, file);
+ write_parsed_elements("TL", 2, v->NbTL, v->TL, file);
+ write_parsed_elements("ST", 3, v->NbST, v->ST, file);
+ write_parsed_elements("VT", 3, v->NbVT, v->VT, file);
+ write_parsed_elements("TT", 3, v->NbTT, v->TT, file);
+ write_parsed_elements("SQ", 4, v->NbSQ, v->SQ, file);
+ write_parsed_elements("VQ", 4, v->NbVQ, v->VQ, file);
+ write_parsed_elements("TQ", 4, v->NbTQ, v->TQ, file);
+ write_parsed_elements("SS", 4, v->NbSS, v->SS, file);
+ write_parsed_elements("VS", 4, v->NbVS, v->VS, file);
+ write_parsed_elements("TS", 4, v->NbTS, v->TS, file);
+ write_parsed_elements("SH", 8, v->NbSH, v->SH, file);
+ write_parsed_elements("VH", 8, v->NbVH, v->VH, file);
+ write_parsed_elements("TH", 8, v->NbTH, v->TH, file);
+ write_parsed_elements("SI", 6, v->NbSI, v->SI, file);
+ write_parsed_elements("VI", 6, v->NbVI, v->VI, file);
+ write_parsed_elements("TI", 6, v->NbTI, v->TI, file);
+ write_parsed_elements("SY", 5, v->NbSY, v->SY, file);
+ write_parsed_elements("VY", 5, v->NbVY, v->VY, file);
+ write_parsed_elements("TY", 5, v->NbTY, v->TY, file);
+ write_parsed_strings(4, v->NbT2, v->T2D, v->T2C, file);
+ write_parsed_strings(5, v->NbT3, v->T3D, v->T3C, file);
+ fprintf(file, "};\n");
+ }
+}
+
+// Write view to file in STL format
+
+static void write_stl(FILE *file, int nbelm, List_T *list, int nbnod)
+{
+ if(nbelm){
+ int nb = List_Nbr(list) / nbelm;
+ for(int i = 0; i < List_Nbr(list); i+=nb){
+ double *x = (double*)List_Pointer(list, i);
+ double n[3];
+ normal3points(x[0], x[3], x[6],
+ x[1], x[4], x[7],
+ x[2], x[5], x[8], n);
+ if(nbnod == 3){
+ fprintf(file, "facet normal %g %g %g\n", n[0], n[1], n[2]);
+ fprintf(file, " outer loop\n");
+ fprintf(file, " vertex %g %g %g\n", x[0], x[3], x[6]);
+ fprintf(file, " vertex %g %g %g\n", x[1], x[4], x[7]);
+ fprintf(file, " vertex %g %g %g\n", x[2], x[5], x[8]);
+ fprintf(file, " endloop\n");
+ fprintf(file, "endfacet\n");
+ }
+ else{
+ fprintf(file, "facet normal %g %g %g\n", n[0], n[1], n[2]);
+ fprintf(file, " outer loop\n");
+ fprintf(file, " vertex %g %g %g\n", x[0], x[4], x[8]);
+ fprintf(file, " vertex %g %g %g\n", x[1], x[5], x[9]);
+ fprintf(file, " vertex %g %g %g\n", x[2], x[6], x[10]);
+ fprintf(file, " endloop\n");
+ fprintf(file, "endfacet\n");
+ fprintf(file, "facet normal %g %g %g\n", n[0], n[1], n[2]);
+ fprintf(file, " outer loop\n");
+ fprintf(file, " vertex %g %g %g\n", x[0], x[4], x[8]);
+ fprintf(file, " vertex %g %g %g\n", x[2], x[6], x[10]);
+ fprintf(file, " vertex %g %g %g\n", x[3], x[7], x[11]);
+ fprintf(file, " endloop\n");
+ fprintf(file, "endfacet\n");
+ }
+ }
+ }
+}
+
+void WriteViewSTL(Post_View *v, FILE *file)
+{
+ if(!v->NbST && !v->NbVT && !v->NbTT &&
+ !v->NbSQ && !v->NbVQ && !v->NbTQ){
+ Msg(GERROR, "No surface elements to save");
+ return;
+ }
+
+ fprintf(file, "solid Created by Gmsh\n");
+ write_stl(file, v->NbST, v->ST, 3);
+ write_stl(file, v->NbVT, v->VT, 3);
+ write_stl(file, v->NbTT, v->TT, 3);
+ write_stl(file, v->NbSQ, v->SQ, 4);
+ write_stl(file, v->NbVQ, v->VQ, 4);
+ write_stl(file, v->NbTQ, v->TQ, 4);
+ fprintf(file, "endsolid Created by Gmsh\n");
+}
+
+// Write view to file in Text format (should change this to have 2
+// choices: "Tabular (By Element)" and "Tabular (By Time Step)")
+
+static void write_txt(FILE *file, int nbelm, List_T *list,
+ int nbnod, int nbcomp, int nbtime)
+{
+ if(nbelm){
+ int nb = List_Nbr(list) / nbelm;
+ for(int i = 0; i < List_Nbr(list); i+=nb){
+ double *x = (double*)List_Pointer(list, i);
+ for(int j = 0; j < nbnod*(3+nbcomp*nbtime); j++)
+ fprintf(file, "%.16g ", x[j]);
+ fprintf(file, "\n");
+ }
+ fprintf(file, "\n");
+ }
+}
+
+void WriteViewTXT(Post_View *v, FILE *file)
+{
+ write_txt(file, v->NbSP, v->SP, 1, 1, v->NbTimeStep);
+ write_txt(file, v->NbVP, v->VP, 1, 3, v->NbTimeStep);
+ write_txt(file, v->NbTP, v->TP, 1, 9, v->NbTimeStep);
+ write_txt(file, v->NbSL, v->SL, 2, 1, v->NbTimeStep);
+ write_txt(file, v->NbVL, v->VL, 2, 3, v->NbTimeStep);
+ write_txt(file, v->NbTL, v->TL, 2, 9, v->NbTimeStep);
+ write_txt(file, v->NbST, v->ST, 3, 1, v->NbTimeStep);
+ write_txt(file, v->NbVT, v->VT, 3, 3, v->NbTimeStep);
+ write_txt(file, v->NbTT, v->TT, 3, 9, v->NbTimeStep);
+ write_txt(file, v->NbSQ, v->SQ, 4, 1, v->NbTimeStep);
+ write_txt(file, v->NbVQ, v->VQ, 4, 3, v->NbTimeStep);
+ write_txt(file, v->NbTQ, v->TQ, 4, 9, v->NbTimeStep);
+ write_txt(file, v->NbSS, v->SS, 4, 1, v->NbTimeStep);
+ write_txt(file, v->NbVS, v->VS, 4, 3, v->NbTimeStep);
+ write_txt(file, v->NbTS, v->TS, 4, 9, v->NbTimeStep);
+ write_txt(file, v->NbSH, v->SH, 8, 1, v->NbTimeStep);
+ write_txt(file, v->NbVH, v->VH, 8, 3, v->NbTimeStep);
+ write_txt(file, v->NbTH, v->TH, 8, 9, v->NbTimeStep);
+ write_txt(file, v->NbSI, v->SI, 6, 1, v->NbTimeStep);
+ write_txt(file, v->NbVI, v->VI, 6, 3, v->NbTimeStep);
+ write_txt(file, v->NbTI, v->TI, 6, 9, v->NbTimeStep);
+ write_txt(file, v->NbSY, v->SY, 5, 1, v->NbTimeStep);
+ write_txt(file, v->NbVY, v->VY, 5, 3, v->NbTimeStep);
+ write_txt(file, v->NbTY, v->TY, 5, 9, v->NbTimeStep);
+}
+
+// Write view to file in MSH format
+
+class Nod{
+ public:
+ int Num;
+ double X, Y, Z;
+ Nod() : Num(0), X(0.), Y(0.), Z(0.) {}
+ Nod(double x, double y, double z) : Num(0), X(x), Y(y), Z(z) {}
+};
+
+class NodCompPos{
+ public:
+ bool operator()(const Nod ent1, const Nod ent2) const
+ {
+ double tol = CTX.lc * 1.e-10 ;
+ if(ent1.X - ent2.X > tol) return true;
+ if(ent1.X - ent2.X < -tol) return false;
+ if(ent1.Y - ent2.Y > tol) return true;
+ if(ent1.Y - ent2.Y < -tol) return false;
+ if(ent1.Z - ent2.Z > tol) return true;
+ return false;
+ }
+};
+
+static void get_nod(int nbelm, List_T *list, int nbnod, int nbcomp,
+ std::set<Nod, NodCompPos> *nodes,
+ int *numelm)
+{
+ if(nbelm){
+ int nb = List_Nbr(list) / nbelm;
+ for(int i = 0; i < List_Nbr(list); i+=nb){
+ double *x = (double *)List_Pointer_Fast(list, i);
+ double *y = (double *)List_Pointer_Fast(list, i + nbnod);
+ double *z = (double *)List_Pointer_Fast(list, i + 2 * nbnod);
+ for(int j = 0; j < nbnod; j++) {
+ Nod n(x[j], y[j], z[j]);
+ std::set<Nod, NodCompPos>::iterator it = nodes->find(n);
+ if(it == nodes->end()){
+ n.Num = nodes->size() + 1;
+ nodes->insert(n);
+ }
+ }
+ (*numelm)++;
+ }
+ }
+}
+
+static void print_elm(FILE *file, int num, int nbnod, Nod nod[8],
+ int nbcomp, double *vals, int dim)
+{
+ // compute average value in elm
+ double d = 0.;
+ for(int k = 0; k < nbnod; k++) {
+ double *v = &vals[nbcomp * k];
+ switch(nbcomp) {
+ case 1: // scalar
+ d += v[0];
+ break;
+ case 3 : // vector
+ d += sqrt(DSQR(v[0]) + DSQR(v[1]) + DSQR(v[2]));
+ break;
+ case 9 : // tensor
+ d += ComputeVonMises(v);
+ break;
+ }
+ }
+ d /= (double)nbnod;
+
+ // assign val as elementary region number
+ int ele = (int)fabs(d), phys = 1;
+
+ switch(dim){
+ case 0:
+ fprintf(file, "%d 15 %d %d 1 %d\n", num, phys, ele, nod[0].Num);
+ break;
+ case 1:
+ fprintf(file, "%d 1 %d %d 2 %d %d\n", num, phys, ele, nod[0].Num, nod[1].Num);
+ break;
+ case 2:
+ if(nbnod == 3)
+ fprintf(file, "%d 2 %d %d 3 %d %d %d\n", num, phys, ele,
+ nod[0].Num, nod[1].Num, nod[2].Num);
+ else
+ fprintf(file, "%d 3 %d %d 4 %d %d %d %d\n", num, phys, ele,
+ nod[0].Num, nod[1].Num, nod[2].Num, nod[3].Num);
+ break;
+ case 3:
+ default:
+ if(nbnod == 4)
+ fprintf(file, "%d 4 %d %d 4 %d %d %d %d\n", num, phys, ele,
+ nod[0].Num, nod[1].Num, nod[2].Num, nod[3].Num);
+ else if(nbnod == 5)
+ fprintf(file, "%d 7 %d %d 5 %d %d %d %d %d\n", num, phys, ele,
+ nod[0].Num, nod[1].Num, nod[2].Num, nod[3].Num, nod[4].Num);
+ else if(nbnod == 6)
+ fprintf(file, "%d 6 %d %d 6 %d %d %d %d %d %d\n", num, phys, ele,
+ nod[0].Num, nod[1].Num, nod[2].Num, nod[3].Num, nod[4].Num,
+ nod[5].Num);
+ else
+ fprintf(file, "%d 5 %d %d 8 %d %d %d %d %d %d %d %d\n", num, phys, ele,
+ nod[0].Num, nod[1].Num, nod[2].Num, nod[3].Num, nod[4].Num,
+ nod[5].Num, nod[6].Num, nod[7].Num);
+ break;
+ }
+}
+
+static void print_elms(FILE *file, int nbelm, List_T *list,
+ int nbnod, int nbcomp, int dim,
+ std::set<Nod, NodCompPos> *nodes,
+ int *numelm)
+{
+ Nod nod[8];
+
+ if(nbelm){
+ int nb = List_Nbr(list) / nbelm;
+ for(int i = 0; i < List_Nbr(list); i+=nb){
+ double *x = (double *)List_Pointer_Fast(list, i);
+ double *y = (double *)List_Pointer_Fast(list, i + nbnod);
+ double *z = (double *)List_Pointer_Fast(list, i + 2 * nbnod);
+ double *v = (double *)List_Pointer_Fast(list, i + 3 * nbnod);
+ for(int j = 0; j < nbnod; j++) {
+ Nod n(x[j], y[j], z[j]);
+ std::set<Nod, NodCompPos>::iterator it = nodes->find(n);
+ if(it == nodes->end()){
+ Msg(GERROR, "Unknown node in element");
+ return;
+ }
+ else{
+ nod[j] = (Nod)(*it);
+ }
+ }
+ (*numelm)++;
+ print_elm(file, *numelm, nbnod, nod, nbcomp, v, dim);
+ }
+ }
+}
+
+void WriteViewMSH(Post_View *v, FILE *file)
+{
+ std::set<Nod, NodCompPos> nodes;
+ int numelm = 0;
+ get_nod(v->NbSP, v->SP, 1, 1, &nodes, &numelm);
+ get_nod(v->NbVP, v->VP, 1, 3, &nodes, &numelm);
+ get_nod(v->NbTP, v->TP, 1, 9, &nodes, &numelm);
+ get_nod(v->NbSL, v->SL, 2, 1, &nodes, &numelm);
+ get_nod(v->NbVL, v->VL, 2, 3, &nodes, &numelm);
+ get_nod(v->NbTL, v->TL, 2, 9, &nodes, &numelm);
+ get_nod(v->NbST, v->ST, 3, 1, &nodes, &numelm);
+ get_nod(v->NbVT, v->VT, 3, 3, &nodes, &numelm);
+ get_nod(v->NbTT, v->TT, 3, 9, &nodes, &numelm);
+ get_nod(v->NbSQ, v->SQ, 4, 1, &nodes, &numelm);
+ get_nod(v->NbVQ, v->VQ, 4, 3, &nodes, &numelm);
+ get_nod(v->NbTQ, v->TQ, 4, 9, &nodes, &numelm);
+ get_nod(v->NbSS, v->SS, 4, 1, &nodes, &numelm);
+ get_nod(v->NbVS, v->VS, 4, 3, &nodes, &numelm);
+ get_nod(v->NbTS, v->TS, 4, 9, &nodes, &numelm);
+ get_nod(v->NbSH, v->SH, 8, 1, &nodes, &numelm);
+ get_nod(v->NbVH, v->VH, 8, 3, &nodes, &numelm);
+ get_nod(v->NbTH, v->TH, 8, 9, &nodes, &numelm);
+ get_nod(v->NbSI, v->SI, 6, 1, &nodes, &numelm);
+ get_nod(v->NbVI, v->VI, 6, 3, &nodes, &numelm);
+ get_nod(v->NbTI, v->TI, 6, 9, &nodes, &numelm);
+ get_nod(v->NbSY, v->SY, 5, 1, &nodes, &numelm);
+ get_nod(v->NbVY, v->VY, 5, 3, &nodes, &numelm);
+ get_nod(v->NbTY, v->TY, 5, 9, &nodes, &numelm);
+
+ fprintf(file, "$NOD\n");
+ fprintf(file, "%d\n", (int)nodes.size());
+ std::set<Nod, NodCompPos>::const_iterator it = nodes.begin();
+ std::set<Nod, NodCompPos>::const_iterator ite = nodes.end();
+ for(; it != ite; ++it){
+ Nod n = (Nod)(*it);
+ fprintf(file, "%d %.16g %.16g %.16g\n", n.Num, n.X, n.Y, n.Z);
+ }
+ fprintf(file, "$ENDNOD\n");
+
+ fprintf(file, "$ELM\n");
+ fprintf(file, "%d\n", numelm);
+ numelm = 0;
+ print_elms(file, v->NbSP, v->SP, 1, 1, 0, &nodes, &numelm);
+ print_elms(file, v->NbVP, v->VP, 1, 3, 0, &nodes, &numelm);
+ print_elms(file, v->NbTP, v->TP, 1, 9, 0, &nodes, &numelm);
+ print_elms(file, v->NbSL, v->SL, 2, 1, 1, &nodes, &numelm);
+ print_elms(file, v->NbVL, v->VL, 2, 3, 1, &nodes, &numelm);
+ print_elms(file, v->NbTL, v->TL, 2, 9, 1, &nodes, &numelm);
+ print_elms(file, v->NbST, v->ST, 3, 1, 2, &nodes, &numelm);
+ print_elms(file, v->NbVT, v->VT, 3, 3, 2, &nodes, &numelm);
+ print_elms(file, v->NbTT, v->TT, 3, 9, 2, &nodes, &numelm);
+ print_elms(file, v->NbSQ, v->SQ, 4, 1, 2, &nodes, &numelm);
+ print_elms(file, v->NbVQ, v->VQ, 4, 3, 2, &nodes, &numelm);
+ print_elms(file, v->NbTQ, v->TQ, 4, 9, 2, &nodes, &numelm);
+ print_elms(file, v->NbSS, v->SS, 4, 1, 3, &nodes, &numelm);
+ print_elms(file, v->NbVS, v->VS, 4, 3, 3, &nodes, &numelm);
+ print_elms(file, v->NbTS, v->TS, 4, 9, 3, &nodes, &numelm);
+ print_elms(file, v->NbSH, v->SH, 8, 1, 3, &nodes, &numelm);
+ print_elms(file, v->NbVH, v->VH, 8, 3, 3, &nodes, &numelm);
+ print_elms(file, v->NbTH, v->TH, 8, 9, 3, &nodes, &numelm);
+ print_elms(file, v->NbSI, v->SI, 6, 1, 3, &nodes, &numelm);
+ print_elms(file, v->NbVI, v->VI, 6, 3, 3, &nodes, &numelm);
+ print_elms(file, v->NbTI, v->TI, 6, 9, 3, &nodes, &numelm);
+ print_elms(file, v->NbSY, v->SY, 5, 1, 3, &nodes, &numelm);
+ print_elms(file, v->NbVY, v->VY, 5, 3, 3, &nodes, &numelm);
+ print_elms(file, v->NbTY, v->TY, 5, 9, 3, &nodes, &numelm);
+ fprintf(file, "$ENDELM\n");
+}
+
+// Write view to file
+
+void WriteView(Post_View *v, char *filename, int format, int append)
+{
+ FILE *file;
+ int binary = (format == 1) ? 1 : 0;
+ int parsed = (format == 2);
+ int stl = (format == 3);
+ int txt = (format == 4);
+ int msh = (format == 5);
+
+ if(filename) {
+ file = fopen(filename, append ? (binary ? "ab" : "a") : (binary ? "wb" : "w"));
+ if(!file) {
+ Msg(GERROR, "Unable to open file '%s'", filename);
+ return;
+ }
+ if(!append)
+ Msg(STATUS2, "Writing '%s'", filename);
+ }
+ else
+ file = stdout;
+
+ if(stl)
+ WriteViewSTL(v, file);
+ else if(txt)
+ WriteViewTXT(v, file);
+ else if(msh)
+ WriteViewMSH(v, file);
+ else
+ WriteViewPOS(v, file, binary, parsed, append);
+
+ if(filename) {
+ fclose(file);
+ Msg(STATUS2, "Wrote '%s'", filename);
+ }
+
+}
diff --git a/Fltk/Bitmaps.h b/Fltk/Bitmaps.h
new file mode 100644
index 0000000000000000000000000000000000000000..b959cb2bc70a4e0a0cf0724111821de18b140e3b
--- /dev/null
+++ b/Fltk/Bitmaps.h
@@ -0,0 +1,92 @@
+#ifndef _BITMAPS_H_
+#define _BITMAPS_H_
+
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+// Standard Gmsh Unix icon
+
+#if !defined(WIN32) && !defined(__APPLE__)
+#define gmsh32x32_width 32
+#define gmsh32x32_height 32
+static char gmsh32x32_bits[] = {
+ 0x00, 0x00, 0x00, 0x00, 0x00, 0x80, 0x01, 0x00, 0x00, 0x40, 0x03, 0x00,
+ 0x00, 0x40, 0x03, 0x00, 0x00, 0x20, 0x07, 0x00, 0x00, 0x20, 0x07, 0x00,
+ 0x00, 0x10, 0x0f, 0x00, 0x00, 0x10, 0x0f, 0x00, 0x00, 0x08, 0x1f, 0x00,
+ 0x00, 0x08, 0x1f, 0x00, 0x00, 0x04, 0x3f, 0x00, 0x00, 0x04, 0x3f, 0x00,
+ 0x00, 0x02, 0x7f, 0x00, 0x00, 0x02, 0x7f, 0x00, 0x00, 0x01, 0xff, 0x00,
+ 0x00, 0x01, 0xff, 0x00, 0x80, 0x00, 0xff, 0x01, 0x80, 0x00, 0xff, 0x01,
+ 0x40, 0x00, 0xff, 0x03, 0x40, 0x00, 0xff, 0x03, 0x20, 0x00, 0xff, 0x07,
+ 0x20, 0x00, 0xff, 0x07, 0x10, 0x00, 0xff, 0x0f, 0x10, 0x00, 0xff, 0x0f,
+ 0x08, 0x00, 0xff, 0x1f, 0x08, 0x00, 0xff, 0x1f, 0x04, 0x40, 0xfd, 0x3f,
+ 0x04, 0xa8, 0xea, 0x3f, 0x02, 0x55, 0x55, 0x7f, 0xa2, 0xaa, 0xaa, 0x7a,
+ 0xff, 0xff, 0xff, 0xff, 0x00, 0x00, 0x00, 0x00 };
+#endif
+
+// 'Abort' bitmap
+// disabled until the mesh thread is back
+/*
+#define abort_width 13
+#define abort_height 13
+static char abort_bits[] = {
+ 0x00,0xe0,0x40,0xe0,0x40,0xe0,0x50,0xe1,0x48,0xe2,0x44,0xe4,0x44,0xe4,0x44,
+ 0xe4,0x04,0xe4,0x04,0xe4,0x08,0xe2,0xf0,0xe1,0x00,0xe0};
+*/
+
+// 'Play button' bitmap
+
+#define start_width 9
+#define start_height 13
+static char start_bits[] = {
+ 0x00,0xfe,0x06,0xfe,0x0a,0xfe,0x12,0xfe,0x22,0xfe,0x42,0xfe,0x82,0xfe,0x42,
+ 0xfe,0x22,0xfe,0x12,0xfe,0x0a,0xfe,0x06,0xfe,0x00,0xfe};
+
+// 'Pause button' bitmap
+
+#define stop_width 9
+#define stop_height 13
+static char stop_bits[] = {
+ 0x00,0xfe,0xee,0xfe,0xaa,0xfe,0xaa,0xfe,0xaa,0xfe,0xaa,0xfe,0xaa,0xfe,0xaa,
+ 0xfe,0xaa,0xfe,0xaa,0xfe,0xaa,0xfe,0xee,0xfe,0x00,0xfe};
+
+// 'Rewind button' bitmap
+
+#define rewind_width 12
+#define rewind_height 13
+static char rewind_bits[] = {
+ 0x00,0xf0,0x0e,0xf6,0x0a,0xf5,0x8a,0xf4,0x4a,0xf4,0x2a,0xf4,0x1a,0xf4,0x2a,
+ 0xf4,0x4a,0xf4,0x8a,0xf4,0x0a,0xf5,0x0e,0xf6,0x00,0xf0};
+
+// 'Orthographic projection' bitmap
+#define ortho_width 13
+#define ortho_height 13
+static unsigned char ortho_bits[] = {
+ 0x00, 0x00, 0xf0, 0x0f, 0x18, 0x0c, 0x14, 0x0a, 0xfe, 0x09, 0x12, 0x09,
+ 0x12, 0x09, 0x12, 0x09, 0xf2, 0x0f, 0x0a, 0x05, 0x06, 0x03, 0xfe, 0x01,
+ 0x00, 0x00 };
+
+// '90 degrees rotation' bitmap
+#define rotate_width 12
+#define rotate_height 13
+static unsigned char rotate_bits[] = {
+ 0x00, 0x00, 0xf0, 0x00, 0x08, 0x01, 0x04, 0x02, 0x02, 0x04, 0x02, 0x04,
+ 0x02, 0x04, 0x02, 0x00, 0x24, 0x00, 0x68, 0x00, 0xf0, 0x00, 0x60, 0x00,
+ 0x20, 0x00 };
+
+#endif
diff --git a/Geo/fourierFace.h b/Geo/fourierFace.h
new file mode 100644
index 0000000000000000000000000000000000000000..0c56386f0ad299ffda90145982ab1524b7eff8d3
--- /dev/null
+++ b/Geo/fourierFace.h
@@ -0,0 +1,95 @@
+#ifndef _FOURIER_MODEL_H_
+#define _FOURIER_MODEL_H_
+
+#include "GModel.h"
+
+#if defined(HAVE_FOURIER_MODEL)
+
+#include "GVertex.h"
+#include "GEdge.h"
+#include "GFace.h"
+#include "Range.h"
+
+class fourierVertex : public GVertex {
+ private:
+ MVertex *_v;
+ public:
+ fourierVertex(GModel *m, MVertex *v) : GVertex(m, v->getNum()), _v(v)
+ {
+ mesh_vertices.push_back(v);
+ }
+ virtual ~fourierVertex() {}
+ virtual GPoint point() const { return GPoint(_v->x(), _v->y(), _v->z(), this); }
+ virtual double x() const { return _v->x(); }
+ virtual double y() const { return _v->y(); }
+ virtual double z() const { return _v->z(); }
+ virtual double prescribedMeshSizeAtVertex() const { return 0.1; }
+ ModelType getNativeType() const { return FourierModel; }
+};
+
+class fourierEdge : public GEdge {
+ public:
+ fourierEdge(GModel *m, int num, GVertex *v1, GVertex *v2);
+ virtual ~fourierEdge() {}
+ double period() const { throw ; }
+ virtual bool periodic(int dim=0) const { return false; }
+ virtual Range<double> parBounds(int i) const { throw; }
+ virtual GeomType geomType() const { throw; }
+ virtual bool degenerate(int) const { return false; }
+ virtual bool continuous(int dim) const { return true; }
+ virtual GPoint point(double p) const { throw; }
+ virtual GPoint closestPoint(const SPoint3 & queryPoint) { throw; }
+ virtual int containsPoint(const SPoint3 &pt) const { throw; }
+ virtual int containsParam(double pt) const { throw; }
+ virtual SVector3 firstDer(double par) const { throw; }
+ virtual SPoint2 reparamOnFace(GFace * face, double epar, int dir) const { throw; }
+ virtual double parFromPoint(const SPoint3 &pt) const { throw; }
+ virtual int minimumMeshSegments () const { throw; }
+ virtual int minimumDrawSegments () const { throw; }
+ ModelType getNativeType() const { return FourierModel; }
+};
+
+class fourierFace : public GFace {
+ private:
+ // flag to know if is the face is already meshed
+ int _discrete;
+ // floag to know if the face is just a plane
+ int _plane;
+ // vertices and edges purely local to the face (not shared with the model)
+ GVertex *_v[4];
+ GEdge *_e[4];
+ public:
+ fourierFace(GModel *m, int num);
+ fourierFace(GFace *f, std::vector<MVertex*> &loop, std::vector<MVertex*> &hole);
+ virtual ~fourierFace();
+
+ void meshBoundary();
+
+ Range<double> parBounds(int i) const;
+ virtual int paramDegeneracies(int dir, double *par) { return 0; }
+
+ virtual GPoint point(double par1, double par2) const;
+ virtual GPoint closestPoint(const SPoint3 & queryPoint) const ;
+
+ virtual int containsPoint(const SPoint3 &pt) const;
+ virtual int containsParam(const SPoint2 &pt) const;
+
+ virtual SVector3 normal(const SPoint2 ¶m) const;
+ virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const {throw;}
+
+ virtual GEntity::GeomType geomType() const;
+ virtual int geomDirection() const { return 1; }
+
+ virtual bool continuous(int dim) const { return true; }
+ virtual bool periodic(int dim) const { return false; }
+ virtual bool degenerate(int dim) const { return false; }
+ virtual double period(int dir) const {throw;}
+ ModelType getNativeType() const { return FourierModel; }
+ void * getNativePtr() const {throw;}
+ virtual bool surfPeriodic(int dim) const {throw;}
+ virtual SPoint2 parFromPoint(const SPoint3 &) const;
+};
+
+#endif
+
+#endif
diff --git a/Plugin/DecomposeInSimplex.cpp b/Plugin/DecomposeInSimplex.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8e6f5e9a69cfeb674c4a56f3fdf8859a989b913e
--- /dev/null
+++ b/Plugin/DecomposeInSimplex.cpp
@@ -0,0 +1,253 @@
+// $Id: DecomposeInSimplex.cpp,v 1.22 2007-05-04 10:45:08 geuzaine Exp $
+//
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "Plugin.h"
+#include "DecomposeInSimplex.h"
+#include "List.h"
+#include "Tree.h"
+#include "Views.h"
+#include "Context.h"
+#include "Malloc.h"
+#include "Numeric.h"
+
+extern Context_T CTX;
+
+StringXNumber DecomposeInSimplexOptions_Number[] = {
+ {GMSH_FULLRC, "iView", NULL, -1.}
+};
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterDecomposeInSimplexPlugin()
+ {
+ return new GMSH_DecomposeInSimplexPlugin();
+ }
+}
+
+GMSH_DecomposeInSimplexPlugin::GMSH_DecomposeInSimplexPlugin()
+{
+ ;
+}
+
+void GMSH_DecomposeInSimplexPlugin::getName(char *name) const
+{
+ strcpy(name, "Decompose in Simplex");
+}
+
+void GMSH_DecomposeInSimplexPlugin::getInfos(char *author, char *copyright,
+ char *help_text) const
+{
+ strcpy(author, "C. Geuzaine");
+ strcpy(copyright, "DGR (www.multiphysics.com)");
+ strcpy(help_text,
+ "Plugin(DecomposeInSimplex) decomposes all\n"
+ "non-simplectic elements (quadrangles, prisms,\n"
+ "hexahedra, pyramids) in the view `iView' into\n"
+ "simplices (triangles, tetrahedra). If `iView' < 0,\n"
+ "the plugin is run on the current view.\n"
+ "\n"
+ "Plugin(DecomposeInSimplex) is executed\n"
+ "in-place.\n");
+}
+
+int GMSH_DecomposeInSimplexPlugin::getNbOptions() const
+{
+ return sizeof(DecomposeInSimplexOptions_Number) / sizeof(StringXNumber);
+}
+
+StringXNumber *GMSH_DecomposeInSimplexPlugin::getOption(int iopt)
+{
+ return &DecomposeInSimplexOptions_Number[iopt];
+}
+
+void GMSH_DecomposeInSimplexPlugin::catchErrorMessage(char *errorMessage) const
+{
+ strcpy(errorMessage, "DecomposeInSimplex failed...");
+}
+
+static void decomposeList(Post_View *v, int nbNod, int nbComp,
+ List_T **listIn, int *nbIn, List_T *listOut, int *nbOut)
+{
+ double xNew[4], yNew[4], zNew[4];
+ double *valNew = new double[v->NbTimeStep * nbComp * nbNod];
+ DecomposeInSimplex dec(nbNod, nbComp, v->NbTimeStep);
+
+ if(!(*nbIn))
+ return;
+
+ v->Changed = 1;
+
+ int nb = List_Nbr(*listIn) / (*nbIn);
+ for(int i = 0; i < List_Nbr(*listIn); i += nb){
+ double *x = (double *)List_Pointer(*listIn, i);
+ double *y = (double *)List_Pointer(*listIn, i + nbNod);
+ double *z = (double *)List_Pointer(*listIn, i + 2 * nbNod);
+ double *val = (double *)List_Pointer(*listIn, i + 3 * nbNod);
+ for(int j = 0; j < dec.numSimplices(); j++){
+ dec.decompose(j, x, y, z, val, xNew, yNew, zNew, valNew);
+ for(int k = 0; k < dec.numSimplexNodes(); k++)
+ List_Add(listOut, &xNew[k]);
+ for(int k = 0; k < dec.numSimplexNodes(); k++)
+ List_Add(listOut, &yNew[k]);
+ for(int k = 0; k < dec.numSimplexNodes(); k++)
+ List_Add(listOut, &zNew[k]);
+ for(int k = 0; k < dec.numSimplexNodes()*v->NbTimeStep*nbComp; k++)
+ List_Add(listOut, &valNew[k]);
+ (*nbOut)++;
+ }
+ }
+
+ delete [] valNew;
+
+ List_Reset(*listIn);
+ *nbIn = 0;
+}
+
+Post_View *GMSH_DecomposeInSimplexPlugin::execute(Post_View * v)
+{
+ int iView = (int)DecomposeInSimplexOptions_Number[0].def;
+
+ if(iView < 0)
+ iView = v ? v->Index : 0;
+
+ if(!List_Pointer_Test(CTX.post.list, iView)) {
+ Msg(GERROR, "View[%d] does not exist", iView);
+ return v;
+ }
+
+ Post_View *v1 = *(Post_View **)List_Pointer(CTX.post.list, iView);
+
+ // Bail out if the view is an alias or if other views duplicate it
+ if(v1->AliasOf || v1->Links) {
+ Msg(GERROR, "DecomposeInSimplex cannot be applied to an aliased view");
+ return 0;
+ }
+
+ // quads
+ decomposeList(v1, 4, 1, &v1->SQ, &v1->NbSQ, v1->ST, &v1->NbST);
+ decomposeList(v1, 4, 3, &v1->VQ, &v1->NbVQ, v1->VT, &v1->NbVT);
+ decomposeList(v1, 4, 9, &v1->TQ, &v1->NbTQ, v1->TT, &v1->NbTT);
+
+ // hexas
+ decomposeList(v1, 8, 1, &v1->SH, &v1->NbSH, v1->SS, &v1->NbSS);
+ decomposeList(v1, 8, 3, &v1->VH, &v1->NbVH, v1->VS, &v1->NbVS);
+ decomposeList(v1, 8, 9, &v1->TH, &v1->NbTH, v1->TS, &v1->NbTS);
+
+ // prisms
+ decomposeList(v1, 6, 1, &v1->SI, &v1->NbSI, v1->SS, &v1->NbSS);
+ decomposeList(v1, 6, 3, &v1->VI, &v1->NbVI, v1->VS, &v1->NbVS);
+ decomposeList(v1, 6, 9, &v1->TI, &v1->NbTI, v1->TS, &v1->NbTS);
+
+ // pyramids
+ decomposeList(v1, 5, 1, &v1->SY, &v1->NbSY, v1->SS, &v1->NbSS);
+ decomposeList(v1, 5, 3, &v1->VY, &v1->NbVY, v1->VS, &v1->NbVS);
+ decomposeList(v1, 5, 9, &v1->TY, &v1->NbTY, v1->TS, &v1->NbTS);
+
+ return v1;
+}
+
+// Utility class
+
+DecomposeInSimplex::DecomposeInSimplex(int numNodes, int numComponents, int numTimeSteps)
+ : _numNodes(numNodes), _numComponents(numComponents), _numTimeSteps(numTimeSteps)
+{
+ ;
+}
+
+int DecomposeInSimplex::numSimplices()
+{
+ switch(_numNodes) {
+ case 4: return 2; // quad -> 2 tris
+ case 5: return 2; // pyramid -> 2 tets
+ case 6: return 3; // prism -> 3 tets
+ case 8: return 6; // hexa -> 6 tets
+ }
+ return 0;
+}
+
+int DecomposeInSimplex::numSimplexNodes()
+{
+ if(_numNodes == 4)
+ return 3; // quad -> tris
+ else
+ return 4; // all others -> tets
+}
+
+void DecomposeInSimplex::reorder(int map[4], int n,
+ double *x, double *y, double *z, double *val,
+ double *xn, double *yn, double *zn, double *valn)
+{
+ for(int i = 0; i < n; i++) {
+ xn[i] = x[map[i]];
+ yn[i] = y[map[i]];
+ zn[i] = z[map[i]];
+ }
+
+ int map2[4] = {map[0], map[1], map[2], map[3]};
+#if 0
+ // make sure to create tets with positive volume?
+ if(n == 4){ // tets
+ double mat[3][3];
+ mat[0][0] = xn[1] - xn[0]; mat[0][1] = xn[2] - xn[0]; mat[0][2] = xn[3] - xn[0];
+ mat[1][0] = yn[1] - yn[0]; mat[1][1] = yn[2] - yn[0]; mat[1][2] = yn[3] - yn[0];
+ mat[2][0] = zn[1] - zn[0]; mat[2][1] = zn[2] - zn[0]; mat[2][2] = zn[3] - zn[0];
+ if(det3x3(mat) < 0.){
+ map2[0] = map[1];
+ map2[1] = map[0];
+ xn[0] = x[map2[0]];
+ yn[0] = y[map2[0]];
+ zn[0] = z[map2[0]];
+ xn[1] = x[map2[1]];
+ yn[1] = y[map2[1]];
+ zn[1] = z[map2[1]];
+ }
+ }
+#endif
+
+ for(int ts = 0; ts < _numTimeSteps; ts++)
+ for(int i = 0; i < n; i++) {
+ for(int j = 0; j < _numComponents; j++)
+ valn[ts*n*_numComponents + i*_numComponents + j] =
+ val[ts*_numNodes*_numComponents + map2[i]*_numComponents + j];
+ }
+}
+
+void DecomposeInSimplex::decompose(int num,
+ double *x, double *y, double *z, double *val,
+ double *xn, double *yn, double *zn, double *valn)
+{
+ int quadTri[2][4] = {{0,1,2,-1}, {0,2,3,-1}};
+ int hexaTet[6][4] = {{0,1,3,7}, {0,4,1,7}, {1,4,5,7}, {1,2,3,7}, {1,6,2,7}, {1,5,6,7}};
+ int prisTet[3][4] = {{0,1,2,4}, {0,2,4,5}, {0,3,4,5}};
+ int pyraTet[2][4] = {{0,1,3,4}, {1,2,3,4}};
+
+ if(num < 0 || num > numSimplices()-1) {
+ Msg(GERROR, "Invalid decomposition");
+ num = 0;
+ }
+
+ switch(_numNodes) {
+ case 4: reorder(quadTri[num], 3, x, y, z, val, xn, yn, zn, valn); break ;
+ case 8: reorder(hexaTet[num], 4, x, y, z, val, xn, yn, zn, valn); break ;
+ case 6: reorder(prisTet[num], 4, x, y, z, val, xn, yn, zn, valn); break ;
+ case 5: reorder(pyraTet[num], 4, x, y, z, val, xn, yn, zn, valn); break ;
+ }
+}
diff --git a/Plugin/DecomposeInSimplex.h b/Plugin/DecomposeInSimplex.h
new file mode 100644
index 0000000000000000000000000000000000000000..881024d2f50fced1ca9835cb9dd35cf31244bfa6
--- /dev/null
+++ b/Plugin/DecomposeInSimplex.h
@@ -0,0 +1,67 @@
+#ifndef _DECOMPOSE_IN_SIMPLEX_H_
+#define _DECOMPOSE_IN_SIMPLEX_H_
+
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "Plugin.h"
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterDecomposeInSimplexPlugin();
+}
+
+class GMSH_DecomposeInSimplexPlugin : public GMSH_Post_Plugin
+{
+public:
+ GMSH_DecomposeInSimplexPlugin();
+ void getName(char *name) const;
+ void getInfos(char *author, char *copyright, char *help_text) const;
+ void catchErrorMessage(char *errorMessage) const;
+ int getNbOptions() const;
+ StringXNumber* getOption(int iopt);
+ Post_View *execute(Post_View *);
+};
+
+class DecomposeInSimplex{
+ private:
+ // how many nodes in the element to decompose
+ int _numNodes;
+ // how many field components
+ int _numComponents;
+ // how many time steps
+ int _numTimeSteps;
+ // create a simplex
+ void reorder(int map[4], int n,
+ double *x, double *y, double *z, double *val,
+ double *xn, double *yn, double *zn, double *valn);
+ public:
+ // default constructor
+ DecomposeInSimplex(int numNodes, int numComponents, int numTimeSteps=1);
+ // the number of simplices into which the element is decomposed
+ int numSimplices();
+ // the number of nodes of the simplex
+ int numSimplexNodes();
+ // returns the i-th simplex in the decomposition
+ void decompose(int num,
+ double *x, double *y, double *z, double *val,
+ double *xn, double *yn, double *zn, double *valn);
+};
+
+#endif
diff --git a/benchmarks/bugs/geom8eu_needs_rand_1e-9.geo b/benchmarks/bugs/geom8eu_needs_rand_1e-9.geo
new file mode 100644
index 0000000000000000000000000000000000000000..3ab453faff7bf4eda3a8b6609a6d5f0810537d62
--- /dev/null
+++ b/benchmarks/bugs/geom8eu_needs_rand_1e-9.geo
@@ -0,0 +1,487 @@
+Point(101) = {0, 0, -6, 20};
+Point(102) = {100, 0, -6, 20};
+Point(103) = {-100, 0, -6, 20};
+Point(104) = {0, 100, -6, 20};
+Point(105) = {0, -100, -6, 20};
+Point(106) = {0, 0, 194, 20};
+Point(107) = {0, 0, -206, 20};
+Point(201) = {-1, -1, -7, 0.3};
+Point(202) = {-1, -1, -5, 0.3};
+Point(203) = {-1, 1, -7, 0.3};
+Point(204) = {-1, 1, -5, 0.3};
+Point(205) = {1, -1, -7, 0.3};
+Point(206) = {1, -1, -5, 0.3};
+Point(207) = {1, 1, -7, 0.3};
+Point(208) = {1, 1, -5, 0.3};
+Point(301) = {0, -1, -4, 0.4};
+Point(302) = {0, -1, 6, 0.4};
+Point(303) = {0, 1, -4, 0.4};
+Point(304) = {0, 1, 6, 0.4};
+Point(305) = {-0.2, -1, -4, 0.4};
+Point(306) = {-0.2, -1, 6, 0.4};
+Point(307) = {-0.2, 1, -4, 0.4};
+Point(308) = {-0.2, 1, 6, 0.4};
+Point(401) = {0, -1, -8, 0.6};
+Point(402) = {0, -1, -18, 0.6};
+Point(403) = {0, 1, -8, 0.6};
+Point(404) = {0, 1, -18, 0.6};
+Point(405) = {-0.2, -1, -8, 0.6};
+Point(406) = {-0.2, -1, -18, 0.6};
+Point(407) = {-0.2, 1, -8, 0.6};
+Point(408) = {-0.2, 1, -18, 0.6};
+Point(1001) = {0, -0.01, 0.9995000000000001, 0.0007};
+Point(1002) = {0, -0.01, 1.0005, 0.0007};
+Point(1003) = {0, 0.01, 0.9995000000000001, 0.0007};
+Point(1004) = {0, 0.01, 1.0005, 0.0007};
+Point(1005) = {0, -0.02, 0.997, 0.00175};
+Point(1006) = {0, -0.02, 1.003, 0.00175};
+Point(1007) = {0, 0.02, 0.997, 0.00175};
+Point(1008) = {0, 0.02, 1.003, 0.00175};
+Point(1009) = {0.005590169943749475, -0.02, 0.997, 0.002625};
+Point(1010) = {0.005590169943749475, -0.02, 1.003, 0.002625};
+Point(1011) = {0.005590169943749475, 0.02, 0.997, 0.002625};
+Point(1012) = {0.005590169943749475, 0.02, 1.003, 0.002625};
+Point(1013) = {0, -0.04, 0.991, 0.004375};
+Point(1014) = {0, -0.04, 1.009, 0.004375};
+Point(1015) = {0, 0.04, 0.991, 0.004375};
+Point(1016) = {0, 0.04, 1.009, 0.004375};
+Point(1017) = {0.01397542485937369, -0.04, 0.991, 0.006562500000000001};
+Point(1018) = {0.01397542485937369, -0.04, 1.009, 0.006562500000000001};
+Point(1019) = {0.01397542485937369, 0.04, 0.991, 0.006562500000000001};
+Point(1020) = {0.01397542485937369, 0.04, 1.009, 0.006562500000000001};
+Point(1021) = {0, -0.08, 0.973, 0.0109375};
+Point(1022) = {0, -0.08, 1.027, 0.0109375};
+Point(1023) = {0, 0.08, 0.973, 0.0109375};
+Point(1024) = {0, 0.08, 1.027, 0.0109375};
+Point(1025) = {0.03493856214843422, -0.08, 0.973, 0.01640625};
+Point(1026) = {0.03493856214843422, -0.08, 1.027, 0.01640625};
+Point(1027) = {0.03493856214843422, 0.08, 0.973, 0.01640625};
+Point(1028) = {0.03493856214843422, 0.08, 1.027, 0.01640625};
+Point(1029) = {0, -0.16, 0.919, 0.02734375};
+Point(1030) = {0, -0.16, 1.081, 0.02734375};
+Point(1031) = {0, 0.16, 0.919, 0.02734375};
+Point(1032) = {0, 0.16, 1.081, 0.02734375};
+Point(1033) = {0.08734640537108554, -0.16, 0.919, 0.04101562500000001};
+Point(1034) = {0.08734640537108554, -0.16, 1.081, 0.04101562500000001};
+Point(1035) = {0.08734640537108554, 0.16, 0.919, 0.04101562500000001};
+Point(1036) = {0.08734640537108554, 0.16, 1.081, 0.04101562500000001};
+Point(1037) = {0, -0.32, 0.7569999999999999, 0.06835937500000001};
+Point(1038) = {0, -0.32, 1.243, 0.06835937500000001};
+Point(1039) = {0, 0.32, 0.7569999999999999, 0.06835937500000001};
+Point(1040) = {0, 0.32, 1.243, 0.06835937500000001};
+Point(1041) = {0.2183660134277138, -0.32, 0.7569999999999999, 0.1025390625};
+Point(1042) = {0.2183660134277138, -0.32, 1.243, 0.1025390625};
+Point(1043) = {0.2183660134277138, 0.32, 0.7569999999999999, 0.1025390625};
+Point(1044) = {0.2183660134277138, 0.32, 1.243, 0.1025390625};
+Point(1045) = {0, -0.64, 0.2709999999999999, 0.1708984375};
+Point(1046) = {0, -0.64, 1.729, 0.1708984375};
+Point(1047) = {0, 0.64, 0.2709999999999999, 0.1708984375};
+Point(1048) = {0, 0.64, 1.729, 0.1708984375};
+Point(1049) = {0.5459150335692846, -0.64, 0.2709999999999999, 0.2563476562500001};
+Point(1050) = {0.5459150335692846, -0.64, 1.729, 0.2563476562500001};
+Point(1051) = {0.5459150335692846, 0.64, 0.2709999999999999, 0.2563476562500001};
+Point(1052) = {0.5459150335692846, 0.64, 1.729, 0.2563476562500001};
+Point(1053) = {-2.5, -2.5, -24, 0.75};
+Point(1054) = {-2.5, -2.5, 12, 0.75};
+Point(1055) = {-2.5, 2.5, -24, 0.75};
+Point(1056) = {-2.5, 2.5, 12, 0.75};
+Point(1057) = {2.5, -2.5, -24, 0.75};
+Point(1058) = {2.5, -2.5, 12, 0.75};
+Point(1059) = {2.5, 2.5, -24, 0.75};
+Point(1060) = {2.5, 2.5, 12, 0.75};
+Point(1061) = {-6.25, -6.25, -33, 1.875};
+Point(1062) = {-6.25, -6.25, 21, 1.875};
+Point(1063) = {-6.25, 6.25, -33, 1.875};
+Point(1064) = {-6.25, 6.25, 21, 1.875};
+Point(1065) = {6.25, -6.25, -33, 1.875};
+Point(1066) = {6.25, -6.25, 21, 1.875};
+Point(1067) = {6.25, 6.25, -33, 1.875};
+Point(1068) = {6.25, 6.25, 21, 1.875};
+Point(1069) = {-15.625, -15.625, -46.5, 4.6875};
+Point(1070) = {-15.625, -15.625, 34.5, 4.6875};
+Point(1071) = {-15.625, 15.625, -46.5, 4.6875};
+Point(1072) = {-15.625, 15.625, 34.5, 4.6875};
+Point(1073) = {15.625, -15.625, -46.5, 4.6875};
+Point(1074) = {15.625, -15.625, 34.5, 4.6875};
+Point(1075) = {15.625, 15.625, -46.5, 4.6875};
+Point(1076) = {15.625, 15.625, 34.5, 4.6875};
+Point(1077) = {-39.0625, -39.0625, -66.75, 11.71875};
+Point(1078) = {-39.0625, -39.0625, 54.75, 11.71875};
+Point(1079) = {-39.0625, 39.0625, -66.75, 11.71875};
+Point(1080) = {-39.0625, 39.0625, 54.75, 11.71875};
+Point(1081) = {39.0625, -39.0625, -66.75, 11.71875};
+Point(1082) = {39.0625, -39.0625, 54.75, 11.71875};
+Point(1083) = {39.0625, 39.0625, -66.75, 11.71875};
+Point(1084) = {39.0625, 39.0625, 54.75, 11.71875};
+Ellipse (1) = {103, 101, 103, 104} Plane{0, 0, 1};
+Ellipse (2) = {106, 101, 104, 104} Plane{0, 0, 1};
+Ellipse (3) = {102, 101, 104, 104} Plane{0, 0, 1};
+Ellipse (4) = {107, 101, 104, 104} Plane{0, 0, 1};
+Ellipse (5) = {105, 101, 106, 106} Plane{0, 0, 1};
+Ellipse (6) = {105, 101, 103, 103} Plane{0, 0, 1};
+Ellipse (7) = {105, 101, 107, 107} Plane{0, 0, 1};
+Ellipse (8) = {105, 101, 102, 102} Plane{0, 0, 1};
+Ellipse (9) = {103, 101, 106, 106} Plane{0, 0, 1};
+Ellipse (10) = {106, 101, 102, 102} Plane{0, 0, 1};
+Ellipse (11) = {102, 101, 107, 107} Plane{0, 0, 1};
+Ellipse (12) = {107, 101, 103, 103} Plane{0, 0, 1};
+Line (10031) = {201, 202};
+Line (10032) = {203, 204};
+Line (10033) = {207, 208};
+Line (10034) = {205, 206};
+Line (10035) = {205, 201};
+Line (10036) = {206, 202};
+Line (10037) = {208, 204};
+Line (10038) = {207, 203};
+Line (10039) = {203, 201};
+Line (10040) = {204, 202};
+Line (10041) = {208, 206};
+Line (10042) = {207, 205};
+Line (20057) = {304, 308};
+Line (20058) = {302, 306};
+Line (20059) = {301, 305};
+Line (20060) = {303, 307};
+Line (20061) = {303, 301};
+Line (20062) = {307, 305};
+Line (20063) = {304, 302};
+Line (20064) = {308, 306};
+Line (20065) = {304, 303};
+Line (20066) = {308, 307};
+Line (20067) = {302, 301};
+Line (20068) = {306, 305};
+Line (30057) = {404, 408};
+Line (30058) = {402, 406};
+Line (30059) = {401, 405};
+Line (30060) = {403, 407};
+Line (30061) = {403, 401};
+Line (30062) = {407, 405};
+Line (30063) = {404, 402};
+Line (30064) = {408, 406};
+Line (30065) = {404, 403};
+Line (30066) = {408, 407};
+Line (30067) = {402, 401};
+Line (30068) = {406, 405};
+Line (30082) = {1004, 1003};
+Line (30083) = {1002, 1001};
+Line (30084) = {1004, 1002};
+Line (30085) = {1003, 1001};
+Line (30088) = {1005, 1006};
+Line (30089) = {1007, 1008};
+Line (30090) = {1009, 1010};
+Line (30091) = {1011, 1012};
+Line (30092) = {1005, 1007};
+Line (30093) = {1006, 1008};
+Line (30094) = {1009, 1011};
+Line (30095) = {1010, 1012};
+Line (30096) = {1005, 1009};
+Line (30097) = {1006, 1010};
+Line (30098) = {1007, 1011};
+Line (30099) = {1008, 1012};
+Line (30114) = {1013, 1014};
+Line (30115) = {1015, 1016};
+Line (30116) = {1017, 1018};
+Line (30117) = {1019, 1020};
+Line (30118) = {1013, 1015};
+Line (30119) = {1014, 1016};
+Line (30120) = {1017, 1019};
+Line (30121) = {1018, 1020};
+Line (30122) = {1013, 1017};
+Line (30123) = {1014, 1018};
+Line (30124) = {1015, 1019};
+Line (30125) = {1016, 1020};
+Line (30140) = {1021, 1022};
+Line (30141) = {1023, 1024};
+Line (30142) = {1025, 1026};
+Line (30143) = {1027, 1028};
+Line (30144) = {1021, 1023};
+Line (30145) = {1022, 1024};
+Line (30146) = {1025, 1027};
+Line (30147) = {1026, 1028};
+Line (30148) = {1021, 1025};
+Line (30149) = {1022, 1026};
+Line (30150) = {1023, 1027};
+Line (30151) = {1024, 1028};
+Line (30166) = {1029, 1030};
+Line (30167) = {1031, 1032};
+Line (30168) = {1033, 1034};
+Line (30169) = {1035, 1036};
+Line (30170) = {1029, 1031};
+Line (30171) = {1030, 1032};
+Line (30172) = {1033, 1035};
+Line (30173) = {1034, 1036};
+Line (30174) = {1029, 1033};
+Line (30175) = {1030, 1034};
+Line (30176) = {1031, 1035};
+Line (30177) = {1032, 1036};
+Line (30192) = {1037, 1038};
+Line (30193) = {1039, 1040};
+Line (30194) = {1041, 1042};
+Line (30195) = {1043, 1044};
+Line (30196) = {1037, 1039};
+Line (30197) = {1038, 1040};
+Line (30198) = {1041, 1043};
+Line (30199) = {1042, 1044};
+Line (30200) = {1037, 1041};
+Line (30201) = {1038, 1042};
+Line (30202) = {1039, 1043};
+Line (30203) = {1040, 1044};
+Line (30218) = {1045, 1046};
+Line (30219) = {1047, 1048};
+Line (30220) = {1049, 1050};
+Line (30221) = {1051, 1052};
+Line (30222) = {1045, 1047};
+Line (30223) = {1046, 1048};
+Line (30224) = {1049, 1051};
+Line (30225) = {1050, 1052};
+Line (30226) = {1045, 1049};
+Line (30227) = {1046, 1050};
+Line (30228) = {1047, 1051};
+Line (30229) = {1048, 1052};
+Line (30244) = {1053, 1054};
+Line (30245) = {1055, 1056};
+Line (30246) = {1057, 1058};
+Line (30247) = {1059, 1060};
+Line (30248) = {1053, 1055};
+Line (30249) = {1054, 1056};
+Line (30250) = {1057, 1059};
+Line (30251) = {1058, 1060};
+Line (30252) = {1053, 1057};
+Line (30253) = {1054, 1058};
+Line (30254) = {1055, 1059};
+Line (30255) = {1056, 1060};
+Line (30269) = {1061, 1062};
+Line (30270) = {1063, 1064};
+Line (30271) = {1065, 1066};
+Line (30272) = {1067, 1068};
+Line (30273) = {1061, 1063};
+Line (30274) = {1062, 1064};
+Line (30275) = {1065, 1067};
+Line (30276) = {1066, 1068};
+Line (30277) = {1061, 1065};
+Line (30278) = {1062, 1066};
+Line (30279) = {1063, 1067};
+Line (30280) = {1064, 1068};
+Line (30294) = {1069, 1070};
+Line (30295) = {1071, 1072};
+Line (30296) = {1073, 1074};
+Line (30297) = {1075, 1076};
+Line (30298) = {1069, 1071};
+Line (30299) = {1070, 1072};
+Line (30300) = {1073, 1075};
+Line (30301) = {1074, 1076};
+Line (30302) = {1069, 1073};
+Line (30303) = {1070, 1074};
+Line (30304) = {1071, 1075};
+Line (30305) = {1072, 1076};
+Line (30319) = {1077, 1078};
+Line (30320) = {1079, 1080};
+Line (30321) = {1081, 1082};
+Line (30322) = {1083, 1084};
+Line (30323) = {1077, 1079};
+Line (30324) = {1078, 1080};
+Line (30325) = {1081, 1083};
+Line (30326) = {1082, 1084};
+Line (30327) = {1077, 1081};
+Line (30328) = {1078, 1082};
+Line (30329) = {1079, 1083};
+Line (30330) = {1080, 1084};
+Line Loop (1000014) = {2, -1, 9};
+Ruled Surface (14) = {1000014};
+Line Loop (1000016) = {4, -1, -12};
+Ruled Surface (16) = {1000016};
+Line Loop (1000018) = {2, -3, -10};
+Ruled Surface (18) = {1000018};
+Line Loop (1000020) = {3, -4, -11};
+Ruled Surface (20) = {1000020};
+Line Loop (1000022) = {5, -9, -6};
+Ruled Surface (22) = {1000022};
+Line Loop (1000024) = {7, 12, -6};
+Ruled Surface (24) = {1000024};
+Line Loop (1000026) = {7, -11, -8};
+Ruled Surface (26) = {1000026};
+Line Loop (1000028) = {8, -10, -5};
+Ruled Surface (28) = {1000028};
+Line Loop (1010044) = {10031, -10040, -10032, 10039};
+Plane Surface (10044) = {1010044};
+Line Loop (1010046) = {10042, 10034, -10041, -10033};
+Plane Surface (10046) = {1010046};
+Line Loop (1010048) = {10035, -10039, -10038, 10042};
+Plane Surface (10048) = {1010048};
+Line Loop (1010050) = {10036, -10040, -10037, 10041};
+Plane Surface (10050) = {1010050};
+Line Loop (1010052) = {10034, 10036, -10031, -10035};
+Plane Surface (10052) = {1010052};
+Line Loop (1010054) = {10032, -10037, -10033, 10038};
+Plane Surface (10054) = {1010054};
+Line Loop (1020070) = {20063, 20058, -20064, -20057};
+Plane Surface (20070) = {1020070};
+Line Loop (1020072) = {20061, 20059, -20062, -20060};
+Plane Surface (20072) = {1020072};
+Line Loop (1020074) = {20065, 20060, -20066, -20057};
+Plane Surface (20074) = {1020074};
+Line Loop (1020076) = {20067, 20059, -20068, -20058};
+Plane Surface (20076) = {1020076};
+Line Loop (1020078) = {20064, 20068, -20062, -20066};
+Plane Surface (20078) = {1020078};
+Line Loop (1030070) = {30063, 30058, -30064, -30057};
+Plane Surface (30070) = {1030070};
+Line Loop (1030072) = {30061, 30059, -30062, -30060};
+Plane Surface (30072) = {1030072};
+Line Loop (1030074) = {30065, 30060, -30066, -30057};
+Plane Surface (30074) = {1030074};
+Line Loop (1030076) = {30067, 30059, -30068, -30058};
+Plane Surface (30076) = {1030076};
+Line Loop (1030078) = {30064, 30068, -30062, -30066};
+Plane Surface (30078) = {1030078};
+Line Loop (1030080) = {30063, 30067, -30061, -30065};
+Plane Surface (30080) = {1030080};
+Line Loop (1030087) = {30084, 30083, -30085, -30082};
+Plane Surface (30087) = {1030087};
+Line Loop (1030090) = {20063, 20067, -20061, -20065, 30222, 30219, -30223, -30218};
+Plane Surface (30090) = {1030090};
+Line Loop (1030106) = {30088, 30093, -30089, -30092, 30082, 30085, -30083, -30084};
+Plane Surface (30106) = {1030106};
+Line Loop (1030107) = {-30090, 30094, 30091, -30095};
+Plane Surface (30107) = {1030107};
+Line Loop (1030108) = {-30088, 30096, 30090, -30097};
+Plane Surface (30108) = {1030108};
+Line Loop (1030109) = {30089, 30099, -30091, -30098};
+Plane Surface (30109) = {1030109};
+Line Loop (1030110) = {30092, 30098, -30094, -30096};
+Plane Surface (30110) = {1030110};
+Line Loop (1030111) = {-30093, 30097, 30095, -30099};
+Plane Surface (30111) = {1030111};
+Line Loop (1030132) = {30114, 30119, -30115, -30118, 30092, 30089, -30093, -30088};
+Plane Surface (30132) = {1030132};
+Line Loop (1030133) = {-30116, 30120, 30117, -30121};
+Plane Surface (30133) = {1030133};
+Line Loop (1030134) = {-30114, 30122, 30116, -30123};
+Plane Surface (30134) = {1030134};
+Line Loop (1030135) = {30115, 30125, -30117, -30124};
+Plane Surface (30135) = {1030135};
+Line Loop (1030136) = {30118, 30124, -30120, -30122};
+Plane Surface (30136) = {1030136};
+Line Loop (1030137) = {-30119, 30123, 30121, -30125};
+Plane Surface (30137) = {1030137};
+Line Loop (1030158) = {30140, 30145, -30141, -30144, 30118, 30115, -30119, -30114};
+Plane Surface (30158) = {1030158};
+Line Loop (1030159) = {-30142, 30146, 30143, -30147};
+Plane Surface (30159) = {1030159};
+Line Loop (1030160) = {-30140, 30148, 30142, -30149};
+Plane Surface (30160) = {1030160};
+Line Loop (1030161) = {30141, 30151, -30143, -30150};
+Plane Surface (30161) = {1030161};
+Line Loop (1030162) = {30144, 30150, -30146, -30148};
+Plane Surface (30162) = {1030162};
+Line Loop (1030163) = {-30145, 30149, 30147, -30151};
+Plane Surface (30163) = {1030163};
+Line Loop (1030184) = {30166, 30171, -30167, -30170, 30144, 30141, -30145, -30140};
+Plane Surface (30184) = {1030184};
+Line Loop (1030185) = {-30168, 30172, 30169, -30173};
+Plane Surface (30185) = {1030185};
+Line Loop (1030186) = {-30166, 30174, 30168, -30175};
+Plane Surface (30186) = {1030186};
+Line Loop (1030187) = {30167, 30177, -30169, -30176};
+Plane Surface (30187) = {1030187};
+Line Loop (1030188) = {30170, 30176, -30172, -30174};
+Plane Surface (30188) = {1030188};
+Line Loop (1030189) = {-30171, 30175, 30173, -30177};
+Plane Surface (30189) = {1030189};
+Line Loop (1030210) = {30192, 30197, -30193, -30196, 30170, 30167, -30171, -30166};
+Plane Surface (30210) = {1030210};
+Line Loop (1030211) = {-30194, 30198, 30195, -30199};
+Plane Surface (30211) = {1030211};
+Line Loop (1030212) = {-30192, 30200, 30194, -30201};
+Plane Surface (30212) = {1030212};
+Line Loop (1030213) = {30193, 30203, -30195, -30202};
+Plane Surface (30213) = {1030213};
+Line Loop (1030214) = {30196, 30202, -30198, -30200};
+Plane Surface (30214) = {1030214};
+Line Loop (1030215) = {-30197, 30201, 30199, -30203};
+Plane Surface (30215) = {1030215};
+Line Loop (1030236) = {30218, 30223, -30219, -30222, 30196, 30193, -30197, -30192};
+Plane Surface (30236) = {1030236};
+Line Loop (1030237) = {-30220, 30224, 30221, -30225};
+Plane Surface (30237) = {1030237};
+Line Loop (1030238) = {-30218, 30226, 30220, -30227};
+Plane Surface (30238) = {1030238};
+Line Loop (1030239) = {30219, 30229, -30221, -30228};
+Plane Surface (30239) = {1030239};
+Line Loop (1030240) = {30222, 30228, -30224, -30226};
+Plane Surface (30240) = {1030240};
+Line Loop (1030241) = {-30223, 30227, 30225, -30229};
+Plane Surface (30241) = {1030241};
+Line Loop (1030257) = {30244, 30249, -30245, -30248};
+Plane Surface (30257) = {1030257};
+Line Loop (1030259) = {-30246, 30250, 30247, -30251};
+Plane Surface (30259) = {1030259};
+Line Loop (1030261) = {-30244, 30252, 30246, -30253};
+Plane Surface (30261) = {1030261};
+Line Loop (1030263) = {30245, 30255, -30247, -30254};
+Plane Surface (30263) = {1030263};
+Line Loop (1030265) = {30248, 30254, -30250, -30252};
+Plane Surface (30265) = {1030265};
+Line Loop (1030267) = {-30249, 30253, 30251, -30255};
+Plane Surface (30267) = {1030267};
+Line Loop (1030282) = {30269, 30274, -30270, -30273};
+Plane Surface (30282) = {1030282};
+Line Loop (1030284) = {-30271, 30275, 30272, -30276};
+Plane Surface (30284) = {1030284};
+Line Loop (1030286) = {-30269, 30277, 30271, -30278};
+Plane Surface (30286) = {1030286};
+Line Loop (1030288) = {30270, 30280, -30272, -30279};
+Plane Surface (30288) = {1030288};
+Line Loop (1030290) = {30273, 30279, -30275, -30277};
+Plane Surface (30290) = {1030290};
+Line Loop (1030292) = {-30274, 30278, 30276, -30280};
+Plane Surface (30292) = {1030292};
+Line Loop (1030307) = {30294, 30299, -30295, -30298};
+Plane Surface (30307) = {1030307};
+Line Loop (1030309) = {-30296, 30300, 30297, -30301};
+Plane Surface (30309) = {1030309};
+Line Loop (1030311) = {-30294, 30302, 30296, -30303};
+Plane Surface (30311) = {1030311};
+Line Loop (1030313) = {30295, 30305, -30297, -30304};
+Plane Surface (30313) = {1030313};
+Line Loop (1030315) = {30298, 30304, -30300, -30302};
+Plane Surface (30315) = {1030315};
+Line Loop (1030317) = {-30299, 30303, 30301, -30305};
+Plane Surface (30317) = {1030317};
+Line Loop (1030332) = {30319, 30324, -30320, -30323};
+Plane Surface (30332) = {1030332};
+Line Loop (1030334) = {-30321, 30325, 30322, -30326};
+Plane Surface (30334) = {1030334};
+Line Loop (1030336) = {-30319, 30327, 30321, -30328};
+Plane Surface (30336) = {1030336};
+Line Loop (1030338) = {30320, 30330, -30322, -30329};
+Plane Surface (30338) = {1030338};
+Line Loop (1030340) = {30323, 30329, -30325, -30327};
+Plane Surface (30340) = {1030340};
+Line Loop (1030342) = {-30324, 30328, 30326, -30330};
+Plane Surface (30342) = {1030342};
+Surface Loop (1000001) = {30282, 30284, 30286, 30288, 30290, 30292, 30257, 30259, 30261, 30263, 30265, 30267};
+Volume (1) = {1000001};
+Surface Loop (1000002) = {30307, 30309, 30311, 30313, 30315, 30317, 30282, 30284, 30286, 30288, 30290, 30292};
+Volume (2) = {1000002};
+Surface Loop (1000003) = {30332, 30334, 30336, 30338, 30340, 30342, 30307, 30309, 30311, 30313, 30315, 30317};
+Volume (3) = {1000003};
+Surface Loop (1000004) = {22, 28, 26, 24, 16, 20, 18, 14, 30332, 30334, 30336, 30338, 30340, 30342};
+Volume (4) = {1000004};
+Surface Loop (1001000) = {30257, 30259, 30261, 30263, 30265, 30267, 10052, 10046, 10048, 10044, 10050, 10054, 30090, 30237, 30238, 30239, 30240, 30241, 20070, 20076, 20072, 20078, 20074, 30080, 30070, 30076, 30072, 30078, 30074};
+Volume (1000) = {1001000};
+Surface Loop (1030113) = {30106, 30107, 30108, 30109, 30110, 30111, 30087};
+Volume (30113) = {1030113};
+Surface Loop (1030139) = {30132, 30133, 30134, 30135, 30136, 30137, 30107, 30108, 30109, 30110, 30111};
+Volume (30139) = {1030139};
+Surface Loop (1030165) = {30158, 30159, 30160, 30161, 30162, 30163, 30133, 30134, 30135, 30136, 30137};
+Volume (30165) = {1030165};
+Surface Loop (1030191) = {30184, 30185, 30186, 30187, 30188, 30189, 30159, 30160, 30161, 30162, 30163};
+Volume (30191) = {1030191};
+Surface Loop (1030217) = {30210, 30211, 30212, 30213, 30214, 30215, 30185, 30186, 30187, 30188, 30189};
+Volume (30217) = {1030217};
+Surface Loop (1030243) = {30236, 30237, 30238, 30239, 30240, 30241, 30211, 30212, 30213, 30214, 30215};
+Volume (30243) = {1030243};
+Physical Surface (1) = {14, 16, 18, 20, 22, 24, 26, 28};
+Physical Surface (2) = {10044, 10046, 10048, 10050, 10052, 10054};
+Physical Surface (3) = {30087};
+Physical Surface (5) = {30070, 30072, 30074, 30076, 30078, 30080};
diff --git a/examples/cube.geo b/examples/cube.geo
deleted file mode 100644
index decdd4a696bb42dae6ae5bda42e91f3bb5b099e1..0000000000000000000000000000000000000000
--- a/examples/cube.geo
+++ /dev/null
@@ -1,83 +0,0 @@
-lcar1 = .0275;
-
-length = 0.25;
-height = 1.0;
-depth = 0.25;
-
-Point(newp) = {length/2,height/2,depth,lcar1}; /* Point 1 */
-Point(newp) = {length/2,height/2,0,lcar1}; /* Point 2 */
-Point(newp) = {-length/2,height/2,depth,lcar1}; /* Point 3 */
-Point(newp) = {-length/2,-height/2,depth,lcar1}; /* Point 4 */
-Point(newp) = {length/2,-height/2,depth,lcar1}; /* Point 5 */
-Point(newp) = {length/2,-height/2,0,lcar1}; /* Point 6 */
-Point(newp) = {-length/2,height/2,0,lcar1}; /* Point 7 */
-Point(newp) = {-length/2,-height/2,0,lcar1}; /* Point 8 */
-Line(1) = {3,1};
-Line(2) = {3,7};
-Line(3) = {7,2};
-Line(4) = {2,1};
-Line(5) = {1,5};
-Line(6) = {5,4};
-Line(7) = {4,8};
-Line(8) = {8,6};
-Line(9) = {6,5};
-Line(10) = {6,2};
-Line(11) = {3,4};
-Line(12) = {8,7};
-Line Loop(13) = {-6,-5,-1,11};
-Plane Surface(14) = {13};
-Line Loop(15) = {4,5,-9,10};
-Plane Surface(16) = {15};
-Line Loop(17) = {-3,-12,8,10};
-Plane Surface(18) = {17};
-Line Loop(19) = {7,12,-2,11};
-Plane Surface(20) = {19};
-Line Loop(21) = {-4,-3,-2,1};
-Plane Surface(22) = {21};
-Line Loop(23) = {8,9,6,7};
-Plane Surface(24) = {23};
-
-Surface Loop(25) = {14,24,-18,22,16,-20};
-Complex Volume(26) = {25};
-
-n = 21;
-
-/*
-sizex = 4;
-sizey = 6;
-sizez = 4;
-*/
-sizex = n*length;
-sizey = n*height;
-sizez = n*depth;
-
-sizex = 9;
-sizey = 33;
-sizez = 9;
-
-
-Transfinite Line {2,4,9,7} = sizez;
-Transfinite Line {11,5,10,12} = sizey;
-Transfinite Line {3,1,6,8} = sizex;
-
-Transfinite Surface {14} = {5,4,3,1};
-Transfinite Surface {16} = {6,2,1,5};
-Transfinite Surface {18} = {6,2,7,8};
-Transfinite Surface {20} = {3,7,8,4};
-Transfinite Surface {22} = {3,7,2,1};
-Transfinite Surface {24} = {4,8,6,5};
-
-Recombine Surface {14,16,18,20,22,24};
-
-Transfinite Volume {26} = {3,7,2,1,4,8,6,5};
-
-Physical Surface(200) = {14,16,18,20,24};
-
-Physical Volume(100) = {26};
-/*
-Mesh.use_cut_plane = 1;
-Mesh.cut_planea = 0;
-Mesh.cut_planeb = 0;
-Mesh.cut_planec = 1;
-Mesh.cut_planed = -0.125;
-*/
\ No newline at end of file