diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 5e3e02784d550c773db628896207971b64ee784b..f4cd1112e414c53ca15cbf92c2e1e43981ca7bd1 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.68 2006-11-27 22:35:01 geuzaine Exp $
+// $Id: BDS.cpp,v 1.69 2006-12-01 16:16:50 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -843,6 +843,9 @@ bool BDS_Mesh::recombine_edge(BDS_Edge * e)
BDS_Point *pts1[4];
e->faces(0)->getNodes(pts1);
+ // FIXME !!!!!!!!!!!!!!!!!
+ // should ensure that orientation is unchanged
+
BDS_Edge *p1_op1 = find_edge(p1, op[0], e->faces(0));
BDS_Edge *op1_p2 = find_edge(op[0], p2, e->faces(0));
BDS_Edge *p1_op2 = find_edge(p1, op[1], e->faces(1));
@@ -1015,20 +1018,22 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
double U = 0;
double V = 0;
+ // double tot_length = 0;
double LC = 0;
std::list < BDS_Edge * >::iterator eit = p->edges.begin();
while(eit != p->edges.end()) {
if ((*eit)->numfaces() == 1) return false;
BDS_Point *op = ((*eit)->p1 == p) ? (*eit)->p2 : (*eit)->p1;
- U += op->u;
+ U += op->u;
V += op->v;
+ // tot_length += (*eit)->length();
LC += op->lc();
++eit;
}
-
- U /= p->edges.size();
- V /= p->edges.size();
+
+ U /= (p->edges.size());
+ V /= (p->edges.size());
LC /= p->edges.size();
std::list < BDS_Face * >ts;
@@ -1050,8 +1055,6 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
p->X = gp.x();
p->Y = gp.y();
p->Z = gp.z();
- // p->radius() = gf->curvature (SPoint2(U,V));
-
eit = p->edges.begin();
while(eit != p->edges.end()) {
(*eit)->update();
diff --git a/Mesh/DivideAndConquer.cpp b/Mesh/DivideAndConquer.cpp
index e67e9744b4b0fa4cd9e109ce9cbf6329e4414990..d4d103347095ea029bc3034a01d3a8bd53b60ed9 100644
--- a/Mesh/DivideAndConquer.cpp
+++ b/Mesh/DivideAndConquer.cpp
@@ -1,4 +1,4 @@
-// $Id: DivideAndConquer.cpp,v 1.6 2006-11-30 13:55:20 geuzaine Exp $
+// $Id: DivideAndConquer.cpp,v 1.7 2006-12-01 16:16:50 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -533,9 +533,6 @@ PointNumero *ConvertDlistToArray(DListPeek * dlist, int *n)
void filldel(Delaunay * deladd, int aa, int bb, int cc,
PointRecord * points)
{
- double qual, newqual;
- MPoint pt2;
-
deladd->t.a = aa;
deladd->t.b = bb;
deladd->t.c = cc;
diff --git a/Mesh/meshGFaceTransfinite.cpp b/Mesh/meshGFaceTransfinite.cpp
index f5c517261991dc555f39dbc72de23bfaa3560722..820afe0ec19533bbe2dffdf6eedbc93c775d41d7 100644
--- a/Mesh/meshGFaceTransfinite.cpp
+++ b/Mesh/meshGFaceTransfinite.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceTransfinite.cpp,v 1.7 2006-11-27 22:22:17 geuzaine Exp $
+// $Id: meshGFaceTransfinite.cpp,v 1.8 2006-12-01 16:16:50 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -29,6 +29,8 @@
#include "Context.h"
#include "Message.h"
+#define TRAN_TRI(c1,c2,c3,s1,s2,s3,u,v) u*c2+(1.-v)*c1+v*c3-(u*(1.-v)*s2+u*v*s3);
+
#define TRAN_QUA(c1,c2,c3,c4,s1,s2,s3,s4,u,v) \
(1.-u)*c4+u*c2+(1.-v)*c1+v*c3-((1.-u)*(1.-v)*s1+u*(1.-v)*s2+u*v*s3+(1.-u)*v*s4)
@@ -62,7 +64,7 @@ int MeshTransfiniteSurface( GFace *gf)
m_vertices.push_back(d_vertices[(I + j) % d_vertices.size()]);
int iCorner = 0;
- int N[4];
+ int N[4] = {0,0,0,0};
std::vector<double> U;
std::vector<double> V;
std::map<std::pair<int,int> , MVertex*> tab;
@@ -90,15 +92,25 @@ int MeshTransfiniteSurface( GFace *gf)
int L = N2 - N1;
int H = N3 - N2;
- int Lb = N4 - N3;
- int Hb = m_vertices.size() - N4;
-
- // FIXME need to reimplement 3-side transfinite interpolation!
-
- if(Lb != L || Hb != H){
- Msg(GERROR,"Surface %d cannot be meshed using the transfinite algo", gf->tag());
- return 0;
- }
+ if (corners.size () == 4)
+ {
+ int Lb = N4 - N3;
+ int Hb = m_vertices.size() - N4;
+
+ if(Lb != L || Hb != H){
+ Msg(GERROR,"Surface %d cannot be meshed using the transfinite algo", gf->tag());
+ return 0;
+ }
+ }
+ else
+ {
+ int Lb = m_vertices.size() - N3;
+ // printf("%d %d %d %d\n",N1,N2,N3,N4);
+ if(Lb != L){
+ Msg(GERROR,"Surface %d cannot be meshed using the transfinite algo %d != %d", gf->tag(),L,Lb);
+ return 0;
+ }
+ }
std::vector<double> lengths_i;
std::vector<double> lengths_j;
@@ -127,7 +139,6 @@ int MeshTransfiniteSurface( GFace *gf)
lengths_j.push_back(L_j);
}
- //Msg(INFO,"L %d H %d -- %d -- %d %d %d %d", L, H, m_vertices.size(), N1, N2, N3, N4);
/*
2L+H
@@ -140,49 +151,91 @@ int MeshTransfiniteSurface( GFace *gf)
0 L
*/
- tab[std::make_pair(0,0)] = m_vertices[0];
- tab[std::make_pair(L,0)] = m_vertices[L];
- tab[std::make_pair(L,H)] = m_vertices[L+H];
- tab[std::make_pair(0,H)] = m_vertices[2*L+H];
-
- for (int i = 1; i < L; i++){
- tab[std::make_pair(i,0)] = m_vertices[i];
- tab[std::make_pair(i,H)] = m_vertices[2*L+H-i];
- }
- for(int i = 1; i < H;i++){
- tab[std::make_pair(L,i)] = m_vertices[L+i];
- tab[std::make_pair(0,i)] = m_vertices[2*L+2*H-i];
- }
+ if (corners.size () == 4)
+ {
+ tab[std::make_pair(0,0)] = m_vertices[0];
+ tab[std::make_pair(L,0)] = m_vertices[L];
+ tab[std::make_pair(L,H)] = m_vertices[L+H];
+ tab[std::make_pair(0,H)] = m_vertices[2*L+H];
+
+ for (int i = 1; i < L; i++){
+ tab[std::make_pair(i,0)] = m_vertices[i];
+ tab[std::make_pair(i,H)] = m_vertices[2*L+H-i];
+ }
+ for(int i = 1; i < H;i++){
+ tab[std::make_pair(L,i)] = m_vertices[L+i];
+ tab[std::make_pair(0,i)] = m_vertices[2*L+2*H-i];
+ }
+ }
+ else
+ {
+ tab[std::make_pair(0,0)] = m_vertices[0];
+ tab[std::make_pair(L,0)] = m_vertices[L];
+ tab[std::make_pair(L,H)] = m_vertices[L+H];
+
+ for (int i = 1; i < L; i++){
+ tab[std::make_pair(i,0)] = m_vertices[i];
+ tab[std::make_pair(i,H)] = m_vertices[2*L+H-i];
+ }
+ for(int i = 1; i < H;i++){
+ tab[std::make_pair(L,i)] = m_vertices[L+i];
+ }
+ }
double UC1 = U[N1];
double UC2 = U[N2];
double UC3 = U[N3];
- double UC4 = U[N4];
double VC1 = V[N1];
double VC2 = V[N2];
double VC3 = V[N3];
- double VC4 = V[N4];
+ Msg(INFO,"L %d H %d -- %d -- %d %d %d %d", L, H, m_vertices.size(), N1, N2, N3, N4);
//create points using transfinite interpolation
- for(int i = 1; i < L; i++){
- double u = lengths_i[i]/L_i;
- for(int j = 1; j < H; j++){
- double v = lengths_j[j]/L_j;
- int iP1 = N1+i;
- int iP2 = N2+j;
- int iP3 = N4-i;
- int iP4 = (N4+(N3-N2)-j)%m_vertices.size();
- double Up = TRAN_QUA ( U[iP1], U[iP2], U[iP3], U[iP4] , UC1, UC2, UC3, UC4, u, v );
- double Vp = TRAN_QUA ( V[iP1], V[iP2], V[iP3], V[iP4] , VC1, VC2, VC3, VC4, u, v );
- GPoint gp = gf->point (SPoint2(Up,Vp));
- MFaceVertex *newv = new MFaceVertex ( gp.x(),gp.y(),gp.z(), gf, Up, Vp );
- gf->mesh_vertices.push_back (newv);
- tab[std::make_pair(i,j)] = newv;
+ if (corners.size() == 4)
+ {
+ double UC4 = U[N4];
+ double VC4 = V[N4];
+ for(int i = 1; i < L; i++){
+ double u = lengths_i[i]/L_i;
+ for(int j = 1; j < H; j++){
+ double v = lengths_j[j]/L_j;
+ int iP1 = N1+i;
+ int iP2 = N2+j;
+ int iP3 = N4-i;
+ int iP4 = (N4+(N3-N2)-j)%m_vertices.size();
+ double Up = TRAN_QUA ( U[iP1], U[iP2], U[iP3], U[iP4] , UC1, UC2, UC3, UC4, u, v );
+ double Vp = TRAN_QUA ( V[iP1], V[iP2], V[iP3], V[iP4] , VC1, VC2, VC3, VC4, u, v );
+ GPoint gp = gf->point (SPoint2(Up,Vp));
+ MFaceVertex *newv = new MFaceVertex ( gp.x(),gp.y(),gp.z(), gf, Up, Vp );
+ gf->mesh_vertices.push_back (newv);
+ tab[std::make_pair(i,j)] = newv;
+ }
+ }
}
- }
-
+ else
+ {
+ for(int i = 1; i < L; i++){
+ double u = lengths_i[i]/L_i;
+ for(int j = 1; j < H; j++){
+ double v = lengths_j[j]/L_j;
+ int iP1 = N1+i;
+ int iP2 = N2+j;
+ int iP3 = ((N3+N2)-i)%m_vertices.size();
+
+
+ double Up = TRAN_TRI ( U[iP1], U[iP2], U[iP3] , UC1, UC2, UC3, u, v );
+ double Vp = TRAN_TRI ( V[iP1], V[iP2], V[iP3] , VC1, VC2, VC3, u, v );
+
+ GPoint gp = gf->point (SPoint2(Up,Vp));
+ MFaceVertex *newv = new MFaceVertex ( gp.x(),gp.y(),gp.z(), gf, Up, Vp );
+ gf->mesh_vertices.push_back (newv);
+ tab[std::make_pair(i,j)] = newv;
+ }
+ }
+ }
+
// elliptic smoother (don't apply this by default)
- if(CTX.mesh.nb_smoothing > 1 && gf->geomType() == GEntity::Plane){
+ if(corners.size() == 4 && CTX.mesh.nb_smoothing > 1 && gf->geomType() == GEntity::Plane){
for (int IT = 0; IT< CTX.mesh.nb_smoothing; IT++){
for(int i = 1; i < L; i++){
for(int j = 1; j < H; j++){
@@ -221,28 +274,62 @@ int MeshTransfiniteSurface( GFace *gf)
}
}
}
-
- // create elements
- for(int i = 0; i < L ; i++){
- for(int j = 0; j < H; j++){
- MVertex *v1 = tab[std::make_pair(i,j)];
- MVertex *v2 = tab[std::make_pair(i+1,j)];
- MVertex *v3 = tab[std::make_pair(i+1,j+1)];
- MVertex *v4 = tab[std::make_pair(i,j+1)];
- if(gf->meshAttributes.recombine)
- gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
- else if(gf->meshAttributes.transfiniteArrangement == 1 ||
- (gf->meshAttributes.transfiniteArrangement == 0 &&
- ((i % 2 == 0 && j % 2 == 1) ||
- (i % 2 == 1 && j % 2 == 0)))){
+
+ if (corners.size() == 4)
+ {
+ // create elements
+ for(int i = 0; i < L ; i++){
+ for(int j = 0; j < H; j++){
+ MVertex *v1 = tab[std::make_pair(i,j)];
+ MVertex *v2 = tab[std::make_pair(i+1,j)];
+ MVertex *v3 = tab[std::make_pair(i+1,j+1)];
+ MVertex *v4 = tab[std::make_pair(i,j+1)];
+ if(gf->meshAttributes.recombine)
+ gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
+ else if(gf->meshAttributes.transfiniteArrangement == 1 ||
+ (gf->meshAttributes.transfiniteArrangement == 0 &&
+ ((i % 2 == 0 && j % 2 == 1) ||
+ (i % 2 == 1 && j % 2 == 0)))){
+ gf->triangles.push_back(new MTriangle(v1, v2, v3));
+ gf->triangles.push_back(new MTriangle(v3, v4, v1));
+ }
+ else{
+ gf->triangles.push_back(new MTriangle(v1, v2, v4));
+ gf->triangles.push_back(new MTriangle(v4, v2, v3));
+ }
+ }
+ }
+ }
+ else
+ {
+ for(int j = 0; j < H; j++){
+ MVertex *v1 = tab[std::make_pair(0,0)];
+ MVertex *v2 = tab[std::make_pair(1,j)];
+ MVertex *v3 = tab[std::make_pair(1,j+1)];
gf->triangles.push_back(new MTriangle(v1, v2, v3));
- gf->triangles.push_back(new MTriangle(v3, v4, v1));
- }
- else{
- gf->triangles.push_back(new MTriangle(v1, v2, v4));
- gf->triangles.push_back(new MTriangle(v4, v2, v3));
- }
+ }
+
+ for(int i = 1; i < L ; i++){
+ for(int j = 0; j < H; j++){
+ MVertex *v1 = tab[std::make_pair(i,j)];
+ MVertex *v2 = tab[std::make_pair(i+1,j)];
+ MVertex *v3 = tab[std::make_pair(i+1,j+1)];
+ MVertex *v4 = tab[std::make_pair(i,j+1)];
+ if(gf->meshAttributes.recombine)
+ gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
+ else if(gf->meshAttributes.transfiniteArrangement == 1 ||
+ (gf->meshAttributes.transfiniteArrangement == 0 &&
+ ((i % 2 == 0 && j % 2 == 1) ||
+ (i % 2 == 1 && j % 2 == 0)))){
+ gf->triangles.push_back(new MTriangle(v1, v2, v3));
+ gf->triangles.push_back(new MTriangle(v3, v4, v1));
+ }
+ else{
+ gf->triangles.push_back(new MTriangle(v1, v2, v4));
+ gf->triangles.push_back(new MTriangle(v4, v2, v3));
+ }
+ }
+ }
}
- }
return 1;
}
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index a7443aca00d0d327e4087e0f496fd94ae4521fec..07d1204a156360b5c15ee2e8dfeaadb240435c3a 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionDelaunayInsertion.cpp,v 1.8 2006-11-30 11:32:26 remacle Exp $
+// $Id: meshGRegionDelaunayInsertion.cpp,v 1.9 2006-12-01 16:16:50 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -284,6 +284,18 @@ static void setLcs ( MTetrahedron *t, std::map<MVertex*,double> &vSizes)
}
}
+// void recover_volumes( GRegion *gr , std::set<MTet4*,compareTet4Ptr> & allTets )
+// {
+// std::set<MTet4*,compareTet4Ptr>::iterator it = allTets.begin();
+// for ( ; it != allTets.end() ; ++it )
+// {
+// MTet4 *t = *allTets.begin();
+// if (!t->isDeleted())
+// {
+// }
+// }
+// }
+
void insertVerticesInRegion (GRegion *gr)
{
@@ -366,7 +378,6 @@ void insertVerticesInRegion (GRegion *gr)
}
}
}
-
while (1)
{
diff --git a/doc/TODO b/doc/TODO
index c944acc0b1ae56a620bc4740f9126213bc9c5d0c..466cf4817a8c120c60c09ce89f46f3b8a5b1d4d5 100644
--- a/doc/TODO
+++ b/doc/TODO
@@ -1,9 +1,12 @@
-$Id: TODO,v 1.29 2006-11-30 15:14:29 geuzaine Exp $
+$Id: TODO,v 1.30 2006-12-01 16:16:50 remacle Exp $
********************************************************************
test isSeam() dans le second ordre -- si oui, interpo lineaireemnt
-
+transfinite 3 curves
+subdivision of extrusion
+bug : ruled surface mesh generation
+reclassify volumes in new 3D mesh
********************************************************************