diff --git a/Common/AdaptiveViews.cpp b/Common/AdaptiveViews.cpp
index 5bd11e8ed737f644c434e6291ab1333252c72012..d237d4db7711a5328901f56d775da220e29bdbe4 100644
--- a/Common/AdaptiveViews.cpp
+++ b/Common/AdaptiveViews.cpp
@@ -35,6 +35,7 @@ void computeShapeFunctions ( Double_Matrix *coeffs, Double_Matrix *eexps , doubl
std::set<_point> _point::all_points;
std::list<_triangle*> _triangle::all_triangles;
std::list<_tet*> _tet::all_tets;
+std::list<_quad*> _quad::all_quads;
_point * _point::New ( double x, double y, double z, Double_Matrix *coeffs, Double_Matrix *eexps)
{
@@ -64,6 +65,18 @@ void _triangle::clean ()
_point::all_points.clear();
}
+void _quad::clean ()
+{
+ std::list<_quad*>::iterator it = all_quads.begin();
+ std::list<_quad*>::iterator ite = all_quads.end();
+ for (;it!=ite;++it)
+ {
+ delete *it;
+ }
+ all_quads.clear();
+ _point::all_points.clear();
+}
+
void _tet::clean ()
{
std::list<_tet*>::iterator it = all_tets.begin();
@@ -91,6 +104,21 @@ void _triangle::Create (int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexp
printf("%d _triangle and %d _point created\n",(int)_triangle::all_triangles.size(),(int)_point::all_points.size());
}
+void _quad::Create (int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps)
+{
+ printf("creating the sub-elements\n");
+ int level = 0;
+ clean();
+ _point *p1 = _point::New ( -1,-1,0, coeffs, eexps);
+ _point *p2 = _point::New ( -1,1,0, coeffs, eexps);
+ _point *p3 = _point::New ( 1,1,0, coeffs, eexps);
+ _point *p4 = _point::New ( 1,-1,0, coeffs, eexps);
+ _quad *q = new _quad(p1,p2,p3,p4);
+ Recur_Create (q, maxlevel,level,coeffs,eexps) ;
+
+ printf("%d _triangle and %d _point created\n",(int)_triangle::all_triangles.size(),(int)_point::all_points.size());
+}
+
void _tet::Create (int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps)
{
Msg(INFO, "creating sub-elements");
@@ -142,6 +170,46 @@ void _triangle::Recur_Create (_triangle *t, int maxlevel, int level , Double_Mat
}
+void _quad::Recur_Create (_quad *q, int maxlevel, int level , Double_Matrix *coeffs, Double_Matrix *eexps)
+{
+ all_quads.push_back(q);
+ if (level++ >= maxlevel)
+ return;
+
+ /*
+
+ p2 p23 p3
+
+
+ p12 pc p34
+
+
+ p1 p14 p4
+
+
+ */
+
+ _point *p1 = q->p[0];
+ _point *p2 = q->p[1];
+ _point *p3 = q->p[2];
+ _point *p4 = q->p[3];
+ _point *p12 = _point::New ( (p1->x+p2->x)*0.5,(p1->y+p2->y)*0.5,0, coeffs, eexps);
+ _point *p23 = _point::New ( (p2->x+p3->x)*0.5,(p2->y+p3->y)*0.5,0, coeffs, eexps);
+ _point *p34 = _point::New ( (p4->x+p3->x)*0.5,(p4->y+p3->y)*0.5,0, coeffs, eexps);
+ _point *p14 = _point::New ( (p1->x+p4->x)*0.5,(p1->y+p4->y)*0.5,0, coeffs, eexps);
+ _point *pc = _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);
+ _quad *q1 = new _quad (p1,p12,pc,p14);
+ Recur_Create (q1, maxlevel,level,coeffs,eexps);
+ _quad *q2 = new _quad (p12,p2,p23,pc);
+ Recur_Create (q2, maxlevel,level,coeffs,eexps);
+ _quad *q3 = new _quad (pc,p23,p3,p34);
+ Recur_Create (q3, maxlevel,level,coeffs,eexps);
+ _quad *q4 = new _quad (p14,pc,p34,p4);
+ Recur_Create (q4, maxlevel,level,coeffs,eexps);
+ q->q[0]=q1;q->q[1]=q2;q->q[2]=q3;q->q[3]=q4;
+
+}
+
void _tet::Recur_Create (_tet *t, int maxlevel, int level , Double_Matrix *coeffs, Double_Matrix *eexps)
{
all_tets.push_back(t);
@@ -188,6 +256,12 @@ void _triangle::Error ( double AVG , double tol )
Recur_Error (t,AVG,tol);
}
+void _quad::Error ( double AVG , double tol )
+{
+ _quad *q = *all_quads.begin();
+ Recur_Error (q,AVG,tol);
+}
+
void _tet::Error ( double AVG , double tol )
{
_tet *t = *all_tets.begin();
@@ -266,6 +340,78 @@ void _triangle::Recur_Error ( _triangle *t, double AVG, double tol )
}
}
+void _quad::Recur_Error ( _quad *q, double AVG, double tol )
+{
+ if(!q->q[0])q->visible = true;
+ else
+ {
+ double vr;
+ if (!q->q[0]->q[0])
+ {
+ double v1 = q->q[0]->V();
+ double v2 = q->q[1]->V();
+ double v3 = q->q[2]->V();
+ double v4 = q->q[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->q[0],AVG,tol);
+ Recur_Error (q->q[1],AVG,tol);
+ Recur_Error (q->q[2],AVG,tol);
+ Recur_Error (q->q[3],AVG,tol);
+ }
+ else
+ q->visible = true;
+ }
+ else
+ {
+ double v11 = q->q[0]->q[0]->V();
+ double v12 = q->q[0]->q[1]->V();
+ double v13 = q->q[0]->q[2]->V();
+ double v14 = q->q[0]->q[3]->V();
+ double v21 = q->q[1]->q[0]->V();
+ double v22 = q->q[1]->q[1]->V();
+ double v23 = q->q[1]->q[2]->V();
+ double v24 = q->q[1]->q[3]->V();
+ double v31 = q->q[2]->q[0]->V();
+ double v32 = q->q[2]->q[1]->V();
+ double v33 = q->q[2]->q[2]->V();
+ double v34 = q->q[2]->q[3]->V();
+ double v41 = q->q[3]->q[0]->V();
+ double v42 = q->q[3]->q[1]->V();
+ double v43 = q->q[3]->q[2]->V();
+ double v44 = q->q[3]->q[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->q[0]->V() - vr1) > AVG * tol ||
+ fabs(q->q[1]->V() - vr2) > AVG * tol ||
+ fabs(q->q[2]->V() - vr3) > AVG * tol ||
+ fabs(q->q[3]->V() - vr4) > AVG * tol ||
+ fabs(q->V() - vr) > AVG * tol )
+ //if ( fabs(t->t[0]->V() - vr1) > (fabs(t->t[0]->V())+fabs(vr1)+AVG * tol)*tol ||
+ // fabs(t->t[1]->V() - vr2) > (fabs(t->t[1]->V())+fabs(vr2)+AVG * tol)*tol ||
+ // fabs(t->t[2]->V() - vr3) > (fabs(t->t[2]->V())+fabs(vr3)+AVG * tol)*tol ||
+ // fabs(t->t[3]->V() - vr4) > (fabs(t->t[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->q[0],AVG,tol);
+ Recur_Error (q->q[1],AVG,tol);
+ Recur_Error (q->q[2],AVG,tol);
+ Recur_Error (q->q[3],AVG,tol);
+ }
+ else
+ q->visible = true;
+ }
+ }
+}
+
void _tet::Recur_Error ( _tet *t, double AVG, double tol )
{
@@ -302,10 +448,10 @@ void _tet::Recur_Error ( _tet *t, double AVG, double tol )
static double t0,t1,t2,t3;
-void Adaptive_Post_View:: zoomElement (Post_View * view ,
- int ielem ,
- int level,
- GMSH_Post_Plugin *plug)
+void Adaptive_Post_View:: zoomTriangle (Post_View * view ,
+ int ielem ,
+ int level,
+ GMSH_Post_Plugin *plug)
{
std::set<_point>::iterator it = _point::all_points.begin();
std::set<_point>::iterator ite = _point::all_points.end();
@@ -396,21 +542,122 @@ void Adaptive_Post_View:: zoomElement (Post_View * view ,
}
+void Adaptive_Post_View:: zoomQuad (Post_View * view ,
+ int ielem ,
+ int level,
+ GMSH_Post_Plugin *plug)
+{
+ std::set<_point>::iterator it = _point::all_points.begin();
+ std::set<_point>::iterator ite = _point::all_points.end();
+
+ double c0 = Cpu();
+
+ const int N = _coefs->size1();
+ Double_Vector val ( N ), res(_point::all_points.size());
+ Double_Matrix xyz (4,3);
+ Double_Matrix XYZ (_point::all_points.size(),3);
+
+ for ( int k=0;k<4;++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);
+
+ double c1 = Cpu();
+
+ int kk=0;
+ for ( ; it !=ite ; ++it)
+ {
+ _point *p = (_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();
+
+ std::list<_quad*>::iterator itt = _quad::all_quads.begin();
+ std::list<_quad*>::iterator itte = _quad::all_quads.end();
+
+ for ( ;itt != itte ; itt++)
+ {
+ (*itt)->visible = false;
+ }
+
+
+ if (plug)
+ plug->assign_specific_visibility ();
+ else
+ {
+ _quad::Error ( max-min, tol );
+ }
+ double c3 = Cpu();
+ itt = _quad::all_quads.begin();
+ for ( ;itt != itte ; itt++)
+ {
+ if ((*itt)->visible)
+ {
+ _point *p1 = (*itt)->p[0];
+ _point *p2 = (*itt)->p[1];
+ _point *p3 = (*itt)->p[2];
+ _point *p4 = (*itt)->p[3];
+ List_Add ( view->SQ , &p1->X );
+ List_Add ( view->SQ , &p2->X );
+ List_Add ( view->SQ , &p3->X );
+ List_Add ( view->SQ , &p4->X );
+ List_Add ( view->SQ , &p1->Y );
+ List_Add ( view->SQ , &p2->Y );
+ List_Add ( view->SQ , &p3->Y );
+ List_Add ( view->SQ , &p4->Y );
+ List_Add ( view->SQ , &p1->Z );
+ List_Add ( view->SQ , &p2->Z );
+ List_Add ( view->SQ , &p3->Z );
+ List_Add ( view->SQ , &p4->Z );
+ List_Add ( view->SQ , &p1->val );
+ List_Add ( view->SQ , &p2->val );
+ List_Add ( view->SQ , &p3->val );
+ List_Add ( view->SQ , &p4->val );
+ view->NbSQ++;
+ }
+ }
+ double c4 = Cpu();
+
+ t0 += c1-c0;
+ t1 += c2-c1;
+ t2 += c3-c2;
+ t3 += c4-c3;
+}
+
void Adaptive_Post_View:: zoomTet (Post_View * view ,
int ielem ,
int level,
- GMSH_Post_Plugin *plug)
+ GMSH_Post_Plugin *plug,
+ Double_Vector & va2l,
+ Double_Vector & re2s,
+ Double_Matrix & XY2Z)
{
std::set<_point>::iterator it = _point::all_points.begin();
std::set<_point>::iterator ite = _point::all_points.end();
+
+
+
double c0 = Cpu();
- const int N = _coefs->size1();
- Double_Vector val ( N ), res(_point::all_points.size());
Double_Matrix xyz (4,3);
- Double_Matrix XYZ (_point::all_points.size(),3);
for ( int k=0;k<4;++k)
{
@@ -418,6 +665,10 @@ void Adaptive_Post_View:: zoomTet (Post_View * view ,
xyz(k,1) = (*_STposY) ( ielem , k );
xyz(k,2) = (*_STposZ) ( ielem , k );
}
+
+ const int N = _coefs->size1();
+ Double_Vector val ( N ), res(_point::all_points.size());
+ Double_Matrix XYZ (_point::all_points.size(),3);
for ( int k=0;k<N;++k)
{
@@ -500,105 +751,153 @@ void Adaptive_Post_View:: zoomTet (Post_View * view ,
void Adaptive_Post_View:: setAdaptiveResolutionLevel (Post_View * view , int level, GMSH_Post_Plugin *plug)
{
- if (presentTol==tol && presentZoomLevel == level && !plug)return;
- if (view->NbST)
+ const int N = _coefs->size1();
+
+ if (presentTol==tol && presentZoomLevel == level && !plug)return;
+
+ if (view->NbST)
{
- _triangle::Create ( level, _coefs, _eexps );
- std::set<_point>::iterator it = _point::all_points.begin();
- std::set<_point>::iterator ite = _point::all_points.end();
-
- const int N = _coefs->size1();
- if (_Interpolate)
- delete _Interpolate;
- if (_Geometry)
- delete _Geometry;
- _Interpolate = new Double_Matrix ( _point::all_points.size(), N);
- _Geometry = new Double_Matrix ( _point::all_points.size(), 3);
-
- int kk=0;
- for ( ; it !=ite ; ++it)
+ _triangle::Create ( level, _coefs, _eexps );
+ std::set<_point>::iterator it = _point::all_points.begin();
+ std::set<_point>::iterator ite = _point::all_points.end();
+
+ if (_Interpolate)
+ delete _Interpolate;
+ if (_Geometry)
+ delete _Geometry;
+ _Interpolate = new Double_Matrix ( _point::all_points.size(), N);
+ _Geometry = new Double_Matrix ( _point::all_points.size(), 3);
+
+ int kk=0;
+ for ( ; it !=ite ; ++it)
{
- _point *p = (_point*) &(*it);
- for ( int k=0;k<N;++k)
+ _point *p = (_point*) &(*it);
+ for ( int k=0;k<N;++k)
{
- (*_Interpolate)(kk,k) = p->shape_functions[k];
+ (*_Interpolate)(kk,k) = p->shape_functions[k];
}
- (*_Geometry)(kk,0) = ( 1.-p->x-p->y);
- (*_Geometry)(kk,1) = p->x;
- (*_Geometry)(kk,2) = p->y;
- kk++;
+ (*_Geometry)(kk,0) = ( 1.-p->x-p->y);
+ (*_Geometry)(kk,1) = p->x;
+ (*_Geometry)(kk,2) = p->y;
+ kk++;
}
-
- List_Delete(view->ST); view->ST = 0;
- view->NbST = 0;
- /// for now, that's all we do, 1 TS
- view->NbTimeStep=1;
- int nbelm = _STposX->size1();
- view->ST = List_Create ( nbelm * 4, nbelm , sizeof(double));
-
-
- t0=t1 = t2 = t3 = 0;
-
- for ( int i=0;i<nbelm;++i)
+
+ List_Delete(view->ST); view->ST = 0;
+ view->NbST = 0;
+ /// for now, that's all we do, 1 TS
+ view->NbTimeStep=1;
+ int nbelm = _STposX->size1();
+ view->ST = List_Create ( nbelm * 4, nbelm , sizeof(double));
+
+
+ t0=t1 = t2 = t3 = 0;
+
+ for ( int i=0;i<nbelm;++i)
{
- zoomElement ( view , i , level, plug);
+ zoomTriangle ( view , i , level, plug);
}
-
- printf("finished %g %g %g %g\n",t0,t1,t2,t3);
+
+ printf("finished %g %g %g %g\n",t0,t1,t2,t3);
}
- else if (view->NbSS)
+ else if (view->NbSS)
{
- _tet::Create ( level, _coefs, _eexps );
- std::set<_point>::iterator it = _point::all_points.begin();
- std::set<_point>::iterator ite = _point::all_points.end();
-
- const int N = _coefs->size1();
- if (_Interpolate)
- delete _Interpolate;
- if (_Geometry)
- delete _Geometry;
- _Interpolate = new Double_Matrix ( _point::all_points.size(), N);
- _Geometry = new Double_Matrix ( _point::all_points.size(), 4);
-
- int kk=0;
- for ( ; it !=ite ; ++it)
+ _tet::Create ( level, _coefs, _eexps );
+ std::set<_point>::iterator it = _point::all_points.begin();
+ std::set<_point>::iterator ite = _point::all_points.end();
+
+ if (_Interpolate)
+ delete _Interpolate;
+ if (_Geometry)
+ delete _Geometry;
+ _Interpolate = new Double_Matrix ( _point::all_points.size(), N);
+ _Geometry = new Double_Matrix ( _point::all_points.size(), 4);
+
+ int kk=0;
+ for ( ; it !=ite ; ++it)
{
- _point *p = (_point*) &(*it);
- for ( int k=0;k<N;++k)
+ _point *p = (_point*) &(*it);
+ for ( int k=0;k<N;++k)
{
- (*_Interpolate)(kk,k) = p->shape_functions[k];
+ (*_Interpolate)(kk,k) = p->shape_functions[k];
}
- (*_Geometry)(kk,0) = ( 1.-p->x-p->y-p->z);
- (*_Geometry)(kk,1) = p->x;
- (*_Geometry)(kk,2) = p->y;
- (*_Geometry)(kk,3) = p->z;
- kk++;
+ (*_Geometry)(kk,0) = ( 1.-p->x-p->y-p->z);
+ (*_Geometry)(kk,1) = p->x;
+ (*_Geometry)(kk,2) = p->y;
+ (*_Geometry)(kk,3) = p->z;
+ kk++;
}
-
- List_Delete(view->SS); view->SS = 0;
- view->NbSS = 0;
- /// for now, that's all we do, 1 TS
- view->NbTimeStep=1;
- int nbelm = _STposX->size1();
- view->SS = List_Create ( nbelm * 4, nbelm , sizeof(double));
-
-
- t0=t1 = t2 = t3 = 0;
-
- for ( int i=0;i<nbelm;++i)
+
+ List_Delete(view->SS); view->SS = 0;
+ view->NbSS = 0;
+ /// for now, that's all we do, 1 TS
+ view->NbTimeStep=1;
+ int nbelm = _STposX->size1();
+ view->SS = List_Create ( nbelm * 4, nbelm , sizeof(double));
+
+
+ t0=t1 = t2 = t3 = 0;
+
+ Double_Vector val ( N ), res(_point::all_points.size());
+ Double_Matrix XYZ (_point::all_points.size(),3);
+
+ for ( int i=0;i<nbelm;++i)
{
- zoomTet ( view , i , level, plug);
+ zoomTet ( view , i , level, plug,val,res,XYZ);
}
-
- printf("finished B %g %g %g %g\n",t0,t1,t2,t3);
+
+ printf("finished B %g %g %g %g\n",t0,t1,t2,t3);
}
- else
- return;
-
- view->Changed = 1;
- presentZoomLevel = level;
- presentTol = tol;
-
+ else if (view->NbSQ)
+ {
+ _quad::Create ( level, _coefs, _eexps );
+ std::set<_point>::iterator it = _point::all_points.begin();
+ std::set<_point>::iterator ite = _point::all_points.end();
+
+ if (_Interpolate)
+ delete _Interpolate;
+ if (_Geometry)
+ delete _Geometry;
+ _Interpolate = new Double_Matrix ( _point::all_points.size(), N);
+ _Geometry = new Double_Matrix ( _point::all_points.size(), 4);
+
+ int kk=0;
+ for ( ; it !=ite ; ++it)
+ {
+ _point *p = (_point*) &(*it);
+ for ( int k=0;k<N;++k)
+ {
+ (*_Interpolate)(kk,k) = p->shape_functions[k];
+ }
+ (*_Geometry)(kk,0) = ( 1.-p->x) * ( 1.-p->y) * .25;
+ (*_Geometry)(kk,1) = ( 1.+p->x) * ( 1.-p->y) * .25;
+ (*_Geometry)(kk,2) = ( 1.+p->x) * ( 1.+p->y) * .25;
+ (*_Geometry)(kk,3) = ( 1.-p->x) * ( 1.+p->y) * .25;
+ kk++;
+ }
+
+ List_Delete(view->SQ); view->SQ = 0;
+ view->NbSQ = 0;
+ /// for now, that's all we do, 1 TS
+ view->NbTimeStep=1;
+ int nbelm = _STposX->size1();
+ view->SQ = List_Create ( nbelm * 4, nbelm , sizeof(double));
+
+
+ t0=t1 = t2 = t3 = 0;
+
+ for ( int i=0;i<nbelm;++i)
+ {
+ zoomQuad ( view , i , level, plug);
+ }
+
+ printf("finished Q %g %g %g %g\n",t0,t1,t2,t3);
+ }
+ else
+ return;
+
+ view->Changed = 1;
+ presentZoomLevel = level;
+ presentTol = tol;
}
void computeShapeFunctions ( Double_Matrix *coeffs, Double_Matrix *eexps , double u, double v, double w,double *sf)
@@ -635,12 +934,19 @@ void Adaptive_Post_View:: initWithLowResolution (Post_View *view)
nbelm = view->NbST;
nbnod = 3;
}
- else
+ 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 return;
min = VAL_INF;
max = -VAL_INF;
diff --git a/Common/Views.h b/Common/Views.h
index 5fc2e2893700b5098aab6a4e3c0a16745979d1f8..3edaa691dd7e1d31e77d3bd49b05fc87f866c07e 100644
--- a/Common/Views.h
+++ b/Common/Views.h
@@ -95,6 +95,38 @@ public:
static std::list<_triangle*> all_triangles;
};
+class _quad
+{
+public:
+ _quad (_point *p1,_point *p2,_point *p3,_point *p4)
+ : visible (false)
+ {
+ p[0] = p1;
+ p[1] = p2;
+ p[2] = p3;
+ p[3] = p4;
+ q[0]=q[1]=q[2]=q[3]=0;
+ }
+
+ 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 clean ();
+ static void Create (int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps) ;
+ static void Recur_Create (_quad *q, int maxlevel, int level , Double_Matrix *coeffs, Double_Matrix *eexps);
+ static void Error ( double AVG , double tol );
+ static void Recur_Error ( _quad *q, double AVG, double tol );
+ bool visible;
+ _point *p[4];
+ _quad *q[4];
+ static std::list<_quad*> all_quads;
+};
+
class _tet
{
public:
@@ -155,11 +187,15 @@ public:
void initWithLowResolution (Post_View *view);
void setTolerance (const double eps) {tol=eps;}
double getTolerance () const {return tol;}
- void zoomElement (Post_View * view ,
+ void zoomQuad (Post_View * view ,
+ int ielem, int level, GMSH_Post_Plugin *plug);
+ void zoomTriangle (Post_View * view ,
int ielem, int level, GMSH_Post_Plugin *plug);
void zoomTet (Post_View * view ,
- int ielem, int level, GMSH_Post_Plugin *plug);
-
+ int ielem, int level, GMSH_Post_Plugin *plug,
+ Double_Vector & val,
+ Double_Vector & res,
+ Double_Matrix & XYZ);
};
class Post_View{
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 472b6ee4ed56fb25423759b4e099e608c89e910b..722fc463a6a9153fcefda35358d842092fa9c839 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -1,4 +1,4 @@
-// $Id: Numeric.cpp,v 1.17 2004-07-21 22:19:56 geuzaine Exp $
+// $Id: Numeric.cpp,v 1.18 2004-11-25 16:22:45 remacle Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -281,3 +281,123 @@ void gradSimplex(double *x, double *y, double *z, double *v, double *grad)
b[2] = v[3] - v[0];
sys3x3(mat, b, grad, &det);
}
+
+void eigenvalue(double mat[3][3], double v[3])
+ {
+ /// characteristic polynomial of T : find v root of
+ /// v^3 - I1 v^2 + I2 T + I3 = 0
+ /// I1 : first invariant , trace(T)
+ /// I2 : second invariant , 1/2 (I1^2 -trace(T^2))
+ /// I3 : third invariant , det T
+
+ double In[4];
+ In[3] = 1.0;
+ In[2] = - trace(mat);
+ In[1] = 0.5 * (In[2]*In[2] - trace2(mat));
+ In[0] = - det(mat);
+
+ // printf (" %lf x^3 + %lf x^2 + %lf x + %lf = 0\n",
+ // I[3],I[2],I[1],I[0]);
+
+ FindCubicRoots(In,v);
+
+ eigsort(v);
+
+ }
+
+void FindCubicRoots(const double coeff[4], double x[3])
+ {
+ double a1 = coeff[2] / coeff[3];
+ double a2 = coeff[1] / coeff[3];
+ double a3 = coeff[0] / coeff[3];
+
+ double Q = (a1 * a1 - 3 * a2) / 9.;
+ double R = (2. * a1 * a1 * a1 - 9. * a1 * a2 + 27. * a3) / 54.;
+ double Qcubed = Q * Q * Q;
+ double d = Qcubed - R * R;
+
+ // printf ("d = %22.15e Q = %12.5E R = %12.5E Qcubed %12.5E\n",d,Q,R,Qcubed);
+
+ /// three roots, 2 equal
+ if(Qcubed == 0.0 || fabs ( Qcubed - R * R ) < 1.e-8 * (fabs ( Qcubed) + fabs( R * R)) )
+ {
+ double theta;
+ if (Qcubed <= 0.0)theta = acos(1.0);
+ else if (R / sqrt(Qcubed) > 1.0)theta = acos(1.0);
+ else if (R / sqrt(Qcubed) < -1.0)theta = acos(-1.0);
+ else theta = acos(R / sqrt(Qcubed));
+ double sqrtQ = sqrt(Q);
+ // printf("sqrtQ = %12.5E teta=%12.5E a1=%12.5E\n",sqrt(Q),theta,a1);
+ x[0] = -2 * sqrtQ * cos( theta / 3) - a1 / 3;
+ x[1] = -2 * sqrtQ * cos((theta + 2 * M_PI) / 3) - a1 / 3;
+ x[2] = -2 * sqrtQ * cos((theta + 4 * M_PI) / 3) - a1 / 3;
+ // return (3);
+ }
+
+ /* Three real roots */
+ if (d >= 0.0) {
+ double theta = acos(R / sqrt(Qcubed));
+ double sqrtQ = sqrt(Q);
+ x[0] = -2 * sqrtQ * cos( theta / 3) - a1 / 3;
+ x[1] = -2 * sqrtQ * cos((theta + 2 * M_PI) / 3) - a1 / 3;
+ x[2] = -2 * sqrtQ * cos((theta + 4 * M_PI) / 3) - a1 / 3;
+ // return (3);
+ }
+
+ /* One real root */
+ else {
+ printf("IMPOSSIBLE !!!\n");
+
+ double e = pow(sqrt(-d) + fabs(R), 1. / 3.);
+ if (R > 0)
+ e = -e;
+ x[0] = (e + Q / e) - a1 / 3.;
+ // return (1);
+ }
+
+ }
+
+double trace(double mat[3][3])
+ {
+ return mat[0][0] + mat[1][1] + mat[2][2];
+ }
+
+
+double det(double mat[3][3])
+ {
+ return mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
+ mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
+ mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
+ }
+
+double trace2 (double mat[3][3])
+ {
+ double a00 = mat[0][0] * mat[0][0] +
+ mat[1][0] * mat[0][1] +
+ mat[2][0] * mat[0][2];
+ double a11 = mat[1][0] * mat[0][1] +
+ mat[1][1] * mat[1][1] +
+ mat[1][2] * mat[2][1];
+ double a22 = mat[2][0] * mat[0][2] +
+ mat[2][1] * mat[1][2] +
+ mat[2][2] * mat[2][2];
+
+ return a00 + a11 + a22;
+ }
+
+void eigsort(double d[3])
+
+{
+ int k,j,i;
+ double p;
+
+ for (i=0; i<3; i++) {
+ p=d[k=i];
+ for (j=i+1;j<3;j++)
+ if (d[j] >= p) p=d[k=j];
+ if (k != i) {
+ d[k]=d[i];
+ d[i]=p;
+ }
+ }
+}
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 6ada0a0bf1cb46b4b84ca0a5c32b7abf4806b112..b3fcd9664b82075ed1008e82ee75f4c7f5416d36 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -84,5 +84,11 @@ void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb,
double *fc, double (*func)(double));
void newt(double x[], int n, int *check,
void (*vecfunc)(int, double [], double []));
+void eigenvalue(double mat[3][3], double v[3]);
+void FindCubicRoots(const double coeff[4], double x[3]);
+double trace(double mat[3][3]);
+double det(double mat[3][3]);
+double trace2 (double mat[3][3]);
+void eigsort(double d[3]);
#endif
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index 933bcb3f8de0701aae517a1f16f0f63da72b6a5a..2234d6c644832bc815556036af13f4e085425437 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -1,7 +1,7 @@
-/* A Bison parser, made by GNU Bison 1.875c. */
+/* A Bison parser, made by GNU Bison 1.875. */
/* Skeleton parser for Yacc-like parsing with Bison,
- Copyright (C) 1984, 1989, 1990, 2000, 2001, 2002, 2003 Free Software Foundation, Inc.
+ Copyright (C) 1984, 1989, 1990, 2000, 2001, 2002 Free Software Foundation, Inc.
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
@@ -434,7 +434,7 @@
/* Copy the first part of user declarations. */
#line 1 "Gmsh.y"
-// $Id: Gmsh.tab.cpp,v 1.210 2004-11-25 02:10:32 geuzaine Exp $
+// $Id: Gmsh.tab.cpp,v 1.211 2004-11-25 16:22:45 remacle Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -538,7 +538,7 @@ typedef union YYSTYPE {
List_T *l;
} YYSTYPE;
/* Line 191 of yacc.c. */
-#line 542 "Gmsh.tab.cpp"
+#line 541 "Gmsh.tab.cpp"
# define yystype YYSTYPE /* obsolescent; will be withdrawn */
# define YYSTYPE_IS_DECLARED 1
# define YYSTYPE_IS_TRIVIAL 1
@@ -550,29 +550,22 @@ typedef union YYSTYPE {
/* Line 214 of yacc.c. */
-#line 554 "Gmsh.tab.cpp"
+#line 553 "Gmsh.tab.cpp"
#if ! defined (yyoverflow) || YYERROR_VERBOSE
-# ifndef YYFREE
-# define YYFREE free
-# endif
-# ifndef YYMALLOC
-# define YYMALLOC malloc
-# endif
-
/* The parser invokes alloca or malloc; define the necessary symbols. */
-# ifdef YYSTACK_USE_ALLOCA
-# if YYSTACK_USE_ALLOCA
-# define YYSTACK_ALLOC alloca
-# endif
+# if YYSTACK_USE_ALLOCA
+# define YYSTACK_ALLOC alloca
# else
-# if defined (alloca) || defined (_ALLOCA_H)
-# define YYSTACK_ALLOC alloca
-# else
-# ifdef __GNUC__
-# define YYSTACK_ALLOC __builtin_alloca
+# ifndef YYSTACK_USE_ALLOCA
+# if defined (alloca) || defined (_ALLOCA_H)
+# define YYSTACK_ALLOC alloca
+# else
+# ifdef __GNUC__
+# define YYSTACK_ALLOC __builtin_alloca
+# endif
# endif
# endif
# endif
@@ -585,15 +578,15 @@ typedef union YYSTYPE {
# include <stdlib.h> /* INFRINGES ON USER NAME SPACE */
# define YYSIZE_T size_t
# endif
-# define YYSTACK_ALLOC YYMALLOC
-# define YYSTACK_FREE YYFREE
+# define YYSTACK_ALLOC malloc
+# define YYSTACK_FREE free
# endif
#endif /* ! defined (yyoverflow) || YYERROR_VERBOSE */
#if (! defined (yyoverflow) \
&& (! defined (__cplusplus) \
- || (defined (YYSTYPE_IS_TRIVIAL) && YYSTYPE_IS_TRIVIAL)))
+ || (YYSTYPE_IS_TRIVIAL)))
/* A type that is properly aligned for any stack member. */
union yyalloc
@@ -614,7 +607,7 @@ union yyalloc
/* Copy COUNT objects from FROM to TO. The source and destination do
not overlap. */
# ifndef YYCOPY
-# if defined (__GNUC__) && 1 < __GNUC__
+# if 1 < __GNUC__
# define YYCOPY(To, From, Count) \
__builtin_memcpy (To, From, (Count) * sizeof (*(From)))
# else
@@ -1176,79 +1169,80 @@ static const unsigned short yyrline[] =
First, the terminals, then, starting at YYNTOKENS, nonterminals. */
static const char *const yytname[] =
{
- "$end", "error", "$undefined", "tDOUBLE", "tSTRING", "tBIGSTR", "tEND",
- "tAFFECT", "tDOTS", "tPi", "tMPI_Rank", "tMPI_Size", "tExp", "tLog",
- "tLog10", "tSqrt", "tSin", "tAsin", "tCos", "tAcos", "tTan", "tRand",
- "tAtan", "tAtan2", "tSinh", "tCosh", "tTanh", "tFabs", "tFloor", "tCeil",
- "tFmod", "tModulo", "tHypot", "tPrintf", "tSprintf", "tStrCat",
- "tStrPrefix", "tBoundingBox", "tDraw", "tPoint", "tCircle", "tEllipse",
- "tLine", "tSurface", "tSpline", "tVolume", "tCharacteristic", "tLength",
- "tParametric", "tElliptic", "tPlane", "tRuled", "tTriangulation",
- "tTransfinite", "tComplex", "tPhysical", "tUsing", "tBump",
- "tProgression", "tPlugin", "tRotate", "tTranslate", "tSymmetry",
- "tDilate", "tExtrude", "tDuplicata", "tLoop", "tRecombine", "tDelete",
- "tCoherence", "tIntersect", "tAttractor", "tLayers", "tScalarPoint",
- "tVectorPoint", "tTensorPoint", "tScalarLine", "tVectorLine",
- "tTensorLine", "tScalarTriangle", "tVectorTriangle", "tTensorTriangle",
- "tScalarQuadrangle", "tVectorQuadrangle", "tTensorQuadrangle",
- "tScalarTetrahedron", "tVectorTetrahedron", "tTensorTetrahedron",
- "tScalarHexahedron", "tVectorHexahedron", "tTensorHexahedron",
- "tScalarPrism", "tVectorPrism", "tTensorPrism", "tScalarPyramid",
- "tVectorPyramid", "tTensorPyramid", "tText2D", "tText3D",
- "tInterpolationScheme", "tCombine", "tBSpline", "tBezier", "tNurbs",
- "tOrder", "tWith", "tBounds", "tKnots", "tColor", "tColorTable", "tFor",
- "tIn", "tEndFor", "tIf", "tEndIf", "tExit", "tReturn", "tCall",
- "tFunction", "tTrimmed", "tShow", "tHide",
- "tB_SPLINE_SURFACE_WITH_KNOTS", "tB_SPLINE_CURVE_WITH_KNOTS",
- "tCARTESIAN_POINT", "tTRUE", "tFALSE", "tUNSPECIFIED", "tU", "tV",
- "tEDGE_CURVE", "tVERTEX_POINT", "tORIENTED_EDGE", "tPLANE",
- "tFACE_OUTER_BOUND", "tEDGE_LOOP", "tADVANCED_FACE", "tVECTOR",
- "tDIRECTION", "tAXIS2_PLACEMENT_3D", "tISO", "tENDISO", "tENDSEC",
- "tDATA", "tHEADER", "tFILE_DESCRIPTION", "tFILE_SCHEMA", "tFILE_NAME",
- "tMANIFOLD_SOLID_BREP", "tCLOSED_SHELL",
- "tADVANCED_BREP_SHAPE_REPRESENTATION", "tFACE_BOUND",
- "tCYLINDRICAL_SURFACE", "tCONICAL_SURFACE", "tCIRCLE", "tTRIMMED_CURVE",
- "tGEOMETRIC_SET", "tCOMPOSITE_CURVE_SEGMENT", "tCONTINUOUS",
- "tCOMPOSITE_CURVE", "tTOROIDAL_SURFACE", "tPRODUCT_DEFINITION",
- "tPRODUCT_DEFINITION_SHAPE", "tSHAPE_DEFINITION_REPRESENTATION",
- "tELLIPSE", "tSolid", "tEndSolid", "tVertex", "tFacet", "tNormal",
- "tOuter", "tLoopSTL", "tEndLoop", "tEndFacet", "tAFFECTDIVIDE",
- "tAFFECTTIMES", "tAFFECTMINUS", "tAFFECTPLUS", "'?'", "tOR", "tAND",
- "tAPPROXEQUAL", "tNOTEQUAL", "tEQUAL", "'<'", "'>'", "tGREATEROREQUAL",
- "tLESSOREQUAL", "'+'", "'-'", "'*'", "'/'", "'%'", "tCROSSPRODUCT",
- "'!'", "UNARYPREC", "tMINUSMINUS", "tPLUSPLUS", "'^'", "'('", "')'",
- "'['", "']'", "'.'", "'#'", "','", "'{'", "'}'", "$accept", "All",
- "SignedDouble", "StlFormatItems", "StlFormatItem", "StepFormatItems",
- "StepFormatItem", "StepSpecial", "StepHeaderItem", "StepDataItem",
- "GeoFormatItems", "GeoFormatItem", "Printf", "View", "Views",
- "ScalarPointValues", "ScalarPoint", "@1", "VectorPointValues",
- "VectorPoint", "@2", "TensorPointValues", "TensorPoint", "@3",
- "ScalarLineValues", "ScalarLine", "@4", "VectorLineValues", "VectorLine",
- "@5", "TensorLineValues", "TensorLine", "@6", "ScalarTriangleValues",
- "ScalarTriangle", "@7", "VectorTriangleValues", "VectorTriangle", "@8",
- "TensorTriangleValues", "TensorTriangle", "@9", "ScalarQuadrangleValues",
- "ScalarQuadrangle", "@10", "VectorQuadrangleValues", "VectorQuadrangle",
- "@11", "TensorQuadrangleValues", "TensorQuadrangle", "@12",
- "ScalarTetrahedronValues", "ScalarTetrahedron", "@13",
- "VectorTetrahedronValues", "VectorTetrahedron", "@14",
- "TensorTetrahedronValues", "TensorTetrahedron", "@15",
- "ScalarHexahedronValues", "ScalarHexahedron", "@16",
- "VectorHexahedronValues", "VectorHexahedron", "@17",
- "TensorHexahedronValues", "TensorHexahedron", "@18", "ScalarPrismValues",
- "ScalarPrism", "@19", "VectorPrismValues", "VectorPrism", "@20",
- "TensorPrismValues", "TensorPrism", "@21", "ScalarPyramidValues",
- "ScalarPyramid", "@22", "VectorPyramidValues", "VectorPyramid", "@23",
- "TensorPyramidValues", "TensorPyramid", "@24", "Text2DValues", "Text2D",
- "@25", "Text3DValues", "Text3D", "@26", "InterpolationMatrix",
- "NumericAffectation", "NumericIncrement", "Affectation", "Shape",
- "Transform", "MultipleShape", "ListOfShapes", "Duplicata", "Delete",
- "Colorify", "Visibility", "Command", "Loop", "Extrude", "@27", "@28",
- "@29", "@30", "@31", "@32", "@33", "@34", "@35", "ExtrudeParameters",
- "ExtrudeParameter", "Transfinite", "Coherence", "BoolExpr", "FExpr",
- "FExpr_Single", "VExpr", "VExpr_Single", "ListOfStrings",
- "RecursiveListOfStrings", "ListOfListOfDouble",
- "RecursiveListOfListOfDouble", "ListOfDouble", "FExpr_Multi",
- "RecursiveListOfDouble", "ColorExpr", "ListOfColor",
+ "$end", "error", "$undefined", "tDOUBLE", "tSTRING", "tBIGSTR", "tEND",
+ "tAFFECT", "tDOTS", "tPi", "tMPI_Rank", "tMPI_Size", "tExp", "tLog",
+ "tLog10", "tSqrt", "tSin", "tAsin", "tCos", "tAcos", "tTan", "tRand",
+ "tAtan", "tAtan2", "tSinh", "tCosh", "tTanh", "tFabs", "tFloor",
+ "tCeil", "tFmod", "tModulo", "tHypot", "tPrintf", "tSprintf", "tStrCat",
+ "tStrPrefix", "tBoundingBox", "tDraw", "tPoint", "tCircle", "tEllipse",
+ "tLine", "tSurface", "tSpline", "tVolume", "tCharacteristic", "tLength",
+ "tParametric", "tElliptic", "tPlane", "tRuled", "tTriangulation",
+ "tTransfinite", "tComplex", "tPhysical", "tUsing", "tBump",
+ "tProgression", "tPlugin", "tRotate", "tTranslate", "tSymmetry",
+ "tDilate", "tExtrude", "tDuplicata", "tLoop", "tRecombine", "tDelete",
+ "tCoherence", "tIntersect", "tAttractor", "tLayers", "tScalarPoint",
+ "tVectorPoint", "tTensorPoint", "tScalarLine", "tVectorLine",
+ "tTensorLine", "tScalarTriangle", "tVectorTriangle", "tTensorTriangle",
+ "tScalarQuadrangle", "tVectorQuadrangle", "tTensorQuadrangle",
+ "tScalarTetrahedron", "tVectorTetrahedron", "tTensorTetrahedron",
+ "tScalarHexahedron", "tVectorHexahedron", "tTensorHexahedron",
+ "tScalarPrism", "tVectorPrism", "tTensorPrism", "tScalarPyramid",
+ "tVectorPyramid", "tTensorPyramid", "tText2D", "tText3D",
+ "tInterpolationScheme", "tCombine", "tBSpline", "tBezier", "tNurbs",
+ "tOrder", "tWith", "tBounds", "tKnots", "tColor", "tColorTable", "tFor",
+ "tIn", "tEndFor", "tIf", "tEndIf", "tExit", "tReturn", "tCall",
+ "tFunction", "tTrimmed", "tShow", "tHide",
+ "tB_SPLINE_SURFACE_WITH_KNOTS", "tB_SPLINE_CURVE_WITH_KNOTS",
+ "tCARTESIAN_POINT", "tTRUE", "tFALSE", "tUNSPECIFIED", "tU", "tV",
+ "tEDGE_CURVE", "tVERTEX_POINT", "tORIENTED_EDGE", "tPLANE",
+ "tFACE_OUTER_BOUND", "tEDGE_LOOP", "tADVANCED_FACE", "tVECTOR",
+ "tDIRECTION", "tAXIS2_PLACEMENT_3D", "tISO", "tENDISO", "tENDSEC",
+ "tDATA", "tHEADER", "tFILE_DESCRIPTION", "tFILE_SCHEMA", "tFILE_NAME",
+ "tMANIFOLD_SOLID_BREP", "tCLOSED_SHELL",
+ "tADVANCED_BREP_SHAPE_REPRESENTATION", "tFACE_BOUND",
+ "tCYLINDRICAL_SURFACE", "tCONICAL_SURFACE", "tCIRCLE", "tTRIMMED_CURVE",
+ "tGEOMETRIC_SET", "tCOMPOSITE_CURVE_SEGMENT", "tCONTINUOUS",
+ "tCOMPOSITE_CURVE", "tTOROIDAL_SURFACE", "tPRODUCT_DEFINITION",
+ "tPRODUCT_DEFINITION_SHAPE", "tSHAPE_DEFINITION_REPRESENTATION",
+ "tELLIPSE", "tSolid", "tEndSolid", "tVertex", "tFacet", "tNormal",
+ "tOuter", "tLoopSTL", "tEndLoop", "tEndFacet", "tAFFECTDIVIDE",
+ "tAFFECTTIMES", "tAFFECTMINUS", "tAFFECTPLUS", "'?'", "tOR", "tAND",
+ "tAPPROXEQUAL", "tNOTEQUAL", "tEQUAL", "'<'", "'>'", "tGREATEROREQUAL",
+ "tLESSOREQUAL", "'+'", "'-'", "'*'", "'/'", "'%'", "tCROSSPRODUCT",
+ "'!'", "UNARYPREC", "tMINUSMINUS", "tPLUSPLUS", "'^'", "'('", "')'",
+ "'['", "']'", "'.'", "'#'", "','", "'{'", "'}'", "$accept", "All",
+ "SignedDouble", "StlFormatItems", "StlFormatItem", "StepFormatItems",
+ "StepFormatItem", "StepSpecial", "StepHeaderItem", "StepDataItem",
+ "GeoFormatItems", "GeoFormatItem", "Printf", "View", "Views",
+ "ScalarPointValues", "ScalarPoint", "@1", "VectorPointValues",
+ "VectorPoint", "@2", "TensorPointValues", "TensorPoint", "@3",
+ "ScalarLineValues", "ScalarLine", "@4", "VectorLineValues",
+ "VectorLine", "@5", "TensorLineValues", "TensorLine", "@6",
+ "ScalarTriangleValues", "ScalarTriangle", "@7", "VectorTriangleValues",
+ "VectorTriangle", "@8", "TensorTriangleValues", "TensorTriangle", "@9",
+ "ScalarQuadrangleValues", "ScalarQuadrangle", "@10",
+ "VectorQuadrangleValues", "VectorQuadrangle", "@11",
+ "TensorQuadrangleValues", "TensorQuadrangle", "@12",
+ "ScalarTetrahedronValues", "ScalarTetrahedron", "@13",
+ "VectorTetrahedronValues", "VectorTetrahedron", "@14",
+ "TensorTetrahedronValues", "TensorTetrahedron", "@15",
+ "ScalarHexahedronValues", "ScalarHexahedron", "@16",
+ "VectorHexahedronValues", "VectorHexahedron", "@17",
+ "TensorHexahedronValues", "TensorHexahedron", "@18",
+ "ScalarPrismValues", "ScalarPrism", "@19", "VectorPrismValues",
+ "VectorPrism", "@20", "TensorPrismValues", "TensorPrism", "@21",
+ "ScalarPyramidValues", "ScalarPyramid", "@22", "VectorPyramidValues",
+ "VectorPyramid", "@23", "TensorPyramidValues", "TensorPyramid", "@24",
+ "Text2DValues", "Text2D", "@25", "Text3DValues", "Text3D", "@26",
+ "InterpolationMatrix", "NumericAffectation", "NumericIncrement",
+ "Affectation", "Shape", "Transform", "MultipleShape", "ListOfShapes",
+ "Duplicata", "Delete", "Colorify", "Visibility", "Command", "Loop",
+ "Extrude", "@27", "@28", "@29", "@30", "@31", "@32", "@33", "@34",
+ "@35", "ExtrudeParameters", "ExtrudeParameter", "Transfinite",
+ "Coherence", "BoolExpr", "FExpr", "FExpr_Single", "VExpr",
+ "VExpr_Single", "ListOfStrings", "RecursiveListOfStrings",
+ "ListOfListOfDouble", "RecursiveListOfListOfDouble", "ListOfDouble",
+ "FExpr_Multi", "RecursiveListOfDouble", "ColorExpr", "ListOfColor",
"RecursiveListOfColor", "StringExpr", 0
};
#endif
@@ -4983,8 +4977,7 @@ static const unsigned short yystos[] =
#define YYACCEPT goto yyacceptlab
#define YYABORT goto yyabortlab
-#define YYERROR goto yyerrorlab
-
+#define YYERROR goto yyerrlab1
/* Like YYERROR except do call yyerror. This remains here temporarily
to ease the transition to the new meaning of YYERROR, for GCC.
@@ -5018,11 +5011,11 @@ while (0)
are run). */
#ifndef YYLLOC_DEFAULT
-# define YYLLOC_DEFAULT(Current, Rhs, N) \
- ((Current).first_line = (Rhs)[1].first_line, \
- (Current).first_column = (Rhs)[1].first_column, \
- (Current).last_line = (Rhs)[N].last_line, \
- (Current).last_column = (Rhs)[N].last_column)
+# define YYLLOC_DEFAULT(Current, Rhs, N) \
+ Current.first_line = Rhs[1].first_line; \
+ Current.first_column = Rhs[1].first_column; \
+ Current.last_line = Rhs[N].last_line; \
+ Current.last_column = Rhs[N].last_column;
#endif
/* YYLEX -- calling `yylex' with the right arguments. */
@@ -5066,7 +5059,7 @@ do { \
/*------------------------------------------------------------------.
| yy_stack_print -- Print the state stack from its BOTTOM up to its |
-| TOP (included). |
+| TOP (cinluded). |
`------------------------------------------------------------------*/
#if defined (__STDC__) || defined (__cplusplus)
@@ -5106,9 +5099,9 @@ yy_reduce_print (yyrule)
#endif
{
int yyi;
- unsigned int yylno = yyrline[yyrule];
+ unsigned int yylineno = yyrline[yyrule];
YYFPRINTF (stderr, "Reducing stack by rule %d (line %u), ",
- yyrule - 1, yylno);
+ yyrule - 1, yylineno);
/* Print the symbols being reduced, and their result. */
for (yyi = yyprhs[yyrule]; 0 <= yyrhs[yyi]; yyi++)
YYFPRINTF (stderr, "%s ", yytname [yyrhs[yyi]]);
@@ -5145,7 +5138,7 @@ int yydebug;
SIZE_MAX < YYSTACK_BYTES (YYMAXDEPTH)
evaluated with infinite-precision integer arithmetic. */
-#if defined (YYMAXDEPTH) && YYMAXDEPTH == 0
+#if YYMAXDEPTH == 0
# undef YYMAXDEPTH
#endif
@@ -10303,8 +10296,8 @@ yyreduce:
}
-/* Line 1000 of yacc.c. */
-#line 10308 "Gmsh.tab.cpp"
+/* Line 991 of yacc.c. */
+#line 10300 "Gmsh.tab.cpp"
yyvsp -= yylen;
yyssp -= yylen;
@@ -10345,33 +10338,18 @@ yyerrlab:
{
YYSIZE_T yysize = 0;
int yytype = YYTRANSLATE (yychar);
- const char* yyprefix;
char *yymsg;
- int yyx;
+ int yyx, yycount;
+ yycount = 0;
/* Start YYX at -YYN if negative to avoid negative indexes in
YYCHECK. */
- int yyxbegin = yyn < 0 ? -yyn : 0;
-
- /* Stay within bounds of both yycheck and yytname. */
- int yychecklim = YYLAST - yyn;
- int yyxend = yychecklim < YYNTOKENS ? yychecklim : YYNTOKENS;
- int yycount = 0;
-
- yyprefix = ", expecting ";
- for (yyx = yyxbegin; yyx < yyxend; ++yyx)
+ for (yyx = yyn < 0 ? -yyn : 0;
+ yyx < (int) (sizeof (yytname) / sizeof (char *)); yyx++)
if (yycheck[yyx + yyn] == yyx && yyx != YYTERROR)
- {
- yysize += yystrlen (yyprefix) + yystrlen (yytname [yyx]);
- yycount += 1;
- if (yycount == 5)
- {
- yysize = 0;
- break;
- }
- }
- yysize += (sizeof ("syntax error, unexpected ")
- + yystrlen (yytname[yytype]));
+ yysize += yystrlen (yytname[yyx]) + 15, yycount++;
+ yysize += yystrlen ("syntax error, unexpected ") + 1;
+ yysize += yystrlen (yytname[yytype]);
yymsg = (char *) YYSTACK_ALLOC (yysize);
if (yymsg != 0)
{
@@ -10380,13 +10358,16 @@ yyerrlab:
if (yycount < 5)
{
- yyprefix = ", expecting ";
- for (yyx = yyxbegin; yyx < yyxend; ++yyx)
+ yycount = 0;
+ for (yyx = yyn < 0 ? -yyn : 0;
+ yyx < (int) (sizeof (yytname) / sizeof (char *));
+ yyx++)
if (yycheck[yyx + yyn] == yyx && yyx != YYTERROR)
{
- yyp = yystpcpy (yyp, yyprefix);
+ const char *yyq = ! yycount ? ", expecting " : " or ";
+ yyp = yystpcpy (yyp, yyq);
yyp = yystpcpy (yyp, yytname[yyx]);
- yyprefix = " or ";
+ yycount++;
}
}
yyerror (yymsg);
@@ -10407,56 +10388,52 @@ yyerrlab:
/* If just tried and failed to reuse lookahead token after an
error, discard it. */
- if (yychar <= YYEOF)
+ /* Return failure if at end of input. */
+ if (yychar == YYEOF)
{
- /* If at end of input, pop the error token,
- then the rest of the stack, then return failure. */
- if (yychar == YYEOF)
- for (;;)
- {
- YYPOPSTACK;
- if (yyssp == yyss)
- YYABORT;
- YYDSYMPRINTF ("Error: popping", yystos[*yyssp], yyvsp, yylsp);
- yydestruct (yystos[*yyssp], yyvsp);
- }
+ /* Pop the error token. */
+ YYPOPSTACK;
+ /* Pop the rest of the stack. */
+ while (yyss < yyssp)
+ {
+ YYDSYMPRINTF ("Error: popping", yystos[*yyssp], yyvsp, yylsp);
+ yydestruct (yystos[*yyssp], yyvsp);
+ YYPOPSTACK;
+ }
+ YYABORT;
}
- else
- {
- YYDSYMPRINTF ("Error: discarding", yytoken, &yylval, &yylloc);
- yydestruct (yytoken, &yylval);
- yychar = YYEMPTY;
- }
+ YYDSYMPRINTF ("Error: discarding", yytoken, &yylval, &yylloc);
+ yydestruct (yytoken, &yylval);
+ yychar = YYEMPTY;
+
}
/* Else will try to reuse lookahead token after shifting the error
token. */
- goto yyerrlab1;
+ goto yyerrlab2;
-/*---------------------------------------------------.
-| yyerrorlab -- error raised explicitly by YYERROR. |
-`---------------------------------------------------*/
-yyerrorlab:
+/*----------------------------------------------------.
+| yyerrlab1 -- error raised explicitly by an action. |
+`----------------------------------------------------*/
+yyerrlab1:
-#ifdef __GNUC__
- /* Pacify GCC when the user code never invokes YYERROR and the label
- yyerrorlab therefore never appears in user code. */
- if (0)
- goto yyerrorlab;
+ /* Suppress GCC warning that yyerrlab1 is unused when no action
+ invokes YYERROR. */
+#if defined (__GNUC_MINOR__) && 2093 <= (__GNUC__ * 1000 + __GNUC_MINOR__) \
+ && !defined __cplusplus
+ __attribute__ ((__unused__))
#endif
- yyvsp -= yylen;
- yyssp -= yylen;
- yystate = *yyssp;
- goto yyerrlab1;
+ goto yyerrlab2;
-/*-------------------------------------------------------------.
-| yyerrlab1 -- common code for both syntax error and YYERROR. |
-`-------------------------------------------------------------*/
-yyerrlab1:
+
+/*---------------------------------------------------------------.
+| yyerrlab2 -- pop states until the error token can be shifted. |
+`---------------------------------------------------------------*/
+yyerrlab2:
yyerrstatus = 3; /* Each real token shifted decrements this. */
for (;;)
@@ -10479,8 +10456,9 @@ yyerrlab1:
YYDSYMPRINTF ("Error: popping", yystos[*yyssp], yyvsp, yylsp);
yydestruct (yystos[yystate], yyvsp);
- YYPOPSTACK;
- yystate = *yyssp;
+ yyvsp--;
+ yystate = *--yyssp;
+
YY_STACK_PRINT (yyss, yyssp);
}
diff --git a/Parser/Gmsh.tab.hpp b/Parser/Gmsh.tab.hpp
index e8749a89361dea57af40778049b68131497ba2e3..835267587c9d03882ae1bf2974e0051f180f8bf5 100644
--- a/Parser/Gmsh.tab.hpp
+++ b/Parser/Gmsh.tab.hpp
@@ -1,7 +1,7 @@
-/* A Bison parser, made by GNU Bison 1.875c. */
+/* A Bison parser, made by GNU Bison 1.875. */
/* Skeleton parser for Yacc-like parsing with Bison,
- Copyright (C) 1984, 1989, 1990, 2000, 2001, 2002, 2003 Free Software Foundation, Inc.
+ Copyright (C) 1984, 1989, 1990, 2000, 2001, 2002 Free Software Foundation, Inc.
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
@@ -418,8 +418,8 @@ typedef union YYSTYPE {
Shape s;
List_T *l;
} YYSTYPE;
-/* Line 1275 of yacc.c. */
-#line 423 "Gmsh.tab.hpp"
+/* Line 1249 of yacc.c. */
+#line 422 "Gmsh.tab.hpp"
# define yystype YYSTYPE /* obsolescent; will be withdrawn */
# define YYSTYPE_IS_DECLARED 1
# define YYSTYPE_IS_TRIVIAL 1
diff --git a/Parser/Gmsh.yy.cpp b/Parser/Gmsh.yy.cpp
index d2fbf72142a2b577608cf60cf96e079edbf307ec..6b725f9d5ab532636e37ba4333e07fbd5c590fe6 100644
--- a/Parser/Gmsh.yy.cpp
+++ b/Parser/Gmsh.yy.cpp
@@ -2,7 +2,7 @@
/* A lexical scanner generated by flex */
/* Scanner skeleton version:
- * $Header: /cvsroot/gmsh/Parser/Gmsh.yy.cpp,v 1.209 2004-11-25 02:10:39 geuzaine Exp $
+ * $Header: /cvsroot/gmsh/Parser/Gmsh.yy.cpp,v 1.210 2004-11-25 16:22:47 remacle Exp $
*/
#define FLEX_SCANNER
@@ -10,7 +10,6 @@
#define YY_FLEX_MINOR_VERSION 5
#include <stdio.h>
-#include <unistd.h>
/* cfront 1.2 defines "c_plusplus" instead of "__cplusplus" */
@@ -24,6 +23,7 @@
#ifdef __cplusplus
#include <stdlib.h>
+#include <unistd.h>
/* Use prototypes in function declarations. */
#define YY_USE_PROTOS
@@ -1029,7 +1029,7 @@ char *yytext;
#line 1 "Gmsh.l"
#define INITIAL 0
#line 2 "Gmsh.l"
-// $Id: Gmsh.yy.cpp,v 1.209 2004-11-25 02:10:39 geuzaine Exp $
+// $Id: Gmsh.yy.cpp,v 1.210 2004-11-25 16:22:47 remacle Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -1238,7 +1238,7 @@ YY_MALLOC_DECL
YY_DECL
{
register yy_state_type yy_current_state;
- register char *yy_cp = NULL, *yy_bp = NULL;
+ register char *yy_cp, *yy_bp;
register int yy_act;
#line 80 "Gmsh.l"
@@ -2750,7 +2750,6 @@ register char *yy_bp;
#endif /* ifndef YY_NO_UNPUT */
-#ifndef YY_NO_INPUT
#ifdef __cplusplus
static int yyinput()
#else
@@ -2822,7 +2821,7 @@ static int input()
return c;
}
-#endif /* YY_NO_INPUT */
+
#ifdef YY_USE_PROTOS
void yyrestart( FILE *input_file )
@@ -2933,6 +2932,11 @@ YY_BUFFER_STATE b;
}
+#ifndef YY_ALWAYS_INTERACTIVE
+#ifndef YY_NEVER_INTERACTIVE
+extern int isatty YY_PROTO(( int ));
+#endif
+#endif
#ifdef YY_USE_PROTOS
void yy_init_buffer( YY_BUFFER_STATE b, FILE *file )
diff --git a/Plugin/Integrate.cpp b/Plugin/Integrate.cpp
index 4071f4d009db60fbd205ddb338623a4760768803..3de2cff56c631dd2cc0ee4d9cd9ae47ce82aa7df 100644
--- a/Plugin/Integrate.cpp
+++ b/Plugin/Integrate.cpp
@@ -1,4 +1,4 @@
-// $Id: Integrate.cpp,v 1.4 2004-11-25 02:10:40 geuzaine Exp $
+// $Id: Integrate.cpp,v 1.5 2004-11-25 16:22:48 remacle Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -93,6 +93,7 @@ static double integrate(int nbList, List_T *list, int dim,
const int levelsetPositive = (int)IntegrateOptions_Number[1].def;
double res = 0.;
+ double res2 = 0.;
int nb = List_Nbr(list) / nbList;
for(int i = 0; i < List_Nbr(list); i += nb) {
double *x = (double *)List_Pointer_Fast(list, i);
@@ -138,22 +139,38 @@ static double integrate(int nbList, List_T *list, int dim,
else if(nbComp == 3) res += q.integrateFlux(v);
}
}
- else if(dim == 3){
- if(nbNod == 4){
- tetrahedron t(x, y, z);
- if(nbComp == 1) res += t.integrate(v);
- }
- else if(nbNod == 8){
- hexahedron h(x, y, z);
- if(nbComp == 1) res += h.integrate(v);
- }
- else if(nbNod == 6){
- prism p(x, y, z);
- if(nbComp == 1) res += p.integrate(v);
- }
- else if(nbNod == 5){
- pyramid p(x, y, z);
- if(nbComp == 1) res += p.integrate(v);
+ else if(dim == 3)
+ {
+ if(nbNod == 4)
+ {
+ tetrahedron t(x, y, z);
+ if(nbComp == 1)
+ {
+ if ( ! levelsetPositive )
+ res += t.integrate(v);
+ else
+ {
+ double ONES[] = { 1.0 , 1.0 , 1.0, 1.0 };
+ const double area = fabs(t.integrate (ONES));
+ const double SUM = v[0] + v[1] + v[2] + v[3];
+ const double SUMABS = fabs(v[0]) + fabs(v[1]) + fabs (v[2]) + fabs (v[3]);
+ const double XI = SUM / SUMABS;
+ res += area * (1 - XI) * 0.5 ;
+ res2 += area * (1 + XI) * 0.5 ;
+ }
+ }
+ }
+ else if(nbNod == 8){
+ hexahedron h(x, y, z);
+ if(nbComp == 1) res += h.integrate(v);
+ }
+ else if(nbNod == 6){
+ prism p(x, y, z);
+ if(nbComp == 1) res += p.integrate(v);
+ }
+ else if(nbNod == 5){
+ pyramid p(x, y, z);
+ if(nbComp == 1) res += p.integrate(v);
}
}
}
diff --git a/Plugin/Levelset.cpp b/Plugin/Levelset.cpp
index f5d0f2ca0a061092387f440743b82b4169687f21..0a079c09fc614603e3ae3026a45ba6748efa8eae 100644
--- a/Plugin/Levelset.cpp
+++ b/Plugin/Levelset.cpp
@@ -1,4 +1,4 @@
-// $Id: Levelset.cpp,v 1.18 2004-11-25 02:10:40 geuzaine Exp $
+// $Id: Levelset.cpp,v 1.19 2004-11-25 16:22:48 remacle Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
@@ -343,6 +343,8 @@ Post_View *GMSH_LevelsetPlugin::execute(Post_View * v)
v->setAdaptiveResolutionLevel ( _recurLevel , this );
if (v->adaptive && v->NbSS)
v->setAdaptiveResolutionLevel ( _recurLevel , this );
+ if (v->adaptive && v->NbSQ)
+ v->setAdaptiveResolutionLevel ( _recurLevel , this );
if(_valueView < 0) {
@@ -621,6 +623,40 @@ static bool recur_sign_change (_tet *t, double val, const GMSH_LevelsetPlugin *p
}
}
+static bool recur_sign_change (_quad *q, double val, const GMSH_LevelsetPlugin *plug)
+{
+
+ if (!q->q[0])
+ {
+ double v1 = plug->levelset (q->p[0]->X,q->p[0]->Y,q->p[0]->Z,q->p[0]->val);
+ double v2 = plug->levelset (q->p[1]->X,q->p[1]->Y,q->p[1]->Z,q->p[1]->val);
+ double v3 = plug->levelset (q->p[2]->X,q->p[2]->Y,q->p[2]->Z,q->p[2]->val);
+ double v4 = plug->levelset (q->p[3]->X,q->p[3]->Y,q->p[3]->Z,q->p[3]->val);
+ if ( v1 * v2 > 0 && v1 * v3 > 0 && v1 * v4 > 0)
+ q->visible = false;
+ else
+ q->visible = true;
+ return q->visible;
+ }
+ else
+ {
+ bool sc1 = recur_sign_change(q->q[0],val,plug);
+ bool sc2 = recur_sign_change(q->q[1],val,plug);
+ bool sc3 = recur_sign_change(q->q[2],val,plug);
+ bool sc4 = recur_sign_change(q->q[3],val,plug);
+ if (sc1 || sc2 || sc3 || sc4 )
+ {
+ if (!sc1) q->q[0]->visible = true;
+ if (!sc2) q->q[1]->visible = true;
+ if (!sc3) q->q[2]->visible = true;
+ if (!sc4) q->q[3]->visible = true;
+ return true;
+ }
+ q->visible = false;
+ return false;
+ }
+}
+
void GMSH_LevelsetPlugin::assign_specific_visibility () const
{
if (_triangle::all_triangles.size())
@@ -633,4 +669,9 @@ void GMSH_LevelsetPlugin::assign_specific_visibility () const
_tet *te = *_tet::all_tets.begin();
te->visible = !recur_sign_change (te, _valueView, this);
}
+ if (_quad::all_quads.size())
+ {
+ _quad *qe = *_quad::all_quads.begin();
+ qe->visible = !recur_sign_change (qe, _valueView, this);
+ }
}
diff --git a/Plugin/Makefile b/Plugin/Makefile
index 6f9c86118b4f588f39e09db3497e8fd74c8d1c7f..13279d2a56a087ebe0fefddd3e488870cf04b491 100644
--- a/Plugin/Makefile
+++ b/Plugin/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.60 2004-11-13 22:52:46 geuzaine Exp $
+# $Id: Makefile,v 1.61 2004-11-25 16:22:48 remacle Exp $
#
# Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
#
@@ -29,7 +29,7 @@ CFLAGS = -DMPICH_SKIP_MPICXX ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = Plugin.cpp\
Levelset.cpp\
CutPlane.cpp CutSphere.cpp CutMap.cpp \
- Smooth.cpp CutParametric.cpp\
+ Smooth.cpp CutParametric.cpp Gradient.cpp Lambda2.cpp\
Octree.cpp OctreeInternals.cpp OctreePost.cpp\
StreamLines.cpp CutGrid.cpp\
Transform.cpp\