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

Transfinite interpolation with triangular patches done

parent 0257c4e6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Mesh/BDS.cpp
+ 10
− 7
View file @ d1aa0471
// $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();
......
Mesh/DivideAndConquer.cpp
+ 1
− 4
View file @ d1aa0471
// $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;
......
Mesh/meshGFaceTransfinite.cpp
+ 152
− 65
View file @ d1aa0471
// $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;
}
Mesh/meshGRegionDelaunayInsertion.cpp
+ 13
− 2
View file @ d1aa0471
// $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)
{
......
doc/TODO
+ 5
− 2
View file @ d1aa0471
$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
********************************************************************
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment