Skip to content
Snippets Groups Projects
Commit e974a475 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

This commit was manufactured by cvs2svn to create tag 'gmsh_2_0_8'.

parent bce72800
No related branches found
No related tags found
No related merge requests found
//
// 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 > ();
}
#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
// $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);
}
}
}
#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
// $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;
}
#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
This diff is collapsed.
#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
This diff is collapsed.
#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
#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 &param) const;
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) 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
// $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 ;
}
}
#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
This diff is collapsed.
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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment