diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index d4a88c20a7deb376149c18d9d5e5e90604e736b7..c68e1d69f1eb40c38b8f21d654d0ef6ad1cf3b79 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.37 2007-09-04 13:47:01 remacle Exp $
+// $Id: MElement.cpp,v 1.38 2007-09-05 13:19:15 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -407,14 +407,14 @@ double MTriangle::getSurfaceXY() const
return s * 0.5;
}
-void MTriangle::circumcenterXYZ(double *p1, double *p2, double *p3,double *res)
+void MTriangle::circumcenterXYZ(double *p1, double *p2, double *p3,double *res, double *uv)
{
double v1[3] = {p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]};
double v2[3] = {p3[0]-p1[0],p3[1]-p1[1],p3[2]-p1[2]};
double vx[3] = {p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]};
double vy[3] = {p3[0]-p1[0],p3[1]-p1[1],p3[2]-p1[2]};
double vz[3]; prodve (vx,vy,vz);prodve (vz,vx,vy);
- norme(vx);norme(vy);norme(vz);
+ norme(vx); norme(vy);norme(vz);
double p1P[2] = {0.0,0.0};
double p2P[2];prosca(v1,vx,&p2P[0]);prosca(v1,vy,&p2P[1]);
double p3P[2];prosca(v2,vx,&p3P[0]);prosca(v2,vy,&p3P[1]);
@@ -422,6 +422,16 @@ void MTriangle::circumcenterXYZ(double *p1, double *p2, double *p3,double *res)
circumcenterXY(p1P, p2P, p3P,resP);
+ if (uv)
+ {
+ double mat[2][2] = {{p2P[0]-p1P[0],p3P[0]-p1P[0]},
+ {p2P[1]-p1P[1],p3P[1]-p1P[1]}};
+ double rhs[2] = {resP[0]-p1P[0],resP[1]-p1P[1]};
+ sys2x2(mat,rhs,uv);
+ }
+
+
+
// double d1 = sqrt((p2P[0] - resP[0]) * (p2P[0] - resP[0]) +
// (p2P[1] - resP[1]) * (p2P[1] - resP[1]));
@@ -437,24 +447,32 @@ void MTriangle::circumcenterXYZ(double *p1, double *p2, double *p3,double *res)
res[0] = p1[0] + resP[0] * vx[0] + resP[1] * vy[0];
res[1] = p1[1] + resP[0] * vx[1] + resP[1] * vy[1];
res[2] = p1[2] + resP[0] * vx[2] + resP[1] * vy[2];
+
+
+// double d2 = sqrt((p1P[0] - resP[0]) * (p1P[0] - resP[0]) +
+// (p1P[1] - resP[1]) * (p1P[1] - resP[1])) ;
+
+// double d3 = sqrt((p3P[0] - resP[0]) * (p3P[0] - resP[0]) +
+// (p3P[1] - resP[1]) * (p3P[1] - resP[1]) );
+
+
+
return;
- double d1 = sqrt((p1[0] - res[0]) * (p1[0] - res[0]) +
- (p1[1] - res[1]) * (p1[1] - res[1]) +
- (p1[2] - res[2]) * (p1[2] - res[2]) );
- double d2 = sqrt((p2[0] - res[0]) * (p2[0] - res[0]) +
- (p2[1] - res[1]) * (p2[1] - res[1]) +
- (p2[2] - res[2]) * (p2[2] - res[2]) );
- double d3 = sqrt((p3[0] - res[0]) * (p3[0] - res[0]) +
- (p3[1] - res[1]) * (p3[1] - res[1]) +
- (p3[2] - res[2]) * (p3[2] - res[2]) );
+ double d1 = sqrt((res[0] - p1[0]) * (res[0] - p1[0]) +
+ (res[1] - p1[1]) * (res[1] - p1[1]) +
+ (res[2] - p1[2]) * (res[2] - p1[2]) );
+ double d2 = sqrt((res[0] - p2[0]) * (res[0] - p2[0]) +
+ (res[1] - p2[1]) * (res[1] - p2[1]) +
+ (res[2] - p2[2]) * (res[2] - p2[2]) );
+ double d3 = sqrt((res[0] - p3[0]) * (res[0] - p3[0]) +
+ (res[1] - p3[1]) * (res[1] - p3[1]) +
+ (res[2] - p3[2]) * (res[2] - p3[2]) );
+ printf(" -- %g %g %g\n",d1,d2,d3);
- printf("%g %g %g\n",d1,d2,d3);
-
-
}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index b01ad669f321ec00e790cb69db6685e4749ce706..08d6d94a817f83bfd17290380505e61f8c4a4828 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -314,7 +314,7 @@ class MTriangle : public MElement {
}
void circumcenterXY(double *res) const;
void circumcenterUV(GFace*,double *res);
- static void circumcenterXYZ(double *p1, double *p2, double *p3,double *res);
+ static void circumcenterXYZ(double *p1, double *p2, double *p3,double *res, double *uv = 0);
static void circumcenterXY (double *p1, double *p2, double *p3,double *res);
double getSurfaceXY() const;
double getSurfaceUV(GFace*);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 70920753006c2acae71fd85adb6b5187f1e873ba..31061b49c66e0ad14d13304326d87b584826c245 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.85 2007-09-05 10:11:30 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.86 2007-09-05 13:19:15 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -40,6 +40,13 @@ extern Context_T CTX;
static double SCALINGU=1,SCALINGV=1;
+bool AlgoDelaunay2D ( GFace *gf )
+{
+ if ( gf->getNativeType() == GEntity::GmshModel && CTX.mesh.algo2d == ALGO_2D_DELAUNAY && gf->geomType() == GEntity::Plane)
+ return true;
+ return false;
+}
+
void computeEdgeLoops(const GFace *gf,
std::vector<MVertex*> &all_mvertices,
std::vector<int> &indices)
@@ -915,7 +922,7 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
m->del_point(m->find_point(-4));
// start mesh generation
- if (CTX.mesh.algo2d != ALGO_2D_DELAUNAY)
+ if (!AlgoDelaunay2D ( gf ))
{
RefineMesh (gf,*m,10);
OptimizeMesh(gf, *m, 2);
@@ -975,10 +982,10 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
// the delaunay algo is based directly on internal gmsh structures
// BDS mesh is passed in order not to recompute local coordinates
// of vertices
- if (CTX.mesh.algo2d == ALGO_2D_DELAUNAY)
+ if (AlgoDelaunay2D ( gf ))
{
insertVerticesInFace (gf,m) ;
- for (int i=0;i<5;i++)laplaceSmoothing (gf);
+ laplaceSmoothing (gf);
}
else if (debug){
char name[256];
@@ -1458,7 +1465,7 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
// goto hhh;
// start mesh generation
- if (CTX.mesh.algo2d != ALGO_2D_DELAUNAY)
+ if (!AlgoDelaunay2D ( gf ))
{
RefineMesh (gf,*m,10);
OptimizeMesh(gf, *m, 2);
@@ -1528,10 +1535,10 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
outputScalarField(m->triangles, name,1);
}
- if (CTX.mesh.algo2d == ALGO_2D_DELAUNAY)
+ if (AlgoDelaunay2D ( gf ))
{
insertVerticesInFace (gf,m) ;
- for (int i=0;i<5;i++)laplaceSmoothing (gf);
+ laplaceSmoothing (gf);
}
delete m;
@@ -1573,8 +1580,7 @@ void meshGFace::operator() (GFace *gf)
break;
case ALGO_2D_DELAUNAY:
// FIXME: Delaunay not available in all cases at the moment
- if(gf->geomType() == GEntity::Plane &&
- (gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty()))
+ if(AlgoDelaunay2D(gf))
algo = "Delaunay";
else
algo = "MeshAdapt+Delaunay";
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index afbbb17937fefd9d785bc278ebed3eafa7df762c..16e81e7d5f705fda5bc5f04e1d761ff4b6a5b6db 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceDelaunayInsertion.cpp,v 1.3 2007-09-04 13:47:02 remacle Exp $
+// $Id: meshGFaceDelaunayInsertion.cpp,v 1.4 2007-09-05 13:19:15 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -99,6 +99,34 @@ int MTri3::inCircumCircle ( const double *p ) const
return (result > 0) ? 1 : 0;
}
+int inCircumCircle ( MTriangle *base,
+ const double *p ,
+ const double *param ,
+ std::vector<double> & Us,
+ std::vector<double> & Vs)
+{
+ double pa[2] = {Us[base->getVertex(0)->getNum()],
+ Vs[base->getVertex(0)->getNum()]};
+ double pb[2] = {Us[base->getVertex(1)->getNum()],
+ Vs[base->getVertex(1)->getNum()]};
+ double pc[2] = {Us[base->getVertex(2)->getNum()],
+ Vs[base->getVertex(2)->getNum()]};
+
+ double result = gmsh::incircle(pa, pb, pc, (double*)param) * gmsh::orient2d(pa, pb, pc);
+ return (result > 0) ? 1 : 0;
+
+
+
+// double pa[3] = {base->getVertex(0)->x(),base->getVertex(0)->y(),base->getVertex(0)->z()};
+// double pb[3] = {base->getVertex(1)->x(),base->getVertex(1)->y(),base->getVertex(1)->z()};
+// double pc[3] = {base->getVertex(2)->x(),base->getVertex(2)->y(),base->getVertex(2)->z()};
+// double fourth[3];
+// fourthPoint ( pa,pb,pc,fourth);
+
+// double result = gmsh::insphere(pa, pb, pc, fourth, (double*)p) * gmsh::orient3d(pa, pb, pc,fourth);
+// return (result > 0) ? 1 : 0;
+}
+
template <class ITER>
void connectTris ( ITER beg, ITER end)
@@ -135,7 +163,10 @@ void connectTriangles ( std::list<MTri3*> & l)
void recurFindCavity (std::list<edgeXface> & shell,
std::list<MTri3*> & cavity,
double *v,
- MTri3 *t)
+ double *param,
+ MTri3 *t,
+ std::vector<double> & Us,
+ std::vector<double> & Vs)
{
t->setDeleted(true);
// the cavity that has to be removed
@@ -149,9 +180,9 @@ void recurFindCavity (std::list<edgeXface> & shell,
shell.push_back ( edgeXface ( t, i ) );
else if (!neigh->isDeleted())
{
- int circ = neigh->inCircumCircle ( v );
+ int circ = inCircumCircle (neigh->tri(), v , param, Us, Vs);
if (circ)
- recurFindCavity ( shell, cavity,v , neigh);
+ recurFindCavity ( shell, cavity,v, param, neigh,Us,Vs);
else
shell.push_back ( edgeXface ( t, i ) );
}
@@ -160,12 +191,25 @@ void recurFindCavity (std::list<edgeXface> & shell,
bool circUV ( MTriangle *t ,
std::vector<double> & Us,
- std::vector<double> & Vs , double *res)
+ std::vector<double> & Vs , double *res, GFace *gf)
{
- double p1 [2] = {Us[t->getVertex(0)->getNum()],Vs[t->getVertex(0)->getNum()]};
- double p2 [2] = {Us[t->getVertex(1)->getNum()],Vs[t->getVertex(1)->getNum()]};
- double p3 [2] = {Us[t->getVertex(2)->getNum()],Vs[t->getVertex(2)->getNum()]};
- t->circumcenterXY ( p1, p2 , p3, res);
+ double u1 [3] = {Us[t->getVertex(0)->getNum()],Vs[t->getVertex(0)->getNum()],0};
+ double u2 [3] = {Us[t->getVertex(1)->getNum()],Vs[t->getVertex(1)->getNum()],0};
+ double u3 [3] = {Us[t->getVertex(2)->getNum()],Vs[t->getVertex(2)->getNum()],0};
+ t->circumcenterXY ( u1, u2 , u3, res);
+ return true;
+ double p1 [3] = {t->getVertex(0)->x(),t->getVertex(0)->y(),t->getVertex(0)->z()};
+ double p2 [3] = {t->getVertex(1)->x(),t->getVertex(1)->y(),t->getVertex(1)->z()};
+ double p3 [3] = {t->getVertex(2)->x(),t->getVertex(2)->y(),t->getVertex(2)->z()};
+ double resxy[3], uv[2];
+ t->circumcenterXYZ ( p1, p2 , p3, resxy,uv);
+ return true;
+ // printf("uv = %g %g\n",uv[0],uv[1]);
+ // printf("%g %g vs ",res[0],res[1]);
+ res[0] = u1[0] + uv[0] * ( u2[0] - u1[0] ) + uv[1] * ( u3[0] - u1[0] );
+ res[1] = u1[1] + uv[0] * ( u2[1] - u1[1] ) + uv[1] * ( u3[1] - u1[1] );
+ // printf("%g %g \n",res[0],res[1]);
+
}
bool invMapUV ( MTriangle *t ,
@@ -221,6 +265,7 @@ double getSurfUV ( MTriangle *t ,
}
bool insertVertex (MVertex *v ,
+ double *param ,
MTri3 *t ,
std::set<MTri3*,compareTri3Ptr> &allTets,
std::vector<double> & vSizes,
@@ -234,7 +279,7 @@ bool insertVertex (MVertex *v ,
double p[3] = {v->x(),v->y(),v->z()};
- recurFindCavity ( shell, cavity, p , t);
+ recurFindCavity ( shell, cavity, p , param, t, Us, Vs);
// check that volume is conserved
double newVolume = 0;
@@ -251,8 +296,9 @@ bool insertVertex (MVertex *v ,
// fprintf(ff,"};\n");
// fclose (ff);
-// sprintf(name,"test%d.pos",III++);
-// ff = fopen (name,"w");
+// char name[256];
+// sprintf(name,"test.pos");
+// FILE *ff = fopen (name,"w");
// fprintf(ff,"View\"test\"{\n");
@@ -284,6 +330,16 @@ bool insertVertex (MVertex *v ,
// Us [t->getVertex(2)->getNum()],
// Vs [t->getVertex(2)->getNum()],
// 0.0);
+// fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0};\n",
+// t->getVertex(0)->x(),
+// t->getVertex(0)->y(),
+// t->getVertex(0)->z(),
+// t->getVertex(1)->x(),
+// t->getVertex(1)->y(),
+// t->getVertex(1)->z(),
+// t->getVertex(2)->x(),
+// t->getVertex(2)->y(),
+// t->getVertex(2)->z());
newTris[k++]=t4;
// all new triangles are pushed front in order to
// ba able to destroy them if the cavity is not
@@ -295,12 +351,12 @@ bool insertVertex (MVertex *v ,
new_cavity.push_back (otherSide);
// if (!it->t1->isDeleted())throw;
double ss = fabs(getSurfUV(t4->tri(),Us,Vs));
- if (ss < 1.e-12)ss = 1256172121;
+ if (ss < 1.e-25)ss = 1256172121;
newVolume += ss;
++it;
}
-// fprintf(ff,"};\n");
-// fclose (ff);
+// fprintf(ff,"};\n");
+// fclose (ff);
// printf("%d %d newVolume %g oldVolume %g\n",cavity.size(),new_cavity.size(),newVolume,oldVolume);
// getchar();
@@ -315,7 +371,7 @@ bool insertVertex (MVertex *v ,
// The cavity is NOT star shaped
else
{
- printf("the cavity (%d members) is not star shaped %g vs %g\n",shell.size(),oldVolume,newVolume);
+ // printf("the cavity (%d members) is not star shaped %g vs %g\n",shell.size(),oldVolume,newVolume);
for (unsigned int i=0;i<shell.size();i++)delete newTris[i];
delete [] newTris;
@@ -366,6 +422,33 @@ static void setLcs ( MTriangle *t, std::map<MVertex*,double> &vSizes)
}
}
+void _printTris (char *name, std::set<MTri3*,compareTri3Ptr>&AllTris, std::vector<double>&Us, std::vector<double>&Vs)
+{
+ FILE *ff = fopen (name,"w");
+ fprintf(ff,"View\"test\"{\n");
+ std::set<MTri3*,compareTri3Ptr>::iterator it = AllTris.begin();
+ while (it != AllTris.end() )
+ {
+ MTri3 *worst = *it;
+ if (!worst->isDeleted())
+ {
+ fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0};\n",
+ Us [(worst)->tri()->getVertex(0)->getNum()],
+ Vs [(worst)->tri()->getVertex(0)->getNum()],
+ 0.0,
+ Us [(worst)->tri()->getVertex(1)->getNum()],
+ Vs [(worst)->tri()->getVertex(1)->getNum()],
+ 0.0,
+ Us [(worst)->tri()->getVertex(2)->getNum()],
+ Vs [(worst)->tri()->getVertex(2)->getNum()],
+ 0.0);
+ }
+ ++it;
+ }
+ fprintf(ff,"};\n");
+ fclose (ff);
+}
+
void insertVerticesInFace (GFace *gf, BDS_Mesh *bds)
{
@@ -423,11 +506,17 @@ void insertVerticesInFace (GFace *gf, BDS_Mesh *bds)
Msg(DEBUG,"%7d points created -- Worst tri radius is %8.3f",vSizes.size(),worst->getRadius());
double center[2],uv[2];
if (worst->getRadius() < 0.5 * sqrt(2.0)) break;
- circUV(worst->tri(),Us,Vs,center);
+ circUV(worst->tri(),Us,Vs,center,gf);
bool inside = invMapUV(worst->tri(),center,Us,Vs,uv,1.e-8);
- if (1 || inside)
+ if (inside)
{
+
+// char name[245];
+// sprintf(name,"param%d.pos",ITER);
+// _printTris (name, AllTris, Us,Vs);
+
+
// we use here local coordinates as real coordinates
// x,y and z will be computed hereafter
// Msg(INFO,"Point is inside");
@@ -446,9 +535,9 @@ void insertVerticesInFace (GFace *gf, BDS_Mesh *bds)
Us.push_back( center[0] );
Vs.push_back( center[1] );
- if (!insertVertex ( v , worst, AllTris,vSizes,vSizesBGM,Us,Vs))
+ if (!insertVertex ( v , center, worst, AllTris,vSizes,vSizesBGM,Us,Vs))
{
- Msg(INFO,"cavity is invalid");
+ Msg(DEBUG,"2D Delaunay : a cavity is not star shaped");
AllTris.erase(AllTris.begin());
worst->forceRadius(-1);
AllTris.insert(worst);
@@ -459,18 +548,19 @@ void insertVerticesInFace (GFace *gf, BDS_Mesh *bds)
}
else
{
- Msg(INFO,"Point %g %g is outside %g %g - %g %g - %g %g (%g %g)",
- center[0],center[1],
- Us [(worst)->tri()->getVertex(0)->getNum()],
- Vs [(worst)->tri()->getVertex(0)->getNum()],
- Us [(worst)->tri()->getVertex(1)->getNum()],
- Vs [(worst)->tri()->getVertex(1)->getNum()],
- Us [(worst)->tri()->getVertex(2)->getNum()],
- Vs [(worst)->tri()->getVertex(2)->getNum()],
- uv[0],uv[1]);
+// Msg(DEBUG,"Point %g %g is outside %g %g - %g %g - %g %g (%g %g)",
+// center[0],center[1],
+// Us [(worst)->tri()->getVertex(0)->getNum()],
+// Vs [(worst)->tri()->getVertex(0)->getNum()],
+// Us [(worst)->tri()->getVertex(1)->getNum()],
+// Vs [(worst)->tri()->getVertex(1)->getNum()],
+// Us [(worst)->tri()->getVertex(2)->getNum()],
+// Vs [(worst)->tri()->getVertex(2)->getNum()],
+// uv[0],uv[1]);
AllTris.erase(AllTris.begin());
worst->forceRadius(0);
- AllTris.insert(worst);
+ AllTris.insert(worst);
+ // break;
}
}
}
@@ -478,10 +568,6 @@ void insertVerticesInFace (GFace *gf, BDS_Mesh *bds)
// optimize the mesh
// fill new gmsh structures with triangles
- char name[245];
- sprintf(name,"param%d.pos",gf->tag());
- FILE *ff = fopen (name,"w");
- fprintf(ff,"View\"test\"{\n");
while (1)
{
if (AllTris.begin() == AllTris.end() ) break;
@@ -493,20 +579,8 @@ void insertVerticesInFace (GFace *gf, BDS_Mesh *bds)
else
{
gf->triangles.push_back(worst->tri());
- fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0};\n",
- Us [(worst)->tri()->getVertex(0)->getNum()],
- Vs [(worst)->tri()->getVertex(0)->getNum()],
- 0.0,
- Us [(worst)->tri()->getVertex(1)->getNum()],
- Vs [(worst)->tri()->getVertex(1)->getNum()],
- 0.0,
- Us [(worst)->tri()->getVertex(2)->getNum()],
- Vs [(worst)->tri()->getVertex(2)->getNum()],
- 0.0);
}
delete worst;
AllTris.erase(AllTris.begin());
}
- fprintf(ff,"};\n");
- fclose (ff);
}
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 63beef62bc94496f379b5c2f7fbf8e9a11a6c012..01ddec5bbfc523b7176ebe9e83e72e59832d5376 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -78,41 +78,44 @@ void laplaceSmoothing (GFace *gf)
v2t_cont adj;
buildVertexToTriangle ( gf , adj );
- v2t_cont :: iterator it = adj.begin();
-
- while (it != adj.end())
+ for (int i=0;i<5;i++)
{
- MVertex *ver= it->first;
- GEntity *ge = ver->onWhat();
- // this vertex in internal to the face
- if (ge->dim() == 2)
+ v2t_cont :: iterator it = adj.begin();
+
+ while (it != adj.end())
{
- double initu,initv;
- ver->getParameter ( 0,initu);
- ver->getParameter ( 1,initv);
-
- const std::vector<MTriangle*> & lt = it->second;
-
- double fact = lt.size() ? 1./(3.*lt.size()):0;
-
- double cu=0,cv=0;
- double pu[3],pv[3];
- for (unsigned int i=0;i<lt.size();i++)
+ MVertex *ver= it->first;
+ GEntity *ge = ver->onWhat();
+ // this vertex in internal to the face
+ if (ge->dim() == 2)
{
- parametricCoordinates ( lt[i] , gf, pu, pv);
- cu += fact*(pu[0]+pu[1]+pu[2]);
- cv += fact*(pv[0]+pv[1]+pv[2]);
- // have to test validity !
+ double initu,initv;
+ ver->getParameter ( 0,initu);
+ ver->getParameter ( 1,initv);
+
+ const std::vector<MTriangle*> & lt = it->second;
+
+ double fact = lt.size() ? 1./(3.*lt.size()):0;
+
+ double cu=0,cv=0;
+ double pu[3],pv[3];
+ for (unsigned int i=0;i<lt.size();i++)
+ {
+ parametricCoordinates ( lt[i] , gf, pu, pv);
+ cu += fact*(pu[0]+pu[1]+pu[2]);
+ cv += fact*(pv[0]+pv[1]+pv[2]);
+ // have to test validity !
+ }
+ ver->setParameter(0,cu);
+ ver->setParameter(1,cv);
+ GPoint pt = gf->point(SPoint2(cu,cv));
+ ver->x() = pt.x();
+ ver->y() = pt.y();
+ ver->z() = pt.z();
}
- ver->setParameter(0,cu);
- ver->setParameter(1,cv);
- GPoint pt = gf->point(SPoint2(cu,cv));
- ver->x() = pt.x();
- ver->y() = pt.y();
- ver->z() = pt.z();
- }
- ++it;
- }
+ ++it;
+ }
+ }
}
extern void fourthPoint (double *p1, double *p2, double *p3, double *p4);