Skip to content
Snippets Groups Projects
Commit a9b14a1d authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

*** empty log message ***

parent 0809ecfa
No related branches found
No related tags found
No related merge requests found
......@@ -80,6 +80,19 @@ void _triangle::Recur_Create (_triangle *t, int maxlevel, int level , Double_Mat
all_triangles.push_back(t);
if (level++ >= maxlevel)
return;
/*
p3
p13 p23
p1 p12 p2
*/
_point *p1 = t->p[0];
_point *p2 = t->p[1];
......@@ -87,13 +100,13 @@ void _triangle::Recur_Create (_triangle *t, int maxlevel, int level , Double_Mat
_point *p12 = _point::New ( (p1->x+p2->x)*0.5,(p1->y+p2->y)*0.5,0, coeffs, eexps);
_point *p13 = _point::New ( (p1->x+p3->x)*0.5,(p1->y+p3->y)*0.5,0, coeffs, eexps);
_point *p23 = _point::New ( (p3->x+p2->x)*0.5,(p3->y+p2->y)*0.5,0, coeffs, eexps);
_triangle *t1 = new _triangle (p1,p12,p13);
_triangle *t1 = new _triangle (p1,p13,p12);
Recur_Create (t1, maxlevel,level,coeffs,eexps);
_triangle *t2 = new _triangle (p12,p23,p2);
Recur_Create (t2, maxlevel,level,coeffs,eexps);
_triangle *t3 = new _triangle (p23,p13,p3);
Recur_Create (t3, maxlevel,level,coeffs,eexps);
_triangle *t4 = new _triangle (p12,p23,p13);
_triangle *t4 = new _triangle (p12,p13,p23);
Recur_Create (t4, maxlevel,level,coeffs,eexps);
t->t[0]=t1;t->t[1]=t2;t->t[2]=t3;t->t[3]=t4;
......@@ -215,18 +228,15 @@ void Adaptive_Post_View:: zoomElement (Post_View * view ,
{
_point *p = (_point*) &(*it);
p->val = res(kk);
// for ( int k=0;k<N;++k)
// {
// p->val += p->shape_functions[k] * (*_STval )( ielem , k );
// }
// p->X = (*_STposX) ( ielem , 0 ) * ( 1.-p->x-p->y) + (*_STposX) ( ielem , 1 ) * p->x + (*_STposX) ( ielem , 2 ) * p->y;
// p->Y = (*_STposY) ( ielem , 0 ) * ( 1.-p->x-p->y) + (*_STposY) ( ielem , 1 ) * p->x + (*_STposY) ( ielem , 2 ) * p->y;
// p->Z = (*_STposZ) ( ielem , 0 ) * ( 1.-p->x-p->y) + (*_STposZ) ( ielem , 1 ) * p->x + (*_STposZ) ( ielem , 2 ) * p->y;
p->X = XYZ(kk,0);
p->Y = XYZ(kk,1);
p->Z = XYZ(kk,2);
/// -----------------------------------------------
/// TEST TEST ZZUT CHIOTTE VIREZ MOI CA VITE FAIT
/// -----------------------------------------------
// p->val = p->X * p->X + p->Y*p->Y - 1;
if (min > p->val) min = p->val;
if (max < p->val) max = p->val;
kk++;
......
// $Id: Integrate.cpp,v 1.1 2004-11-13 23:57:28 geuzaine Exp $
// $Id: Integrate.cpp,v 1.2 2004-11-22 11:35:26 remacle Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
......@@ -30,7 +30,8 @@
extern Context_T CTX;
StringXNumber IntegrateOptions_Number[] = {
{GMSH_FULLRC, "iView", NULL, -1.}
{GMSH_FULLRC, "iView", NULL, -1.},
{GMSH_FULLRC, "computeLevelsetPositive", NULL, 0.}
};
extern "C"
......@@ -63,7 +64,11 @@ void GMSH_IntegratePlugin::getInfos(char *author, char *copyright,
"line/surface elements. If `iView' < 0, the plugin\n"
"is run on the current view.\n"
"\n"
"Plugin(Integrate) creates one new view.\n");
"If computeLevelsetPositive is not 0, then the plugin \n"
"computes the positive area (volume)of the map\n"
"\n"
"Plugin(Integrate) creates one new view.\n"
);
}
int GMSH_IntegratePlugin::getNbOptions() const
......@@ -86,6 +91,8 @@ static double integrate(int nbList, List_T *list, int dim,
{
if(!nbList) return 0.;
const int levelsetPositive = (int)IntegrateOptions_Number[1].def;
double res = 0.;
int nb = List_Nbr(list) / nbList;
for(int i = 0; i < List_Nbr(list); i += nb) {
......@@ -110,7 +117,20 @@ static double integrate(int nbList, List_T *list, int dim,
else if(dim == 2){
if(nbNod == 3){
triangle t(x, y, z);
if(nbComp == 1) res += t.integrate(v);
if(nbComp == 1)
{
if ( ! levelsetPositive )
res += t.integrate(v);
else
{
double ONES[] = { 1.0 , 1.0 , 1.0 };
const double area = t.integrate (ONES);
const double SUM = v[0] + v[1] + v[2];
const double SUMABS = fabs(v[0]) + fabs(v[1]) + fabs (v[2]);
const double XI = SUM / SUMABS;
res += area * (1 - XI) * 0.5 ;
}
}
else if(nbComp == 3) res += t.integrateFlux(v);
}
else if(nbNod == 4){
......@@ -138,6 +158,7 @@ static double integrate(int nbList, List_T *list, int dim,
}
}
}
printf("integration res = %22.15E\n",res);
return res;
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment