Skip to content
Snippets Groups Projects
Commit c62553c0 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files 

cleaned up InterpolateSurface/Curve (get the negative curve dircetly in InterpolateCurve)

reintroduced fix for transfinite mesh on 3-sided transfinite surfaces: the problem was still there
parent 8077fd2e
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/GeoInterpolation.cpp
+ 45
− 134
View file @ c62553c0
// $Id: GeoInterpolation.cpp,v 1.11 2007-01-06 22:44:19 geuzaine Exp $
// $Id: GeoInterpolation.cpp,v 1.12 2007-01-07 10:52:46 geuzaine Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
......@@ -30,9 +30,14 @@
Vertex InterpolateCurve(Curve * c, double u, int derivee)
{
// we still call this in graphics
//if(c->Num < 0)
// Msg(GERROR, "FIXME: should never interpolate negative curves!");
if(c->Num < 0) {
Curve *C0 = FindCurve(-c->Num);
if(!C0){
Msg(GERROR, "Unknown curve %d", -c->Num);
return Vertex(0., 0., 0.);
}
return InterpolateCurve(C0, C0->ubeg + (C0->uend - C0->ubeg) * (1.-u), derivee);
}
Vertex V;
V.u = u;
......@@ -50,7 +55,7 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee)
int N, i, j;
Vertex *v[5];
double T[4], W, teta, t1, t2, t;
double T[4], W, theta, t1, t2, t;
double vec[4];
Vertex temp1, temp2;
......@@ -91,17 +96,14 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee)
V.u = 1. - u;
u = V.u;
}
teta = c->Circle.t1 - (c->Circle.t1 - c->Circle.t2) * u;
/* pour les ellipses */
teta -= c->Circle.incl;
theta = c->Circle.t1 - (c->Circle.t1 - c->Circle.t2) * u;
theta -= c->Circle.incl; // for ellipses
V.Pos.X =
c->Circle.f1 * cos(teta) * cos(c->Circle.incl) -
c->Circle.f2 * sin(teta) * sin(c->Circle.incl);
c->Circle.f1 * cos(theta) * cos(c->Circle.incl) -
c->Circle.f2 * sin(theta) * sin(c->Circle.incl);
V.Pos.Y =
c->Circle.f1 * cos(teta) * sin(c->Circle.incl) +
c->Circle.f2 * sin(teta) * cos(c->Circle.incl);
c->Circle.f1 * cos(theta) * sin(c->Circle.incl) +
c->Circle.f2 * sin(theta) * cos(c->Circle.incl);
V.Pos.Z = 0.0;
Projette(&V, c->Circle.invmat);
V.Pos.X += c->Circle.v[1]->Pos.X;
......@@ -113,10 +115,10 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee)
case MSH_SEGM_BSPLN:
case MSH_SEGM_BEZIER:
return InterpolateUBS(c, u, derivee);
return InterpolateUBS(c, u, 0);
case MSH_SEGM_NURBS:
return InterpolateNurbs(c, u, derivee);
return InterpolateNurbs(c, u, 0);
case MSH_SEGM_SPLN:
N = List_Nbr(c->Control_Points);
......@@ -174,18 +176,10 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee)
List_Read(c->Control_Points, i + 2, &v[3]);
}
if(derivee) {
T[3] = 0.;
T[2] = 1.;
T[1] = 2. * t;
T[0] = 3. * t * t;
}
else {
T[3] = 1.;
T[2] = t;
T[1] = t * t;
T[0] = t * t * t;
}
V.Pos.X = V.Pos.Y = V.Pos.Z = W = 0.0;
for(i = 0; i < 4; i++) {
......@@ -236,23 +230,10 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee)
for(j = 0; j < 4; j++) {
W += T[j] * vec[j];
}
if(derivee) {
V.Pos.X /= (t2 - t1);
V.Pos.Y /= (t2 - t1);
V.Pos.Z /= (t2 - t1);
}
else {
// V.Pos.X /= W;
// V.Pos.Y /= W;
// V.Pos.Z /= W;
}
return V;
default:
Msg(GERROR, "Unknown curve type in interpolation");
V.Pos.X = V.Pos.Y = V.Pos.Z = 0.0;
V.w = V.lc = 1.0;
return V;
}
......@@ -266,10 +247,6 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee)
#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)
#define DUTRAN_QUA(c1,c2,c3,c4,s1,s2,s3,s4,u,v) \
(-1.)*c4 +c2 -( (-1.)*(1.-v)*s1+ (1.-v)*s2+ v*s3 -v*s4)
#define DVTRAN_QUA(c1,c2,c3,c4,s1,s2,s3,s4,u,v) \
(-1.)*c1+ c3-( (-1.)*(1.-u)*s1 -u*s2+ u*s3+ (1.-u)*s4)
Vertex TransfiniteQua(Vertex c1, Vertex c2, Vertex c3, Vertex c4,
Vertex s1, Vertex s2, Vertex s3, Vertex s4,
......@@ -293,8 +270,6 @@ Vertex TransfiniteQua(Vertex c1, Vertex c2, Vertex c3, Vertex c4,
f(u,v) = u c2 (v) + (1-v) c1(u) + v c3(u) - u(1-v) s2 - uv s3 */
#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 DUTRAN_TRI(c1,c2,c3,s1,s2,s3,u,v) c2+ -( (1.-v)*s2+ v*s3);
#define DVTRAN_TRI(c1,c2,c3,s1,s2,s3,u,v) (-1.)*c1+ c3-(u*(-1) *s2+ u*s3);
Vertex TransfiniteTri(Vertex c1, Vertex c2, Vertex c3,
Vertex s1, Vertex s2, Vertex s3, double u, double v)
......@@ -331,8 +306,7 @@ void TransfiniteSph(Vertex S, Vertex center, Vertex * T)
T->Pos.Z = center.Pos.Z + r * dirz;
}
Vertex InterpolateRuledSurface(Surface * s, double u, double v,
int derivee, int u_v)
Vertex InterpolateRuledSurface(Surface * s, double u, double v)
{
Vertex *c1, *c2, T, D[4], V[4], *S[4];
Curve *C[4];
......@@ -365,32 +339,9 @@ Vertex InterpolateRuledSurface(Surface * s, double u, double v,
S[2] = C[2]->beg;
S[3] = C[3]->beg;
if (C[0]->Num < 0){
Curve *C0 = FindCurve(-C[0]->Num);
V[0] = InterpolateCurve(C0, C0->ubeg + (C0->uend - C0->ubeg) * (1.-u), 0);
}
else
V[0] = InterpolateCurve(C[0], C[0]->ubeg + (C[0]->uend - C[0]->ubeg) * u, 0);
if (C[1]->Num < 0){
Curve *C1 = FindCurve(-C[1]->Num);
V[1] = InterpolateCurve(C1, C1->ubeg + (C1->uend - C1->ubeg) * (1.-v), 0);
}
else
V[1] = InterpolateCurve(C[1], C[1]->ubeg + (C[1]->uend - C[1]->ubeg) * v, 0);
if (C[2]->Num < 0){
Curve *C2 = FindCurve(-C[2]->Num);
V[2] = InterpolateCurve(C2, C2->ubeg + (C2->uend - C2->ubeg) * u, 0);
}
else
V[2] = InterpolateCurve(C[2], C[2]->ubeg + (C[2]->uend - C[2]->ubeg) * (1. - u), 0);
if (C[3]->Num < 0){
Curve *C3 = FindCurve(-C[3]->Num);
V[3] = InterpolateCurve(C3, C3->ubeg + (C3->uend - C3->ubeg) * v, 0);
}
else
V[3] = InterpolateCurve(C[3], C[3]->ubeg + (C[3]->uend - C[3]->ubeg) * (1. - v), 0);
T = TransfiniteQua(V[0], V[1], V[2], V[3], *S[0], *S[1], *S[2], *S[3], u, v);
......@@ -421,25 +372,8 @@ Vertex InterpolateRuledSurface(Surface * s, double u, double v,
S[1] = C[1]->beg;
S[2] = C[2]->beg;
if (C[0]->Num < 0){
Curve *C0 = FindCurve(-C[0]->Num);
V[0] = InterpolateCurve(C0, C0->ubeg + (C0->uend - C0->ubeg) * (1.-u), 0);
}
else
V[0] = InterpolateCurve(C[0], C[0]->ubeg + (C[0]->uend - C[0]->ubeg) * u, 0);
if (C[1]->Num < 0){
Curve *C1 = FindCurve(-C[1]->Num);
V[1] = InterpolateCurve(C1, C1->ubeg + (C1->uend - C1->ubeg) * (1.-v), 0);
}
else
V[1] = InterpolateCurve(C[1], C[1]->ubeg + (C[1]->uend - C[1]->ubeg) * v, 0);
if (C[2]->Num < 0){
Curve *C2 = FindCurve(-C[2]->Num);
V[2] = InterpolateCurve(C2, C2->ubeg + (C2->uend - C2->ubeg) * u, 0);
}
else
V[2] = InterpolateCurve(C[2], C[2]->ubeg + (C[2]->uend - C[2]->ubeg) * (1.-u), 0);
T = TransfiniteTri(V[0], V[1], V[2], *S[0], *S[1], *S[2], u, v);
......@@ -449,17 +383,17 @@ Vertex InterpolateRuledSurface(Surface * s, double u, double v,
}
}
Vertex InterpolateExtrudedSurface(Surface * s, double u, double v,
int derivee, int u_v)
Vertex InterpolateExtrudedSurface(Surface * s, double u, double v)
{
Curve *c = FindCurve(s->Extrude->geo.Source);
// find position of c in the list of generatrices
Curve *C[4];
int num = -1;
for(int i = 0; i < List_Nbr(s->Generatrices); i++){
List_Read(s->Generatrices, i, &C[i]);
if(c == C[i]) num = i;
if(c == *(Curve**)List_Pointer(s->Generatrices, i)){
num = i;
break;
}
}
if(num < 0)
......@@ -469,38 +403,18 @@ Vertex InterpolateExtrudedSurface(Surface * s, double u, double v,
switch(num){
case 0:
if (c->Num < 0){
Curve *C0 = FindCurve(-c->Num);
T = InterpolateCurve(C0, C0->ubeg + (C0->uend - C0->ubeg) * (1. - u), 0);
}
else
T = InterpolateCurve(c, c->ubeg + (c->uend - c->ubeg) * u, 0);
s->Extrude->Extrude(v, T.Pos.X, T.Pos.Y, T.Pos.Z);
return T;
case 1:
if (c->Num < 0){
Curve *C0 = FindCurve(-c->Num);
T = InterpolateCurve(C0, C0->ubeg + (C0->uend - C0->ubeg) * (1. - v), 0);
}
else
T = InterpolateCurve(c, c->ubeg + (c->uend - c->ubeg) * v, 0);
s->Extrude->Extrude(1. - u, T.Pos.X, T.Pos.Y, T.Pos.Z);
return T;
case 2:
if (c->Num < 0){
Curve *C0 = FindCurve(-c->Num);
T = InterpolateCurve(C0, C0->ubeg + (C0->uend - C0->ubeg) * u, 0);
}
else
T = InterpolateCurve(c, c->ubeg + (c->uend - c->ubeg) * (1. - u), 0);
s->Extrude->Extrude(1. - v, T.Pos.X, T.Pos.Y, T.Pos.Z);
return T;
default:
if (c->Num < 0){
Curve *C0 = FindCurve(-c->Num);
T = InterpolateCurve(C0, C0->ubeg + (C0->uend - C0->ubeg) * v, 0);
}
else
T = InterpolateCurve(c, c->ubeg + (c->uend - c->ubeg) * (1. - v), 0);
s->Extrude->Extrude(u, T.Pos.X, T.Pos.Y, T.Pos.Z);
return T;
......@@ -522,7 +436,7 @@ Vertex InterpolateSurface(Surface * s, double u, double v, int derivee, int u_v)
D[1] = InterpolateSurface(s, u, v, 0, 0);
}
}
else if(u_v == 2) {
else {
if(v - eps < 0.0) {
D[0] = InterpolateSurface(s, u, v, 0, 0);
D[1] = InterpolateSurface(s, u, v + eps, 0, 0);
......@@ -532,9 +446,6 @@ Vertex InterpolateSurface(Surface * s, double u, double v, int derivee, int u_v)
D[1] = InterpolateSurface(s, u, v, 0, 0);
}
}
else {
Msg(WARNING, "Arbitrary InterpolateSurface for derivative not done");
}
T.Pos.X = (D[1].Pos.X - D[0].Pos.X) / eps;
T.Pos.Y = (D[1].Pos.Y - D[0].Pos.Y) / eps;
T.Pos.Z = (D[1].Pos.Z - D[0].Pos.Z) / eps;
......@@ -544,12 +455,12 @@ Vertex InterpolateSurface(Surface * s, double u, double v, int derivee, int u_v)
// use the exact extrusion formula if the surface is extruded
if(s->Extrude && s->Extrude->geo.Mode == EXTRUDED_ENTITY &&
s->Typ != MSH_SURF_PLAN)
return InterpolateExtrudedSurface(s, u, v, derivee, u_v);
return InterpolateExtrudedSurface(s, u, v);
switch (s->Typ) {
case MSH_SURF_REGL:
case MSH_SURF_TRIC:
return InterpolateRuledSurface(s, u, v, derivee, u_v);
return InterpolateRuledSurface(s, u, v);
case MSH_SURF_NURBS:
return InterpolateNurbsSurface(s, u, v);
case MSH_SURF_PLAN:
......
Mesh/meshGFaceTransfinite.cpp
+ 41
− 8
View file @ c62553c0
// $Id: meshGFaceTransfinite.cpp,v 1.15 2006-12-21 17:10:15 geuzaine Exp $
// $Id: meshGFaceTransfinite.cpp,v 1.16 2007-01-07 10:52:46 geuzaine Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
......@@ -29,17 +29,30 @@
#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)
/*
s4 +-----c3-----+ s3
| |
| |
c4 c2
| |
| |
s1 +-----c1-----+ s2
*/
// f(u,v) = (1-u) c4(v) + u c2(v) + (1-v) c1(u) + v c3(u)
// - [ (1-u)(1-v) s1 + u(1-v) s2 + uv s3 + (1-u)v s4 ]
#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)
// s1=s4=c4
// f(u,v) = u c2 (v) + (1-v) c1(u) + v c3(u) - u(1-v) s2 - uv s3
#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)
int MeshTransfiniteSurface(GFace *gf)
{
if(gf->meshAttributes.Method != TRANSFINI) return 0;
std::vector <MVertex *> corners;
std::vector <MVertex *> d_vertices;
std::vector <MVertex *> corners, d_vertices;
std::vector <int> indices;
for(unsigned int i = 0; i < gf->meshAttributes.corners.size(); i++)
......@@ -249,8 +262,28 @@ int MeshTransfiniteSurface( GFace *gf)
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);
double Up, Vp;
if(gf->geomType() != GEntity::RuledSurface){
Up = TRAN_TRI(U[iP1], U[iP2], U[iP3], UC1, UC2, UC3, u, v);
Vp = TRAN_TRI(V[iP1], V[iP2], V[iP3], VC1, VC2, VC3, u, v);
}
else{
// FIXME: to get nice meshes we would need to make the u,v
// coords match with the (degenerate) coordinates of the
// underlying ruled surface; so instead we just interpolate
// in real space
double xp = TRAN_TRI(m_vertices[iP1]->x(), m_vertices[iP2]->x(),
m_vertices[iP3]->x(), m_vertices[N1]->x(),
m_vertices[N2]->x(), m_vertices[N3]->x(), u, v);
double yp = TRAN_TRI(m_vertices[iP1]->y(), m_vertices[iP2]->y(),
m_vertices[iP3]->y(), m_vertices[N1]->y(),
m_vertices[N2]->y(), m_vertices[N3]->y(), u, v);
double zp = TRAN_TRI(m_vertices[iP1]->z(), m_vertices[iP2]->z(),
m_vertices[iP3]->z(), m_vertices[N1]->z(),
m_vertices[N2]->z(), m_vertices[N3]->z(), u, v);
// xp,yp,zp can be off the surface so we cannot use parFromPoint
gf->XYZtoUV(xp, yp, zp, Up, Vp, 1.0, false);
}
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);
......
benchmarks/2d/transfinite_tri_ugly_fix.geo
+ 3
− 3
View file @ c62553c0
......@@ -7,12 +7,12 @@ Line(1) = {2,3};
Line(2) = {2,1};
Line(3) = {1,3};
Line Loop(4) = {1,-3,-2};
Plane Surface(5) = {4};
Ruled Surface(5) = {4};
Transfinite Line{1,2,3} = 20;
//Transfinite Surface{5} = {3,1,2};
//Transfinite Surface{5} = {2,1,3};
Transfinite Surface{5} = {1,2,3};
Transfinite Surface{5} = {2,1,3};
//Transfinite Surface{5} = {1,2,3};
/*
Point(10)={3+0,0,0,1};
......
benchmarks/bugs/bug_ruled_new_algo.geo deleted 100644 → 0
+ 0
− 14
View file @ 8077fd2e
Point(2) = {0, 0, 1, 0.8};
Point(5) = {0, 0, 0, 0.8};
Point(6) = {0, 3, 0, 0.8};
Point(7) = {-4, 0, 0, 0.8};
Point(8) = {-4, 5, 0, 0.8};
Point(10) = {-4, 4, 0, 0.8};
Point(11) = {-4, 4, 1, 0.8};
Point(17) = {-4, 0, 1, 0.8};
Circle (3) = {6, 7, 8};
Ellipse (4) = {2, 5, 6, 6};
Circle (8) = {2, 17, 11};
Circle (20) = {8, 10, 11};
Line Loop (1) = {3, 20, -8, 4};
Ruled Surface (1) = {1};
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment