From f55145d23a42a9d63b7eccc5936b6b4b2292a36c Mon Sep 17 00:00:00 2001 From: Christophe Geuzaine <cgeuzaine@ulg.ac.be> Date: Sun, 27 Nov 2005 01:20:55 +0000 Subject: [PATCH] forgot to indent these two --- Common/AdaptiveViews.cpp | 1675 ++++++++++++++++++++------------------ Common/AdaptiveViews.h | 216 ++--- 2 files changed, 987 insertions(+), 904 deletions(-) diff --git a/Common/AdaptiveViews.cpp b/Common/AdaptiveViews.cpp index fd2743c33d..25b1457a92 100644 --- a/Common/AdaptiveViews.cpp +++ b/Common/AdaptiveViews.cpp @@ -18,9 +18,6 @@ // // Please report all bugs and problems to <gmsh@geuz.org>. -// Don't indent this file -// *INDENT-OFF* - #include <stdio.h> #include <math.h> #include <list> @@ -30,14 +27,15 @@ // A recursive effective implementation -void computeShapeFunctions ( Double_Matrix *coeffs, Double_Matrix *eexps , double u, double v, double w,double *sf); +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; +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; @@ -45,715 +43,787 @@ 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 *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); - } + 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); + return (adapt_point *) & (*it); } -void adapt_edge::Create (int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps) +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) ; + 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) +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) ; + 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) +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) ; + 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) +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) ; + 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) +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) ; + 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) +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) + 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; + 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) +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) + 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,p13,p12); - Recur_Create (t1, maxlevel,level,coeffs,eexps); - adapt_triangle *t2 = new adapt_triangle (p12,p23,p2); - Recur_Create (t2, maxlevel,level,coeffs,eexps); - adapt_triangle *t3 = new adapt_triangle (p23,p13,p3); - Recur_Create (t3, maxlevel,level,coeffs,eexps); - adapt_triangle *t4 = new adapt_triangle (p12,p13,p23); - Recur_Create (t4, maxlevel,level,coeffs,eexps); - t->e[0]=t1;t->e[1]=t2;t->e[2]=t3;t->e[3]=t4; - + 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, p13, p12); + Recur_Create(t1, maxlevel, level, coeffs, eexps); + adapt_triangle *t2 = new adapt_triangle(p12, p23, p2); + Recur_Create(t2, maxlevel, level, coeffs, eexps); + adapt_triangle *t3 = new adapt_triangle(p23, p13, p3); + Recur_Create(t3, maxlevel, level, coeffs, eexps); + adapt_triangle *t4 = new adapt_triangle(p12, p13, p23); + 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) +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) + if(level++ >= maxlevel) return; /* - - p2 p23 p3 - - - p12 pc p34 - - - p1 p14 p4 - - - */ - - 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 ( (p4->x+p3->x)*0.5,(p4->y+p3->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 (p12,p2,p23,pc); - Recur_Create (q2, maxlevel,level,coeffs,eexps); - adapt_quad *q3 = new adapt_quad (pc,p23,p3,p34); - Recur_Create (q3, maxlevel,level,coeffs,eexps); - adapt_quad *q4 = new adapt_quad (p14,pc,p34,p4); - Recur_Create (q4, maxlevel,level,coeffs,eexps); - q->e[0]=q1;q->e[1]=q2;q->e[2]=q3;q->e[3]=q4; - + p2 p23 p3 + + p12 pc p34 + + p1 p14 p4 + */ + + 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((p4->x + p3->x) * 0.5, (p4->y + p3->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(p12, p2, p23, pc); + Recur_Create(q2, maxlevel, level, coeffs, eexps); + adapt_quad *q3 = new adapt_quad(pc, p23, p3, p34); + Recur_Create(q3, maxlevel, level, coeffs, eexps); + adapt_quad *q4 = new adapt_quad(p14, pc, p34, p4); + 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) +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) + 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; + 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) +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) + 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; + 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 ) +void adapt_edge::Error(double AVG, double tol) { adapt_edge *e = *all_elems.begin(); - Recur_Error (e,AVG,tol); + Recur_Error(e, AVG, tol); } -void adapt_triangle::Error ( double AVG , double tol ) +void adapt_triangle::Error(double AVG, double tol) { adapt_triangle *t = *all_elems.begin(); - Recur_Error (t,AVG,tol); + Recur_Error(t, AVG, tol); } -void adapt_quad::Error ( double AVG , double tol ) +void adapt_quad::Error(double AVG, double tol) { adapt_quad *q = *all_elems.begin(); - Recur_Error (q,AVG,tol); + Recur_Error(q, AVG, tol); } -void adapt_tet::Error ( double AVG , double tol ) +void adapt_tet::Error(double AVG, double tol) { adapt_tet *t = *all_elems.begin(); - Recur_Error (t,AVG,tol); + Recur_Error(t, AVG, tol); } -void adapt_hex::Error ( double AVG , double tol ) +void adapt_hex::Error(double AVG, double tol) { adapt_hex *h = *all_elems.begin(); - Recur_Error (h,AVG,tol); + Recur_Error(h, AVG, tol); } -void adapt_edge::Recur_Error ( adapt_edge *e, double AVG, double 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; - } + 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 - { - 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; - } + e->visible = true; } + } } -void adapt_triangle::Recur_Error ( adapt_triangle *t, double AVG, double tol ) +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; - } + 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 - { - 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; - } + 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 ) +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; - } + 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 - { - 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; - } + 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 ) +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; - } + 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 ) +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; + 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; +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) +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; - 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(); + 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_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 ); - } + Double_Vector val(N), res(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); - _Geometry->mult(xyz,XYZ); + for(int k = 0; k < N; ++k) { + val(k) = (*_STval) (ielem, k); + } + _Interpolate->mult(val, res); + _Geometry->mult(xyz, XYZ); double c1 = Cpu(); - int kk=0; - for ( ; it !=ite ; ++it) - { - adapt_point *p = (adapt_point*) &(*it); - p->val = res(kk); - p->X = XYZ(kk,0); - p->Y = XYZ(kk,1); - p->Z = XYZ(kk,2); - if (min > p->val) min = p->val; - if (max < p->val) max = p->val; - kk++; - } + int kk = 0; + for(; it != ite; ++it) { + adapt_point *p = (adapt_point *) & (*it); + p->val = res(kk); + p->X = XYZ(kk, 0); + p->Y = XYZ(kk, 1); + p->Z = XYZ(kk, 2); + if(min > p->val) + min = p->val; + if(max < p->val) + max = 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(); + 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; - } + for(; itt != itte; itt++) { + (*itt)->visible = false; + } - if (!plug || tol != 0.0) - { - ELEM::Error ( max-min, tol ); + if(!plug || tol != 0.0) { + ELEM::Error(max - min, tol); } - if (plug) - plug->assign_specific_visibility (); + if(plug) + plug->assign_specific_visibility(); double c3 = Cpu(); - itt = ELEM::all_elems.begin(); + itt = ELEM::all_elems.begin(); - for ( ;itt != itte ; itt++) - { - if ((*itt)->visible && !(*itt)->e[0] && level != levelmax )return 0; - } + for(; itt != itte; itt++) { + if((*itt)->visible && !(*itt)->e[0] && level != levelmax) + return 0; + } - itt = ELEM::all_elems.begin(); + 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 ); - for (int k=0;k<nbNod;++k)List_Add ( theList , &p[k]->val ); - (*counter)++; - } + 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); + 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; - + t0 += c1 - c0; + t1 += c2 - c1; + t2 += c3 - c2; + t3 += c4 - c3; + return 1; } @@ -765,258 +835,263 @@ int Adaptive_Post_View:: zoomElement (Post_View * view , 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) +void Adaptive_Post_View::setAdaptiveResolutionLevel(Post_View * view, + int level, + GMSH_Post_Plugin * plug) { - if (presentTol==tol && presentZoomLevel == level && !plug)return; + 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); + 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 *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 = 5; - List_Delete(view->SL); - view->NbSL = 0; - view->SL =List_Create ( nbelm * 8, 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->NbSH) - { - TYP = 4; - List_Delete(view->SH); - view->NbSH = 0; - view->SH =List_Create ( nbelm * 4, nbelm , sizeof(double)); - } - - view->NbTimeStep=1; - - + if(view->NbSL) { + TYP = 5; + List_Delete(view->SL); + view->NbSL = 0; + view->SL = List_Create(nbelm * 8, 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->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 == 5)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 == 2)setAdaptiveResolutionLevel_TEMPL <adapt_quad> ( view,level_act,level, plug,&(view->SQ),&(view->NbSQ),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 ++; - } + + while(1) { + if(TYP == 5) + 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 == 2) + setAdaptiveResolutionLevel_TEMPL < adapt_quad > (view, level_act, level, + plug, &(view->SQ), + &(view->NbSQ), 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++; + } view->Changed = 1; presentZoomLevel = level; - presentTol = tol; + presentTol = tol; - printf("finished %g %g %g %g\n",t0,t1,t2,t3); - delete [] done; + 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) +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) + 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) + 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++; + _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); - } + 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) +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]; - } + 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) +void Adaptive_Post_View::initWithLowResolution(Post_View * view) { - List_T *myList; int nbelm; int nbnod; - if (view->NbSL) - { - myList = view->SL; - nbelm = view->NbSL; - nbnod = 2; - } - else if (view->NbST) - { - myList = view->ST; - nbelm = view->NbST; - nbnod = 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(view->NbSL) { + myList = view->SL; + nbelm = view->NbSL; + nbnod = 2; + } + else if(view->NbST) { + myList = view->ST; + nbelm = view->NbST; + nbnod = 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; min = VAL_INF; max = -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 ); + _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); /// 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); - for (int j=0;j<nb-3*nbnod;j++){ - (*_STval)(k,j)=val[j]; - } - k++; + 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); + for(int j = 0; j < nb - 3 * nbnod; j++) { + (*_STval) (k, j) = val[j]; } - setAdaptiveResolutionLevel(view,0); + k++; + } + setAdaptiveResolutionLevel(view, 0); } -Adaptive_Post_View:: Adaptive_Post_View (Post_View *view, List_T *_c , List_T *_pol) - : tol(1.e-3) +Adaptive_Post_View::Adaptive_Post_View(Post_View * view, List_T * _c, + List_T * _pol) + :tol(1.e-3) { _Interpolate = _Geometry = 0; - _coefs = new Double_Matrix ( List_Nbr (_c) , List_Nbr (_c) ); - _eexps = new Double_Matrix ( List_Nbr (_c) , 3 ); - - 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; - } + _coefs = new Double_Matrix(List_Nbr(_c), List_Nbr(_c)); + _eexps = new Double_Matrix(List_Nbr(_c), 3); + + 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; } - initWithLowResolution (view); + } + initWithLowResolution(view); } Adaptive_Post_View::~Adaptive_Post_View() @@ -1027,11 +1102,13 @@ Adaptive_Post_View::~Adaptive_Post_View() delete _STposY; delete _STposZ; delete _STval; - if(_Interpolate)delete _Interpolate; - if(_Geometry)delete _Geometry; - cleanElement<adapt_edge>(); - cleanElement<adapt_triangle> (); - cleanElement<adapt_tet> (); - cleanElement<adapt_hex> (); - cleanElement<adapt_quad> (); + if(_Interpolate) + delete _Interpolate; + if(_Geometry) + delete _Geometry; + cleanElement < adapt_edge > (); + cleanElement < adapt_triangle > (); + cleanElement < adapt_tet > (); + cleanElement < adapt_hex > (); + cleanElement < adapt_quad > (); } diff --git a/Common/AdaptiveViews.h b/Common/AdaptiveViews.h index de6f8bbcaf..49970d7662 100644 --- a/Common/AdaptiveViews.h +++ b/Common/AdaptiveViews.h @@ -29,29 +29,30 @@ 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... +// 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; + double x, y, z; + double X, Y, Z, val; double shape_functions[128]; - static adapt_point * New ( double x, double y, double z, Double_Matrix *coeffs, Double_Matrix *eexps); - void print ()const + 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); + printf("p %g %g\n" ,x,y); } - bool operator < ( const adapt_point & other ) const + 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; + 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; @@ -60,13 +61,13 @@ public : 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; - } + 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.; @@ -78,12 +79,13 @@ public: } void print () { - printf ("p1 %g %g p2 %g %g \n",p[0]->x,p[0]->y,p[1]->x,p[1]->y); + 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 ); + 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]; @@ -91,20 +93,17 @@ public: static int nbNod; }; - - class adapt_triangle { public: - adapt_triangle (adapt_point *p1,adapt_point *p2,adapt_point *p3) + 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; + e[0] = e[1] = e[2] = e[3] = 0; } - inline double V () const { return (p[0]->val + p[1]->val + p[2]->val)/3.; @@ -117,15 +116,17 @@ public: } 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); + 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 ); + 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]; + adapt_point *p[3]; + adapt_triangle *e[4]; static std::list<adapt_triangle*> all_elems; static int nbNod; }; @@ -133,38 +134,39 @@ public: class adapt_quad { public: - adapt_quad (adapt_point *p1,adapt_point *p2,adapt_point *p3,adapt_point *p4) + 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; + 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.; + 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); + 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); + 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 ); + 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]; + adapt_point *p[4]; + adapt_quad *e[4]; static std::list<adapt_quad*> all_elems; static int nbNod; }; @@ -172,18 +174,17 @@ public: class adapt_tet { public: - adapt_tet (adapt_point *p1,adapt_point *p2,adapt_point *p3,adapt_point *p4) + 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; + 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[]) + inline static void GSF (const double u, const double v, double w, double sf[]) { sf[0] = 1-u-v-w; sf[1] = u; @@ -192,19 +193,21 @@ public: } inline double V () const { - return (p[0]->val + p[1]->val + p[2]->val+ p[3]->val)/4.; + 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); + 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 ); + 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]; + adapt_point *p[4]; + adapt_tet *e[8]; static std::list<adapt_tet*> all_elems; static int nbNod; }; @@ -212,7 +215,8 @@ public: 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) + 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; @@ -223,35 +227,37 @@ public: 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; + 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[]) + 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); + 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.; + 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); + 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 ); + 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_point *p[8]; adapt_hex *e[8]; static std::list<adapt_hex*> all_elems; static int nbNod; @@ -273,21 +279,24 @@ class Adaptive_Post_View Double_Matrix * _Geometry; public: Adaptive_Post_View (Post_View *view, List_T *_coeffs, List_T *_eexps); - ~Adaptive_Post_View(); - int getGlobalResolutionLevel ( ) const {return presentZoomLevel;} - void setGlobalResolutionLevel ( Post_View * view , int level ) - { - setAdaptiveResolutionLevel ( view , level ); - } + ~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 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;} + 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); + int zoomElement (Post_View * view, int ielem, int level, int levelmax, + GMSH_Post_Plugin *plug, List_T *theList, int *counter); }; template <class ELEM> @@ -295,14 +304,11 @@ 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; - } + for (; it != ite; ++it){ + delete *it; + } ELEM::all_elems.clear(); adapt_point::all_points.clear(); - } - #endif -- GitLab