From 153881abdcbfae51efb54bf47db9dca51f3d29e9 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Mon, 27 Jun 2005 15:03:46 +0000
Subject: [PATCH] *** empty log message ***
---
Common/GmshMatrix.h | 29 +
Fltk/Callbacks.cpp | 16 +-
Fltk/Callbacks.h | 4 +
Fltk/GUI.cpp | 52 +-
Graphics/Geom.cpp | 1234 +++++++++++++++++++-------------------
Mesh/BDS.cpp | 407 ++++++++++---
Mesh/BDS.h | 97 +--
Mesh/DiscreteSurface.cpp | 24 +-
Parser/OpenFile.cpp | 5 +-
9 files changed, 1094 insertions(+), 774 deletions(-)
diff --git a/Common/GmshMatrix.h b/Common/GmshMatrix.h
index f8dcf55bf0..4c02c183ec 100644
--- a/Common/GmshMatrix.h
+++ b/Common/GmshMatrix.h
@@ -1,6 +1,8 @@
#ifndef _GMSH_BOOSTMATRIX_
#define _GMSH_BOOSTMATRIX_
+#include <assert.h>
+
template <class SCALAR>
class Gmsh_Vector
{
@@ -42,6 +44,10 @@ public:
{
for (int i=0;i<r;++i)data[i]=other.data[i];
}
+ inline void lu_solve (const Gmsh_Vector<SCALAR>& rhs, Gmsh_Vector<SCALAR> & result)
+ {
+ throw;
+ }
};
template <class SCALAR>
@@ -97,6 +103,7 @@ public:
};
#ifdef HAVE_GSL
+#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
class GSL_Vector
@@ -161,6 +168,28 @@ public:
gsl_blas_dgemm (CblasNoTrans,CblasNoTrans, 1.0, data, x.data, 1.0, b.data);
}
+ inline void least_squares (const GSL_Vector & rhs, GSL_Vector & result)
+ {
+ assert (r > c);
+ assert (rhs.size() == r);
+ assert (result.size() == c);
+ GSL_Matrix *ls = new GSL_Matrix(c, c);
+ GSL_Vector *ls_rhs = new GSL_Vector(c);
+ gsl_blas_dgemm (CblasTrans,CblasNoTrans, 1.0, data, data, 1.0, ls->data);
+ gsl_blas_dgemv (CblasTrans, 1.0, data, rhs.data, 1.0, ls_rhs->data);
+ ls->lu_solve (*ls_rhs,result);
+ delete ls;
+ delete ls_rhs;
+ }
+ inline void lu_solve (const GSL_Vector & rhs, GSL_Vector & result)
+ {
+ int s;
+ gsl_permutation * p = gsl_permutation_alloc (size1());
+ gsl_linalg_LU_decomp ( data, p, &s);
+ gsl_linalg_LU_solve ( data , p, rhs.data, result.data ) ;
+ gsl_permutation_free (p);
+ }
+
inline void mult (const GSL_Vector & x, GSL_Vector & b )
{
gsl_blas_dgemv (CblasNoTrans, 1.0, data, x.data, 1.0, b.data);
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index 328dd2fc6e..8f7f6b7331 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.356 2005-06-10 20:59:14 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.357 2005-06-27 15:03:45 remacle Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -25,6 +25,7 @@
#include <signal.h>
#include <map>
+#include "BDS.h"
#include "Gmsh.h"
#include "GmshUI.h"
#include "Geo.h"
@@ -829,6 +830,19 @@ void options_restore_defaults_cb(CALLBACK_ARGS)
Draw();
}
+void wizard_update_edges_cb(CALLBACK_ARGS)
+{
+ extern void BDS_To_Mesh(Mesh *m);
+ if (THEM && THEM->bds && WID)
+ {
+ const double angle = WID->swiz_value[0]->value() * M_PI / 180;
+ printf ("value = %g\n",angle);
+ THEM->bds->classify (angle);
+ BDS_To_Mesh (THEM);
+ Draw();
+ }
+}
+
void options_ok_cb(CALLBACK_ARGS)
{
general_options_ok_cb(w, data);
diff --git a/Fltk/Callbacks.h b/Fltk/Callbacks.h
index 15a4a52603..6d4b6fe31d 100644
--- a/Fltk/Callbacks.h
+++ b/Fltk/Callbacks.h
@@ -291,5 +291,9 @@ void solver_kill_cb(CALLBACK_ARGS);
void solver_choose_executable_cb(CALLBACK_ARGS);
void solver_ok_cb(CALLBACK_ARGS);
+// SURFACE MESH WIZARD
+
+void wizard_update_edges_cb(CALLBACK_ARGS);
+
#endif
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index 3c0caf0b73..7dd1f5e0fe 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.443 2005-06-09 22:28:24 geuzaine Exp $
+// $Id: GUI.cpp,v 1.444 2005-06-27 15:03:45 remacle Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -812,6 +812,7 @@ GUI::GUI(int argc, char **argv)
int i;
// initialize static windows
+ swiz_window = NULL;
m_window = NULL;
g_window = NULL;
opt_window = NULL;
@@ -901,6 +902,9 @@ GUI::GUI(int argc, char **argv)
}
call_for_solver_plugin(-1);
+ // create the surface mesh wizard
+ create_surface_mesh_wizard();
+
// Draw the scene
g_opengl_window->redraw();
}
@@ -1612,7 +1616,6 @@ void GUI::create_surface_mesh_wizard()
{
int width = 42 * fontsize;
int height = 12 * BH + 5 * WB;
- int L = 105 + WB;
if(swiz_window) {
swiz_window->show();
@@ -1621,16 +1624,47 @@ void GUI::create_surface_mesh_wizard()
swiz_window = new Dialog_Window(width, height);
swiz_window->box(GMSH_WINDOW_BOX);
- swiz_wiz = new Fl_Wizard(L, 0, width, height, "Surface Mesh Wizard");
+ swiz_wiz = new Fl_Wizard(0, 0, width, height, "Surface Mesh Wizard");
{
- Fl_Tabs *o = new Fl_Tabs(L + WB, WB, width - 2 * WB, height - 2 * WB);
+ Fl_Tabs *o = new Fl_Tabs(WB, WB, width - 2 * WB, height - 2 * WB);
{
- Fl_Group *o = new Fl_Group(L + WB, WB + BH, width - 2 * WB, height - 2 * WB - BH, "Detect Edges");
- swiz_value[0] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 1 * BH, IW, BH, "Angle");
- swiz_value[0]->minimum(0.1);
- swiz_value[0]->maximum(50);
- swiz_value[0]->step(0.1);
+ Fl_Group *o = new Fl_Group(WB, WB + BH, width - 2 * WB, height - 2 * WB - BH, "Model Edge Detection");
+ swiz_value[0] = new Fl_Value_Input(2 * WB, 2 * WB + 1 * BH, IW, BH, "Treshold Dihedral Angle");
+ swiz_value[0]->value(180/8);
+ swiz_value[0]->minimum(0);
+ swiz_value[0]->maximum(90);
+ swiz_value[0]->step(1);
swiz_value[0]->align(FL_ALIGN_RIGHT);
+ {
+ Fl_Return_Button *b = new Fl_Return_Button(width - 3 * BB - 3 * WB, height - BH - WB, BB, BH, "Apply");
+ b->callback(wizard_update_edges_cb);
+ }
+ {
+ Fl_Button *b = new Fl_Button(width - 2 * BB - 2 * WB, height - BH - WB, BB, BH, "Next");
+// b->callback(options_save_cb);
+ }
+ {
+ Fl_Button *b = new Fl_Button(width - BB - WB, height - BH - WB, BB, BH, "Cancel");
+ b->callback(cancel_cb, (void *)swiz_window);
+ }
+
+ o->end();
+ }
+ {
+ Fl_Group *o = new Fl_Group(WB, WB + BH, width - 2 * WB, height - 2 * WB - BH, "Reverse Engineering the CAD");
+ {
+ Fl_Return_Button *b = new Fl_Return_Button(width - 3 * BB - 3 * WB, height - BH - WB, BB, BH, "Apply");
+// b->callback(wizard_update_edges_cb);
+ }
+ {
+ Fl_Button *b = new Fl_Button(width - 2 * BB - 2 * WB, height - BH - WB, BB, BH, "Next");
+// b->callback(options_save_cb);
+ }
+ {
+ Fl_Button *b = new Fl_Button(width - BB - WB, height - BH - WB, BB, BH, "Cancel");
+ b->callback(cancel_cb, (void *)swiz_window);
+ }
+
o->end();
}
o->end();
diff --git a/Graphics/Geom.cpp b/Graphics/Geom.cpp
index 6830c1fa52..9fcf1c378d 100644
--- a/Graphics/Geom.cpp
+++ b/Graphics/Geom.cpp
@@ -1,4 +1,4 @@
-// $Id: Geom.cpp,v 1.84 2005-05-16 00:50:10 geuzaine Exp $
+// $Id: Geom.cpp,v 1.85 2005-06-27 15:03:45 remacle Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -41,304 +41,299 @@ extern Mesh *THEM;
void Draw_Geo_Point(void *a, void *b)
{
- char Num[100];
-
- Vertex *v = *(Vertex **) a;
-
- if(!(v->Visible & VIS_GEOM))
- return;
-
- if(CTX.render_mode == GMSH_SELECT) {
- glLoadName(0);
- glPushName(v->Num);
- }
-
- if(v->Frozen) {
- glPointSize(CTX.geom.point_sel_size);
- gl2psPointSize(CTX.geom.point_sel_size * CTX.print.eps_point_size_factor);
- glColor4ubv((GLubyte *) & CTX.color.geom.point_sel);
- }
- else {
- glPointSize(CTX.geom.point_size);
- gl2psPointSize(CTX.geom.point_size * CTX.print.eps_point_size_factor);
- glColor4ubv((GLubyte *) & CTX.color.geom.point);
- }
-
- if(CTX.geom.points) {
- if(CTX.geom.point_type == 1) {
- if(v->Frozen)
- Draw_Sphere(CTX.geom.point_sel_size, v->Pos.X, v->Pos.Y, v->Pos.Z,
- CTX.geom.light);
- else
- Draw_Sphere(CTX.geom.point_size, v->Pos.X, v->Pos.Y, v->Pos.Z,
- CTX.geom.light);
+ char Num[100];
+
+ Vertex *v = *(Vertex **) a;
+
+ if(!(v->Visible & VIS_GEOM))
+ return;
+
+ if(CTX.render_mode == GMSH_SELECT) {
+ glLoadName(0);
+ glPushName(v->Num);
}
- else if(CTX.geom.point_type == 2) {
- GMSH_Solve_Plugin *sp = GMSH_PluginManager::instance()->findSolverPlugin();
- if(sp) {
- sp-> GL_enhancePoint (v);
- }
- glBegin(GL_POINTS);
- glVertex3d(v->Pos.X, v->Pos.Y, v->Pos.Z);
- glEnd();
+
+ if(v->Frozen) {
+ glPointSize(CTX.geom.point_sel_size);
+ gl2psPointSize(CTX.geom.point_sel_size * CTX.print.eps_point_size_factor);
+ glColor4ubv((GLubyte *) & CTX.color.geom.point_sel);
}
else {
- glBegin(GL_POINTS);
- glVertex3d(v->Pos.X, v->Pos.Y, v->Pos.Z);
- glEnd();
+ glPointSize(CTX.geom.point_size);
+ gl2psPointSize(CTX.geom.point_size * CTX.print.eps_point_size_factor);
+ glColor4ubv((GLubyte *) & CTX.color.geom.point);
}
- }
+ if(CTX.geom.points) {
+ if(CTX.geom.point_type == 1) {
+ if(v->Frozen)
+ Draw_Sphere(CTX.geom.point_sel_size, v->Pos.X, v->Pos.Y, v->Pos.Z,
+ CTX.geom.light);
+ else
+ Draw_Sphere(CTX.geom.point_size, v->Pos.X, v->Pos.Y, v->Pos.Z,
+ CTX.geom.light);
+ }
+ else if(CTX.geom.point_type == 2) {
+ GMSH_Solve_Plugin *sp = GMSH_PluginManager::instance()->findSolverPlugin();
+ if(sp) {
+ sp-> GL_enhancePoint (v);
+ }
+ glBegin(GL_POINTS);
+ glVertex3d(v->Pos.X, v->Pos.Y, v->Pos.Z);
+ glEnd();
+ }
+ else {
+ glBegin(GL_POINTS);
+ glVertex3d(v->Pos.X, v->Pos.Y, v->Pos.Z);
+ glEnd();
+ }
+
+ }
- if(CTX.geom.points_num) {
- sprintf(Num, "%d", v->Num);
- double offset = (0.5 * CTX.geom.point_size + 0.3 * CTX.gl_fontsize) * CTX.pixel_equiv_x;
- glRasterPos3d(v->Pos.X + offset / CTX.s[0],
- v->Pos.Y + offset / CTX.s[1],
- v->Pos.Z + offset / CTX.s[2]);
- Draw_String(Num);
- }
+ if(CTX.geom.points_num) {
+ sprintf(Num, "%d", v->Num);
+ double offset = (0.5 * CTX.geom.point_size + 0.3 * CTX.gl_fontsize) * CTX.pixel_equiv_x;
+ glRasterPos3d(v->Pos.X + offset / CTX.s[0],
+ v->Pos.Y + offset / CTX.s[1],
+ v->Pos.Z + offset / CTX.s[2]);
+ Draw_String(Num);
+ }
- if(CTX.render_mode == GMSH_SELECT) {
- glPopName();
- }
+ if(CTX.render_mode == GMSH_SELECT) {
+ glPopName();
+ }
}
// Curves
void Draw_Curve(void *a, void *b)
{
- int N;
- double mod, x[2], y[2], z[2];
- char Num[100];
- Vertex v, dv;
-
- Curve *c = *(Curve **) a;
-
- if(c->Num < 0 || !(c->Visible & VIS_GEOM))
- return;
-
- if(CTX.render_mode == GMSH_SELECT) {
- glLoadName(1);
- glPushName(c->Num);
- }
-
- if(c->ipar[3]) {
- glLineWidth(CTX.geom.line_sel_width);
- gl2psLineWidth(CTX.geom.line_sel_width * CTX.print.eps_line_width_factor);
- glColor4ubv((GLubyte *) & CTX.color.geom.line_sel);
- }
- else {
- glLineWidth(CTX.geom.line_width);
- gl2psLineWidth(CTX.geom.line_width * CTX.print.eps_line_width_factor);
- glColor4ubv((GLubyte *) & CTX.color.geom.line);
- }
-
- if(CTX.geom.lines) {
- int n = List_Nbr(c->Control_Points);
- switch (c->Typ) {
- case MSH_SEGM_LINE:
- N = n;
- break;
- case MSH_SEGM_CIRC:
- case MSH_SEGM_CIRC_INV:
- case MSH_SEGM_ELLI:
- case MSH_SEGM_ELLI_INV:
- N = CTX.geom.circle_points;
- break;
- default:
- N = 10 * n;
- break;
- }
- if(c->bds) {
- std::set<BDS_Edge*, EdgeLessThan>::iterator it = c->bds->edges.begin();
- std::set<BDS_Edge*, EdgeLessThan>::iterator ite = c->bds->edges.end();
- while (it!=ite)
- {
- if ((*it)->g && (*it)->g->classif_degree == 1&& c->Num == (*it)->g->classif_tag)
- {
- BDS_Edge *e = (*it);
- glBegin(GL_LINES);
- glVertex3d(e->p1->X,e->p1->Y,e->p1->Z);
- glVertex3d(e->p2->X,e->p2->Y,e->p2->Z);
- glEnd();
- }
- ++it;
- }
- }
- else if(c->Typ == MSH_SEGM_DISCRETE) {
- // do nothing
+ int N;
+ double mod, x[2], y[2], z[2];
+ char Num[100];
+ Vertex v, dv;
+
+ Curve *c = *(Curve **) a;
+
+ if(c->Num < 0 || !(c->Visible & VIS_GEOM))
+ return;
+
+ if(CTX.render_mode == GMSH_SELECT) {
+ glLoadName(1);
+ glPushName(c->Num);
+ }
+
+ if(c->ipar[3]) {
+ glLineWidth(CTX.geom.line_sel_width);
+ gl2psLineWidth(CTX.geom.line_sel_width * CTX.print.eps_line_width_factor);
+ glColor4ubv((GLubyte *) & CTX.color.geom.line_sel);
}
else {
- if(CTX.geom.line_type >= 1) {
- for(int i = 0; i < N - 1; i++) {
- v = InterpolateCurve(c, (double)i / (double)(N - 1), 0);
- dv = InterpolateCurve(c, (double)(i + 1) / (double)(N - 1), 0);
- x[0] = v.Pos.X;
- y[0] = v.Pos.Y;
- z[0] = v.Pos.Z;
- x[1] = dv.Pos.X;
- y[1] = dv.Pos.Y;
- z[1] = dv.Pos.Z;
- Draw_Cylinder(c->ipar[3] ? CTX.geom.line_sel_width : CTX.geom.line_width,
- x, y, z, CTX.geom.light);
- }
- if(CTX.geom.line_type == 2) {
- GMSH_Solve_Plugin *sp = GMSH_PluginManager::instance()->findSolverPlugin();
- if(sp) {
- int NN=(N>1)?N:1;
- const double eps=0.e-2;
- for(int i = 0; i < NN - 1; i++) {
- v = InterpolateCurve(c, (double)i / (double)(NN - 1)-eps, 0);
- dv = InterpolateCurve(c, (double)(i + 1) / (double)(NN - 1)+eps, 0);
- sp-> GL_enhanceLine (c->Num,&v,&dv);
+ glLineWidth(CTX.geom.line_width);
+ gl2psLineWidth(CTX.geom.line_width * CTX.print.eps_line_width_factor);
+ glColor4ubv((GLubyte *) & CTX.color.geom.line);
+ }
+
+ if(CTX.geom.lines) {
+ int n = List_Nbr(c->Control_Points);
+ switch (c->Typ) {
+ case MSH_SEGM_LINE:
+ N = n;
+ break;
+ case MSH_SEGM_CIRC:
+ case MSH_SEGM_CIRC_INV:
+ case MSH_SEGM_ELLI:
+ case MSH_SEGM_ELLI_INV:
+ N = CTX.geom.circle_points;
+ break;
+ default:
+ N = 10 * n;
+ break;
+ }
+ if(c->Typ == MSH_SEGM_DISCRETE && c->bds) {
+
+ BDS_GeomEntity *g = c->bds->get_geom ( c->Num,1);
+ std::list<BDS_Edge*>::iterator it = g->e.begin();
+ std::list<BDS_Edge*>::iterator ite = g->e.end();
+ while (it!=ite)
+ {
+ BDS_Edge *e = (*it);
+ glBegin(GL_LINES);
+ glVertex3d(e->p1->X,e->p1->Y,e->p1->Z);
+ glVertex3d(e->p2->X,e->p2->Y,e->p2->Z);
+ glEnd();
+ ++it;
+ }
+ }
+ else {
+ if(CTX.geom.line_type >= 1) {
+ for(int i = 0; i < N - 1; i++) {
+ v = InterpolateCurve(c, (double)i / (double)(N - 1), 0);
+ dv = InterpolateCurve(c, (double)(i + 1) / (double)(N - 1), 0);
+ x[0] = v.Pos.X;
+ y[0] = v.Pos.Y;
+ z[0] = v.Pos.Z;
+ x[1] = dv.Pos.X;
+ y[1] = dv.Pos.Y;
+ z[1] = dv.Pos.Z;
+ Draw_Cylinder(c->ipar[3] ? CTX.geom.line_sel_width : CTX.geom.line_width,
+ x, y, z, CTX.geom.light);
+ }
+ if(CTX.geom.line_type == 2) {
+ GMSH_Solve_Plugin *sp = GMSH_PluginManager::instance()->findSolverPlugin();
+ if(sp) {
+ int NN=(N>1)?N:1;
+ const double eps=0.e-2;
+ for(int i = 0; i < NN - 1; i++) {
+ v = InterpolateCurve(c, (double)i / (double)(NN - 1)-eps, 0);
+ dv = InterpolateCurve(c, (double)(i + 1) / (double)(NN - 1)+eps, 0);
+ sp-> GL_enhanceLine (c->Num,&v,&dv);
+ }
+ }
+ }
+ }
+ else {
+ glBegin(GL_LINE_STRIP);
+ for(int i = 0; i < N; i++) {
+ v = InterpolateCurve(c, (double)i / (double)(N - 1), 0);
+ glVertex3d(v.Pos.X, v.Pos.Y, v.Pos.Z);
+ }
+ glEnd();
}
- }
}
- }
- else {
- glBegin(GL_LINE_STRIP);
- for(int i = 0; i < N; i++) {
- v = InterpolateCurve(c, (double)i / (double)(N - 1), 0);
- glVertex3d(v.Pos.X, v.Pos.Y, v.Pos.Z);
- }
- glEnd();
- }
- }
- }
-
- if(CTX.geom.lines_num) {
- v = InterpolateCurve(c, 0.5, 0);
- sprintf(Num, "%d", c->Num);
- double offset = (0.5 * CTX.geom.line_width + 0.3 * CTX.gl_fontsize) * CTX.pixel_equiv_x;
- glRasterPos3d(v.Pos.X + offset / CTX.s[0],
- v.Pos.Y + offset / CTX.s[1],
- v.Pos.Z + offset / CTX.s[2]);
- Draw_String(Num);
- }
-
- if(CTX.geom.tangents) {
- v = InterpolateCurve(c, 0.5, 0);
- dv = InterpolateCurve(c, 0.5, 1);
- mod = sqrt(dv.Pos.X * dv.Pos.X + dv.Pos.Y * dv.Pos.Y + dv.Pos.Z * dv.Pos.Z);
- dv.Pos.X = dv.Pos.X / mod * CTX.geom.tangents * CTX.pixel_equiv_x / CTX.s[0];
- dv.Pos.Y = dv.Pos.Y / mod * CTX.geom.tangents * CTX.pixel_equiv_x / CTX.s[1];
- dv.Pos.Z = dv.Pos.Z / mod * CTX.geom.tangents * CTX.pixel_equiv_x / CTX.s[2];
- glColor4ubv((GLubyte *) & CTX.color.geom.tangents);
- Draw_Vector(CTX.vector_type, 0, CTX.arrow_rel_head_radius,
- CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius,
- v.Pos.X, v.Pos.Y, v.Pos.Z,
- dv.Pos.X, dv.Pos.Y, dv.Pos.Z, CTX.geom.light);
- }
-
- if(CTX.render_mode == GMSH_SELECT) {
- glPopName();
- }
+ }
+
+ if(CTX.geom.lines_num) {
+ v = InterpolateCurve(c, 0.5, 0);
+ sprintf(Num, "%d", c->Num);
+ double offset = (0.5 * CTX.geom.line_width + 0.3 * CTX.gl_fontsize) * CTX.pixel_equiv_x;
+ glRasterPos3d(v.Pos.X + offset / CTX.s[0],
+ v.Pos.Y + offset / CTX.s[1],
+ v.Pos.Z + offset / CTX.s[2]);
+ Draw_String(Num);
+ }
+
+ if(CTX.geom.tangents) {
+ v = InterpolateCurve(c, 0.5, 0);
+ dv = InterpolateCurve(c, 0.5, 1);
+ mod = sqrt(dv.Pos.X * dv.Pos.X + dv.Pos.Y * dv.Pos.Y + dv.Pos.Z * dv.Pos.Z);
+ dv.Pos.X = dv.Pos.X / mod * CTX.geom.tangents * CTX.pixel_equiv_x / CTX.s[0];
+ dv.Pos.Y = dv.Pos.Y / mod * CTX.geom.tangents * CTX.pixel_equiv_x / CTX.s[1];
+ dv.Pos.Z = dv.Pos.Z / mod * CTX.geom.tangents * CTX.pixel_equiv_x / CTX.s[2];
+ glColor4ubv((GLubyte *) & CTX.color.geom.tangents);
+ Draw_Vector(CTX.vector_type, 0, CTX.arrow_rel_head_radius,
+ CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius,
+ v.Pos.X, v.Pos.Y, v.Pos.Z,
+ dv.Pos.X, dv.Pos.Y, dv.Pos.Z, CTX.geom.light);
+ }
+
+ if(CTX.render_mode == GMSH_SELECT) {
+ glPopName();
+ }
}
// Surfaces
void putZ(Vertex * v, Surface * s)
{
- Vertex V;
- V.Pos.X = s->a;
- V.Pos.Y = s->b;
- V.Pos.Z = s->c;
- Projette(&V, s->plan);
- if(V.Pos.Z != 0.0)
- v->Pos.Z = (s->d - V.Pos.X * v->Pos.X - V.Pos.Y * v->Pos.Y) / V.Pos.Z;
- else
- v->Pos.Z = 0.0;
-
- Projette(v, s->invplan);
+ Vertex V;
+ V.Pos.X = s->a;
+ V.Pos.Y = s->b;
+ V.Pos.Z = s->c;
+ Projette(&V, s->plan);
+ if(V.Pos.Z != 0.0)
+ v->Pos.Z = (s->d - V.Pos.X * v->Pos.X - V.Pos.Y * v->Pos.Y) / V.Pos.Z;
+ else
+ v->Pos.Z = 0.0;
+
+ Projette(v, s->invplan);
}
int isPointOnPlanarSurface(Surface * S, double X, double Y, double Z,
double n[3])
{
- double Angle = 0.0;
- Vertex V;
- V.Pos.X = X;
- V.Pos.Y = Y;
- V.Pos.Z = Z;
+ double Angle = 0.0;
+ Vertex V;
+ V.Pos.X = X;
+ V.Pos.Y = Y;
+ V.Pos.Z = Z;
- for(int i = 0; i < List_Nbr(S->Generatrices); i++) {
+ for(int i = 0; i < List_Nbr(S->Generatrices); i++) {
- Curve *C;
- List_Read(S->Generatrices, i, &C);
+ Curve *C;
+ List_Read(S->Generatrices, i, &C);
- int N = (C->Typ == MSH_SEGM_LINE) ? 1 : 10;
+ int N = (C->Typ == MSH_SEGM_LINE) ? 1 : 10;
- for(int j = 0; j < N; j++) {
- double u1 = (double)j / (double)(N);
- double u2 = (double)(j + 1) / (double)(N);
- Vertex P1 = InterpolateCurve(C, u1, 0);
- Vertex P2 = InterpolateCurve(C, u2, 0);
- Angle += angle_plan(&V, &P1, &P2, n);
- }
+ for(int j = 0; j < N; j++) {
+ double u1 = (double)j / (double)(N);
+ double u2 = (double)(j + 1) / (double)(N);
+ Vertex P1 = InterpolateCurve(C, u1, 0);
+ Vertex P2 = InterpolateCurve(C, u2, 0);
+ Angle += angle_plan(&V, &P1, &P2, n);
+ }
- }
+ }
- // if Angle == 2*pi, we're inside
+ // if Angle == 2*pi, we're inside
- if(fabs(Angle) > 2*M_PI-0.5 && fabs(Angle) < 2*M_PI+0.5)
- return 1;
- return 0;
+ if(fabs(Angle) > 2*M_PI-0.5 && fabs(Angle) < 2*M_PI+0.5)
+ return 1;
+ return 0;
}
void getPlaneSurfaceNormal(Surface *s, double x, double y, double z, double n[3])
{
- double t1[3], t2[3];
- Curve *c;
-
- for(int i = 0; i < 3; i++)
- n[i] = 0.0;
-
- if(s->Typ != MSH_SURF_PLAN)
- return;
-
- // We cannot use s->plan or InterpolateSurface here (they both rely
- // on the mean plane, which borks the orientation). So we need to
- // use this trick:
-
- if(List_Nbr(s->Generatrices) > 1){
- List_Read(s->Generatrices, 0, &c);
- if(c->beg){
- t1[0] = x - c->beg->Pos.X;
- t1[1] = y - c->beg->Pos.Y;
- t1[2] = z - c->beg->Pos.Z;
- for(int i = 1; i < List_Nbr(s->Generatrices); i++){
- List_Read(s->Generatrices, i, &c);
- t2[0] = x - c->beg->Pos.X;
- t2[1] = y - c->beg->Pos.Y;
- t2[2] = z - c->beg->Pos.Z;
- prodve(t1, t2, n);
- if(norme(n))
- break;
- }
- }
- }
- if(List_Nbr(s->Generatrices) < 2 || !norme(n)){
- Msg(WARNING, "Reverting to mean plane normal for surface %d", s->Num);
- n[0] = s->a;
- n[1] = s->b;
- n[2] = s->c;
- norme(n);
- }
+ double t1[3], t2[3];
+ Curve *c;
+
+ for(int i = 0; i < 3; i++)
+ n[i] = 0.0;
+
+ if(s->Typ != MSH_SURF_PLAN)
+ return;
+
+ // We cannot use s->plan or InterpolateSurface here (they both rely
+ // on the mean plane, which borks the orientation). So we need to
+ // use this trick:
+
+ if(List_Nbr(s->Generatrices) > 1){
+ List_Read(s->Generatrices, 0, &c);
+ if(c->beg){
+ t1[0] = x - c->beg->Pos.X;
+ t1[1] = y - c->beg->Pos.Y;
+ t1[2] = z - c->beg->Pos.Z;
+ for(int i = 1; i < List_Nbr(s->Generatrices); i++){
+ List_Read(s->Generatrices, i, &c);
+ t2[0] = x - c->beg->Pos.X;
+ t2[1] = y - c->beg->Pos.Y;
+ t2[2] = z - c->beg->Pos.Z;
+ prodve(t1, t2, n);
+ if(norme(n))
+ break;
+ }
+ }
+ }
+ if(List_Nbr(s->Generatrices) < 2 || !norme(n)){
+ Msg(WARNING, "Reverting to mean plane normal for surface %d", s->Num);
+ n[0] = s->a;
+ n[1] = s->b;
+ n[2] = s->c;
+ norme(n);
+ }
}
void Draw_Polygonal_Surface(Surface * s)
{
- if(CTX.geom.surfaces) {
- if(CTX.geom.light) glEnable(GL_LIGHTING);
- glEnable(GL_POLYGON_OFFSET_FILL); // always!
-
- std::list<BDS_Triangle*>::iterator it = s->bds->triangles.begin();
- std::list<BDS_Triangle*>::iterator ite = s->bds->triangles.end();
- while (it!=ite)
- {
- if (!(*it)->g || s->Num == (*it)->g->classif_tag)
- {
+ if(CTX.geom.surfaces) {
+ if(CTX.geom.light) glEnable(GL_LIGHTING);
+ glEnable(GL_POLYGON_OFFSET_FILL); // always!
+
+ BDS_GeomEntity *g = s->bds->get_geom ( s->Num,2);
+ std::list<BDS_Triangle*>::iterator it = g->t.begin();
+ std::list<BDS_Triangle*>::iterator ite = g->t.end();
+ while (it!=ite)
+ {
glBegin(GL_TRIANGLES);
BDS_Point *n[3];
BDS_Triangle *t = (*it);
@@ -352,290 +347,289 @@ void Draw_Polygonal_Surface(Surface * s)
if(CTX.geom.light) glNormal3dv(c);
glVertex3d(n[2]->X,n[2]->Y,n[2]->Z);
glEnd();
+ ++it;
}
- ++it;
- }
- glDisable(GL_POLYGON_OFFSET_FILL);
- glDisable(GL_LIGHTING);
- }
+ glDisable(GL_POLYGON_OFFSET_FILL);
+ glDisable(GL_LIGHTING);
+ }
}
void Draw_Plane_Surface(Surface * s)
{
- Curve *c;
- Vertex V[4], vv, vv1, vv2;
- double n[3];
+ Curve *c;
+ Vertex V[4], vv, vv1, vv2;
+ double n[3];
- if(!CTX.threads_lock && List_Nbr(s->Orientations) < 1) {
- CTX.threads_lock = 1;
+ if(!CTX.threads_lock && List_Nbr(s->Orientations) < 1) {
+ CTX.threads_lock = 1;
- List_T *points = List_Create(10, 10, sizeof(Vertex *));
- for(int i = 0; i < List_Nbr(s->Generatrices); i++) {
- List_Read(s->Generatrices, i, &c);
- for(int j = 0; j < List_Nbr(c->Control_Points); j++) {
- List_Add(points, List_Pointer(c->Control_Points, j));
- }
- }
+ List_T *points = List_Create(10, 10, sizeof(Vertex *));
+ for(int i = 0; i < List_Nbr(s->Generatrices); i++) {
+ List_Read(s->Generatrices, i, &c);
+ for(int j = 0; j < List_Nbr(c->Control_Points); j++) {
+ List_Add(points, List_Pointer(c->Control_Points, j));
+ }
+ }
- if(!List_Nbr(points)){
- List_Delete(points);
- CTX.threads_lock = 0;
- return;
- }
- MeanPlane(points, s);
- List_Delete(points);
+ if(!List_Nbr(points)){
+ List_Delete(points);
+ CTX.threads_lock = 0;
+ return;
+ }
+ MeanPlane(points, s);
+ List_Delete(points);
- // compute (rough) bounding box of surface
- double minx = 0., miny = 0., maxx = 0., maxy = 0.;
- for(int i = 0; i < List_Nbr(s->Generatrices); i++) {
- List_Read(s->Generatrices, i, &c);
- Vertex P1 = InterpolateCurve(c, 0.0, 0);
- Vertex P2 = InterpolateCurve(c, 0.5, 0);
- Vertex P3 = InterpolateCurve(c, 1.0, 0);
- Projette(&P1, s->plan);
- Projette(&P2, s->plan);
- Projette(&P3, s->plan);
- if(!i) {
- minx = maxx = P1.Pos.X;
- miny = maxy = P1.Pos.Y;
- }
- minx = DMIN(DMIN(DMIN(minx, P1.Pos.X), P2.Pos.X), P3.Pos.X);
- miny = DMIN(DMIN(DMIN(miny, P1.Pos.Y), P2.Pos.Y), P3.Pos.Y);
- maxx = DMAX(DMAX(DMAX(maxx, P1.Pos.X), P2.Pos.X), P3.Pos.X);
- maxy = DMAX(DMAX(DMAX(maxy, P1.Pos.Y), P2.Pos.Y), P3.Pos.Y);
- }
-
- V[0].Pos.X = minx;
- V[0].Pos.Y = miny;
- V[1].Pos.X = maxx;
- V[1].Pos.Y = miny;
- V[2].Pos.X = maxx;
- V[2].Pos.Y = maxy;
- V[3].Pos.X = minx;
- V[3].Pos.Y = maxy;
-
- for(int i = 0; i < 4; i++) {
- V[i].Pos.Z = 0.0;
- putZ(&V[i], s);
- }
-
- n[0] = s->plan[2][0];
- n[1] = s->plan[2][1];
- n[2] = s->plan[2][2];
- norme(n);
-
- // compute which parts of the "middle cross" lie inside the surface
- const int numPoints = 200;
- for(int cross = 0; cross < 2; cross++) {
- int k = 0;
- for(int i = 0; i < numPoints; i++) {
- double t = (double)i / (double)(numPoints-1);
- if(!cross){
- vv.Pos.X = t * 0.5 * (V[0].Pos.X + V[1].Pos.X) +
- (1. - t) * 0.5 * (V[2].Pos.X + V[3].Pos.X);
- vv.Pos.Y = t * 0.5 * (V[0].Pos.Y + V[1].Pos.Y) +
- (1. - t) * 0.5 * (V[2].Pos.Y + V[3].Pos.Y);
- vv.Pos.Z = t * 0.5 * (V[0].Pos.Z + V[1].Pos.Z) +
- (1. - t) * 0.5 * (V[2].Pos.Z + V[3].Pos.Z);
+ // compute (rough) bounding box of surface
+ double minx = 0., miny = 0., maxx = 0., maxy = 0.;
+ for(int i = 0; i < List_Nbr(s->Generatrices); i++) {
+ List_Read(s->Generatrices, i, &c);
+ Vertex P1 = InterpolateCurve(c, 0.0, 0);
+ Vertex P2 = InterpolateCurve(c, 0.5, 0);
+ Vertex P3 = InterpolateCurve(c, 1.0, 0);
+ Projette(&P1, s->plan);
+ Projette(&P2, s->plan);
+ Projette(&P3, s->plan);
+ if(!i) {
+ minx = maxx = P1.Pos.X;
+ miny = maxy = P1.Pos.Y;
+ }
+ minx = DMIN(DMIN(DMIN(minx, P1.Pos.X), P2.Pos.X), P3.Pos.X);
+ miny = DMIN(DMIN(DMIN(miny, P1.Pos.Y), P2.Pos.Y), P3.Pos.Y);
+ maxx = DMAX(DMAX(DMAX(maxx, P1.Pos.X), P2.Pos.X), P3.Pos.X);
+ maxy = DMAX(DMAX(DMAX(maxy, P1.Pos.Y), P2.Pos.Y), P3.Pos.Y);
}
- else{
- vv.Pos.X = t * .5 * (V[0].Pos.X + V[3].Pos.X) +
- (1. - t) * .5 * (V[2].Pos.X + V[1].Pos.X);
- vv.Pos.Y = t * .5 * (V[0].Pos.Y + V[3].Pos.Y) +
- (1. - t) * .5 * (V[2].Pos.Y + V[1].Pos.Y);
- vv.Pos.Z = t * .5 * (V[0].Pos.Z + V[3].Pos.Z) +
- (1. - t) * .5 * (V[2].Pos.Z + V[1].Pos.Z);
+
+ V[0].Pos.X = minx;
+ V[0].Pos.Y = miny;
+ V[1].Pos.X = maxx;
+ V[1].Pos.Y = miny;
+ V[2].Pos.X = maxx;
+ V[2].Pos.Y = maxy;
+ V[3].Pos.X = minx;
+ V[3].Pos.Y = maxy;
+
+ for(int i = 0; i < 4; i++) {
+ V[i].Pos.Z = 0.0;
+ putZ(&V[i], s);
}
- if(isPointOnPlanarSurface(s, vv.Pos.X, vv.Pos.Y, vv.Pos.Z, n)) {
- if(!k) {
- List_Add(s->Orientations, &vv);
- k = 1;
- }
+
+ n[0] = s->plan[2][0];
+ n[1] = s->plan[2][1];
+ n[2] = s->plan[2][2];
+ norme(n);
+
+ // compute which parts of the "middle cross" lie inside the surface
+ const int numPoints = 200;
+ for(int cross = 0; cross < 2; cross++) {
+ int k = 0;
+ for(int i = 0; i < numPoints; i++) {
+ double t = (double)i / (double)(numPoints-1);
+ if(!cross){
+ vv.Pos.X = t * 0.5 * (V[0].Pos.X + V[1].Pos.X) +
+ (1. - t) * 0.5 * (V[2].Pos.X + V[3].Pos.X);
+ vv.Pos.Y = t * 0.5 * (V[0].Pos.Y + V[1].Pos.Y) +
+ (1. - t) * 0.5 * (V[2].Pos.Y + V[3].Pos.Y);
+ vv.Pos.Z = t * 0.5 * (V[0].Pos.Z + V[1].Pos.Z) +
+ (1. - t) * 0.5 * (V[2].Pos.Z + V[3].Pos.Z);
+ }
+ else{
+ vv.Pos.X = t * .5 * (V[0].Pos.X + V[3].Pos.X) +
+ (1. - t) * .5 * (V[2].Pos.X + V[1].Pos.X);
+ vv.Pos.Y = t * .5 * (V[0].Pos.Y + V[3].Pos.Y) +
+ (1. - t) * .5 * (V[2].Pos.Y + V[1].Pos.Y);
+ vv.Pos.Z = t * .5 * (V[0].Pos.Z + V[3].Pos.Z) +
+ (1. - t) * .5 * (V[2].Pos.Z + V[1].Pos.Z);
+ }
+ if(isPointOnPlanarSurface(s, vv.Pos.X, vv.Pos.Y, vv.Pos.Z, n)) {
+ if(!k) {
+ List_Add(s->Orientations, &vv);
+ k = 1;
+ }
+ }
+ else {
+ if(k) {
+ List_Add(s->Orientations, &vv);
+ k = 0;
+ }
+ }
+ }
+ if(k)
+ List_Add(s->Orientations, &vv);
}
- else {
- if(k) {
+
+ Msg(STATUS2, "Plane Surface %d (%d cross points)",
+ s->Num, List_Nbr(s->Orientations));
+
+ if(!List_Nbr(s->Orientations)){ // add dummy
List_Add(s->Orientations, &vv);
- k = 0;
- }
+ //printf("surf %d: min=%g %g max=%g %g\n", s->Num, minx, miny, maxx, maxy);
}
- }
- if(k)
- List_Add(s->Orientations, &vv);
- }
- Msg(STATUS2, "Plane Surface %d (%d cross points)",
- s->Num, List_Nbr(s->Orientations));
-
- if(!List_Nbr(s->Orientations)){ // add dummy
- List_Add(s->Orientations, &vv);
- //printf("surf %d: min=%g %g max=%g %g\n", s->Num, minx, miny, maxx, maxy);
+ CTX.threads_lock = 0;
}
-
- CTX.threads_lock = 0;
- }
- if(List_Nbr(s->Orientations) > 1) {
+ if(List_Nbr(s->Orientations) > 1) {
+
+ if(CTX.geom.surfaces) {
+ glEnable(GL_LINE_STIPPLE);
+ glLineStipple(1, 0x1F1F);
+ gl2psEnable(GL2PS_LINE_STIPPLE);
+ glBegin(GL_LINES);
+ for(int i = 0; i < List_Nbr(s->Orientations); i++) {
+ List_Read(s->Orientations, i, &vv);
+ glVertex3d(vv.Pos.X, vv.Pos.Y, vv.Pos.Z);
+ }
+ glEnd();
+ glDisable(GL_LINE_STIPPLE);
+ gl2psDisable(GL2PS_LINE_STIPPLE);
+ }
- if(CTX.geom.surfaces) {
- glEnable(GL_LINE_STIPPLE);
- glLineStipple(1, 0x1F1F);
- gl2psEnable(GL2PS_LINE_STIPPLE);
- glBegin(GL_LINES);
- for(int i = 0; i < List_Nbr(s->Orientations); i++) {
- List_Read(s->Orientations, i, &vv);
- glVertex3d(vv.Pos.X, vv.Pos.Y, vv.Pos.Z);
- }
- glEnd();
- glDisable(GL_LINE_STIPPLE);
- gl2psDisable(GL2PS_LINE_STIPPLE);
- }
+ if(CTX.geom.surfaces_num) {
+ List_Read(s->Orientations, 0, &vv1);
+ List_Read(s->Orientations, 1, &vv2);
+ char Num[100];
+ sprintf(Num, "%d", s->Num);
+ double offset = 0.3 * CTX.gl_fontsize * CTX.pixel_equiv_x;
+ glRasterPos3d((vv2.Pos.X + vv1.Pos.X) / 2. + offset / CTX.s[0],
+ (vv2.Pos.Y + vv1.Pos.Y) / 2. + offset / CTX.s[1],
+ (vv2.Pos.Z + vv1.Pos.Z) / 2. + offset / CTX.s[2]);
+ Draw_String(Num);
+ }
- if(CTX.geom.surfaces_num) {
- List_Read(s->Orientations, 0, &vv1);
- List_Read(s->Orientations, 1, &vv2);
- char Num[100];
- sprintf(Num, "%d", s->Num);
- double offset = 0.3 * CTX.gl_fontsize * CTX.pixel_equiv_x;
- glRasterPos3d((vv2.Pos.X + vv1.Pos.X) / 2. + offset / CTX.s[0],
- (vv2.Pos.Y + vv1.Pos.Y) / 2. + offset / CTX.s[1],
- (vv2.Pos.Z + vv1.Pos.Z) / 2. + offset / CTX.s[2]);
- Draw_String(Num);
- }
+ if(CTX.geom.normals) {
+ List_Read(s->Orientations, 0, &vv1);
+ List_Read(s->Orientations, 1, &vv2);
+ double x = (vv1.Pos.X + vv2.Pos.X) / 2.;
+ double y = (vv1.Pos.Y + vv2.Pos.Y) / 2.;
+ double z = (vv1.Pos.Z + vv2.Pos.Z) / 2.;
+ getPlaneSurfaceNormal(s, x, y, z, n);
+ n[0] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[0];
+ n[1] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[1];
+ n[2] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[2];
+ glColor4ubv((GLubyte *) & CTX.color.geom.normals);
+ Draw_Vector(CTX.vector_type, 0, CTX.arrow_rel_head_radius,
+ CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius,
+ x, y, z, n[0], n[1], n[2], CTX.geom.light);
+ }
- if(CTX.geom.normals) {
- List_Read(s->Orientations, 0, &vv1);
- List_Read(s->Orientations, 1, &vv2);
- double x = (vv1.Pos.X + vv2.Pos.X) / 2.;
- double y = (vv1.Pos.Y + vv2.Pos.Y) / 2.;
- double z = (vv1.Pos.Z + vv2.Pos.Z) / 2.;
- getPlaneSurfaceNormal(s, x, y, z, n);
- n[0] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[0];
- n[1] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[1];
- n[2] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[2];
- glColor4ubv((GLubyte *) & CTX.color.geom.normals);
- Draw_Vector(CTX.vector_type, 0, CTX.arrow_rel_head_radius,
- CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius,
- x, y, z, n[0], n[1], n[2], CTX.geom.light);
- }
-
- }
+ }
}
void Draw_NonPlane_Surface(Surface * s)
{
- if(CTX.geom.surfaces) {
- if(s->Typ == MSH_SURF_NURBS) {
- if(CTX.geom.light) glEnable(GL_LIGHTING);
- if(CTX.polygon_offset) glEnable(GL_POLYGON_OFFSET_FILL);
- GLUnurbsObj *nurb;
- nurb = gluNewNurbsRenderer();
+ if(CTX.geom.surfaces) {
+ if(s->Typ == MSH_SURF_NURBS) {
+ if(CTX.geom.light) glEnable(GL_LIGHTING);
+ if(CTX.polygon_offset) glEnable(GL_POLYGON_OFFSET_FILL);
+ GLUnurbsObj *nurb;
+ nurb = gluNewNurbsRenderer();
#if defined(GLU_VERSION_1_3)
- gluNurbsProperty(nurb, (GLenum) GLU_SAMPLING_TOLERANCE, 50.0);
- gluNurbsProperty(nurb, (GLenum) GLU_DISPLAY_MODE, GLU_FILL);
+ gluNurbsProperty(nurb, (GLenum) GLU_SAMPLING_TOLERANCE, 50.0);
+ gluNurbsProperty(nurb, (GLenum) GLU_DISPLAY_MODE, GLU_FILL);
#endif
- gluBeginSurface(nurb);
- gluNurbsSurface(nurb, s->Nu + s->OrderU + 1, s->ku,
- s->Nv + s->OrderV + 1, s->kv, 4, 4 * s->Nu, s->cp,
- s->OrderU + 1, s->OrderV + 1, GL_MAP2_VERTEX_4);
- gluEndSurface(nurb);
- gluDeleteNurbsRenderer(nurb);
- glDisable(GL_POLYGON_OFFSET_FILL);
- glDisable(GL_LIGHTING);
+ gluBeginSurface(nurb);
+ gluNurbsSurface(nurb, s->Nu + s->OrderU + 1, s->ku,
+ s->Nv + s->OrderV + 1, s->kv, 4, 4 * s->Nu, s->cp,
+ s->OrderU + 1, s->OrderV + 1, GL_MAP2_VERTEX_4);
+ gluEndSurface(nurb);
+ gluDeleteNurbsRenderer(nurb);
+ glDisable(GL_POLYGON_OFFSET_FILL);
+ glDisable(GL_LIGHTING);
+ }
+ else{
+ glEnable(GL_LINE_STIPPLE);
+ glLineStipple(1, 0x1F1F);
+ gl2psEnable(GL2PS_LINE_STIPPLE);
+ int N = 50;
+ glBegin(GL_LINE_STRIP);
+ for(int i = 0; i < N + 1; i++) {
+ double u = (double)i / (double)N;
+ Vertex v = InterpolateSurface(s, u, 0.5, 0, 0);
+ glVertex3d(v.Pos.X, v.Pos.Y, v.Pos.Z);
+ }
+ glEnd();
+ glBegin(GL_LINE_STRIP);
+ for(int i = 0; i < N + 1; i++) {
+ double u = (double)i / (double)N;
+ Vertex v = InterpolateSurface(s, 0.5, u, 0, 0);
+ glVertex3d(v.Pos.X, v.Pos.Y, v.Pos.Z);
+ }
+ glEnd();
+ glDisable(GL_LINE_STIPPLE);
+ gl2psDisable(GL2PS_LINE_STIPPLE);
+ }
+ }
+
+ if(CTX.geom.surfaces_num) {
+ Vertex v = InterpolateSurface(s, 0.5, 0.5, 0, 0);
+ char Num[100];
+ sprintf(Num, "%d", s->Num);
+ double offset = 0.3 * CTX.gl_fontsize * CTX.pixel_equiv_x;
+ glRasterPos3d(v.Pos.X + offset / CTX.s[0],
+ v.Pos.Y + offset / CTX.s[1],
+ v.Pos.Z + offset / CTX.s[2]);
+ Draw_String(Num);
+ }
+
+ if(CTX.geom.normals) {
+ Vertex v1 = InterpolateSurface(s, 0.5, 0.5, 0, 0);
+ Vertex v2 = InterpolateSurface(s, 0.6, 0.5, 0, 0);
+ Vertex v3 = InterpolateSurface(s, 0.5, 0.6, 0, 0);
+ double n[3];
+ normal3points(v1.Pos.X, v1.Pos.Y, v1.Pos.Z,
+ v2.Pos.X, v2.Pos.Y, v2.Pos.Z,
+ v3.Pos.X, v3.Pos.Y, v3.Pos.Z, n);
+ n[0] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[0];
+ n[1] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[1];
+ n[2] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[2];
+ glColor4ubv((GLubyte *) & CTX.color.geom.normals);
+ Draw_Vector(CTX.vector_type, 0, CTX.arrow_rel_head_radius,
+ CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius,
+ v1.Pos.X, v1.Pos.Y, v1.Pos.Z, n[0], n[1], n[2],
+ CTX.geom.light);
}
- else{
- glEnable(GL_LINE_STIPPLE);
- glLineStipple(1, 0x1F1F);
- gl2psEnable(GL2PS_LINE_STIPPLE);
- int N = 50;
- glBegin(GL_LINE_STRIP);
- for(int i = 0; i < N + 1; i++) {
- double u = (double)i / (double)N;
- Vertex v = InterpolateSurface(s, u, 0.5, 0, 0);
- glVertex3d(v.Pos.X, v.Pos.Y, v.Pos.Z);
- }
- glEnd();
- glBegin(GL_LINE_STRIP);
- for(int i = 0; i < N + 1; i++) {
- double u = (double)i / (double)N;
- Vertex v = InterpolateSurface(s, 0.5, u, 0, 0);
- glVertex3d(v.Pos.X, v.Pos.Y, v.Pos.Z);
- }
- glEnd();
- glDisable(GL_LINE_STIPPLE);
- gl2psDisable(GL2PS_LINE_STIPPLE);
- }
- }
-
- if(CTX.geom.surfaces_num) {
- Vertex v = InterpolateSurface(s, 0.5, 0.5, 0, 0);
- char Num[100];
- sprintf(Num, "%d", s->Num);
- double offset = 0.3 * CTX.gl_fontsize * CTX.pixel_equiv_x;
- glRasterPos3d(v.Pos.X + offset / CTX.s[0],
- v.Pos.Y + offset / CTX.s[1],
- v.Pos.Z + offset / CTX.s[2]);
- Draw_String(Num);
- }
-
- if(CTX.geom.normals) {
- Vertex v1 = InterpolateSurface(s, 0.5, 0.5, 0, 0);
- Vertex v2 = InterpolateSurface(s, 0.6, 0.5, 0, 0);
- Vertex v3 = InterpolateSurface(s, 0.5, 0.6, 0, 0);
- double n[3];
- normal3points(v1.Pos.X, v1.Pos.Y, v1.Pos.Z,
- v2.Pos.X, v2.Pos.Y, v2.Pos.Z,
- v3.Pos.X, v3.Pos.Y, v3.Pos.Z, n);
- n[0] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[0];
- n[1] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[1];
- n[2] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[2];
- glColor4ubv((GLubyte *) & CTX.color.geom.normals);
- Draw_Vector(CTX.vector_type, 0, CTX.arrow_rel_head_radius,
- CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius,
- v1.Pos.X, v1.Pos.Y, v1.Pos.Z, n[0], n[1], n[2],
- CTX.geom.light);
- }
}
void Draw_Surface(void *a, void *b)
{
- Surface *s = *(Surface **) a;
-
- if(!s || !s->Support || !(s->Visible & VIS_GEOM))
- return;
-
- if(CTX.render_mode == GMSH_SELECT) {
- glLoadName(2);
- glPushName(s->Num);
- }
-
- if(s->ipar[4]) {
- glLineWidth(CTX.geom.line_sel_width / 2.);
- gl2psLineWidth(CTX.geom.line_sel_width / 2. *
- CTX.print.eps_line_width_factor);
- glColor4ubv((GLubyte *) & CTX.color.geom.surface_sel);
- }
- else {
- glLineWidth(CTX.geom.line_width / 2.);
- gl2psLineWidth(CTX.geom.line_width / 2. * CTX.print.eps_line_width_factor);
- glColor4ubv((GLubyte *) & CTX.color.geom.surface);
- }
-
- if(s->bds){
- Draw_Polygonal_Surface(s);
- }
- else if(s->Typ == MSH_SURF_DISCRETE){
- // do nothing
- }
- else if(s->Typ == MSH_SURF_PLAN){
- Draw_Plane_Surface(s);
- }
- else{
- Draw_NonPlane_Surface(s);
- }
-
- if(CTX.render_mode == GMSH_SELECT) {
- glPopName();
- }
+ Surface *s = *(Surface **) a;
+
+ if(!s || !s->Support || !(s->Visible & VIS_GEOM))
+ return;
+
+ if(CTX.render_mode == GMSH_SELECT) {
+ glLoadName(2);
+ glPushName(s->Num);
+ }
+
+ if(s->ipar[4]) {
+ glLineWidth(CTX.geom.line_sel_width / 2.);
+ gl2psLineWidth(CTX.geom.line_sel_width / 2. *
+ CTX.print.eps_line_width_factor);
+ glColor4ubv((GLubyte *) & CTX.color.geom.surface_sel);
+ }
+ else {
+ glLineWidth(CTX.geom.line_width / 2.);
+ gl2psLineWidth(CTX.geom.line_width / 2. * CTX.print.eps_line_width_factor);
+ glColor4ubv((GLubyte *) & CTX.color.geom.surface);
+ }
+
+ if(s->bds){
+ Draw_Polygonal_Surface(s);
+ }
+ else if(s->Typ == MSH_SURF_DISCRETE){
+ // do nothing
+ }
+ else if(s->Typ == MSH_SURF_PLAN){
+ Draw_Plane_Surface(s);
+ }
+ else{
+ Draw_NonPlane_Surface(s);
+ }
+
+ if(CTX.render_mode == GMSH_SELECT) {
+ glPopName();
+ }
}
// Volumes
@@ -648,152 +642,152 @@ void DrawVolumes(Mesh * m)
void Draw_Geom(Mesh * m)
{
- if(m->status == -1)
- return;
-
- if(CTX.geom.points || CTX.geom.points_num)
- Tree_Action(m->Points, Draw_Geo_Point);
- if(CTX.geom.lines || CTX.geom.lines_num || CTX.geom.tangents)
- Tree_Action(m->Curves, Draw_Curve);
- if(CTX.geom.surfaces || CTX.geom.surfaces_num || CTX.geom.normals)
- Tree_Action(m->Surfaces, Draw_Surface);
- if(CTX.geom.volumes || CTX.geom.volumes_num)
- DrawVolumes(m);
+ if(m->status == -1)
+ return;
+
+ if(CTX.geom.points || CTX.geom.points_num)
+ Tree_Action(m->Points, Draw_Geo_Point);
+ if(CTX.geom.lines || CTX.geom.lines_num || CTX.geom.tangents)
+ Tree_Action(m->Curves, Draw_Curve);
+ if(CTX.geom.surfaces || CTX.geom.surfaces_num || CTX.geom.normals)
+ Tree_Action(m->Surfaces, Draw_Surface);
+ if(CTX.geom.volumes || CTX.geom.volumes_num)
+ DrawVolumes(m);
}
// Highlight routines
void HighlightEntity(Vertex * v, Curve * c, Surface * s, int permanent)
{
- Curve *cc;
- char Message[256], temp[256];
+ Curve *cc;
+ char Message[256], temp[256];
- if(permanent){
- // we want to draw incrementally (in-between to "Draw()" calls):
- // we need to make sure that the opengl context is set correctly
- SetOpenglContext();
- }
-
- if(v) {
- if(permanent){
- v->Frozen = 1;
- Draw_Geo_Point(&v,NULL);
- }
- else{
- Msg(STATUS1N, "Point %d {%.5g,%.5g,%.5g} (%.5g)", v->Num, v->Pos.X,
- v->Pos.Y, v->Pos.Z, v->lc);
- }
- }
- else if(c) {
if(permanent){
- c->ipar[3] = 1;
- Draw_Curve(&c,NULL);
- }
- else{
- if(c->beg && c->end)
- Msg(STATUS1N, "Curve %d {%d->%d}", c->Num, c->beg->Num, c->end->Num);
- else
- Msg(STATUS1N, "Curve %d", c->Num);
+ // we want to draw incrementally (in-between to "Draw()" calls):
+ // we need to make sure that the opengl context is set correctly
+ SetOpenglContext();
}
- }
- else if(s) {
- if(permanent){
- s->ipar[4] = 1;
- Draw_Surface(&s,NULL);
+
+ if(v) {
+ if(permanent){
+ v->Frozen = 1;
+ Draw_Geo_Point(&v,NULL);
+ }
+ else{
+ Msg(STATUS1N, "Point %d {%.5g,%.5g,%.5g} (%.5g)", v->Num, v->Pos.X,
+ v->Pos.Y, v->Pos.Z, v->lc);
+ }
}
- else{
- int nbg = List_Nbr(s->Generatrices);
- if(nbg){
- sprintf(Message, "Surface %d {", s->Num);
- if(nbg < 10) {
- for(int i = 0; i < nbg; i++) {
- List_Read(s->Generatrices, i, &cc);
- if(!i)
- sprintf(temp, "%d", cc->Num);
+ else if(c) {
+ if(permanent){
+ c->ipar[3] = 1;
+ Draw_Curve(&c,NULL);
+ }
+ else{
+ if(c->beg && c->end)
+ Msg(STATUS1N, "Curve %d {%d->%d}", c->Num, c->beg->Num, c->end->Num);
else
- sprintf(temp, ",%d", cc->Num);
- strcat(Message, temp);
- }
+ Msg(STATUS1N, "Curve %d", c->Num);
}
- else {
- strcat(Message, "...");
+ }
+ else if(s) {
+ if(permanent){
+ s->ipar[4] = 1;
+ Draw_Surface(&s,NULL);
}
- strcat(Message, "}");
- Msg(STATUS1N, Message);
- }
- else{
- Msg(STATUS1N, "Surface %d", s->Num);
- }
- }
- }
- else{
- if(!permanent)
- Msg(STATUS1N, " ");
- }
+ else{
+ int nbg = List_Nbr(s->Generatrices);
+ if(nbg){
+ sprintf(Message, "Surface %d {", s->Num);
+ if(nbg < 10) {
+ for(int i = 0; i < nbg; i++) {
+ List_Read(s->Generatrices, i, &cc);
+ if(!i)
+ sprintf(temp, "%d", cc->Num);
+ else
+ sprintf(temp, ",%d", cc->Num);
+ strcat(Message, temp);
+ }
+ }
+ else {
+ strcat(Message, "...");
+ }
+ strcat(Message, "}");
+ Msg(STATUS1N, Message);
+ }
+ else{
+ Msg(STATUS1N, "Surface %d", s->Num);
+ }
+ }
+ }
+ else{
+ if(!permanent)
+ Msg(STATUS1N, " ");
+ }
}
void HighlightEntityNum(int v, int c, int s, int permanent)
{
- if(v) {
- Vertex *pv = FindPoint(v, THEM);
- if(pv)
- HighlightEntity(pv, NULL, NULL, permanent);
- }
- if(c) {
- Curve *pc = FindCurve(c, THEM);
- if(pc)
- HighlightEntity(NULL, pc, NULL, permanent);
- }
- if(s) {
- Surface *ps = FindSurface(s, THEM);
- if(ps)
- HighlightEntity(NULL, NULL, ps, permanent);
- }
+ if(v) {
+ Vertex *pv = FindPoint(v, THEM);
+ if(pv)
+ HighlightEntity(pv, NULL, NULL, permanent);
+ }
+ if(c) {
+ Curve *pc = FindCurve(c, THEM);
+ if(pc)
+ HighlightEntity(NULL, pc, NULL, permanent);
+ }
+ if(s) {
+ Surface *ps = FindSurface(s, THEM);
+ if(ps)
+ HighlightEntity(NULL, NULL, ps, permanent);
+ }
}
void ZeroHighlightPoint(void *a, void *b)
{
- Vertex *v;
- v = *(Vertex **) a;
- v->Frozen = 0;
+ Vertex *v;
+ v = *(Vertex **) a;
+ v->Frozen = 0;
}
void ZeroHighlightCurve(void *a, void *b)
{
- Curve *c;
- c = *(Curve **) a;
- c->ipar[3] = 0;
+ Curve *c;
+ c = *(Curve **) a;
+ c->ipar[3] = 0;
}
void ZeroHighlightSurface(void *a, void *b)
{
- Surface *s;
- s = *(Surface **) a;
- s->ipar[4] = 0;
+ Surface *s;
+ s = *(Surface **) a;
+ s->ipar[4] = 0;
}
void ZeroHighlight(Mesh * m)
{
- Tree_Action(m->Points, ZeroHighlightPoint);
- Tree_Action(m->Curves, ZeroHighlightCurve);
- Tree_Action(m->Surfaces, ZeroHighlightSurface);
+ Tree_Action(m->Points, ZeroHighlightPoint);
+ Tree_Action(m->Curves, ZeroHighlightCurve);
+ Tree_Action(m->Surfaces, ZeroHighlightSurface);
}
void ZeroHighlightEntityNum(int v, int c, int s)
{
- if(v) {
- Vertex *pv = FindPoint(v, THEM);
- if(pv)
- ZeroHighlightPoint(&pv, NULL);
- }
- if(c) {
- Curve *pc = FindCurve(c, THEM);
- if(pc)
- ZeroHighlightCurve(&pc, NULL);
- }
- if(s) {
- Surface *ps = FindSurface(s, THEM);
- if(ps)
- ZeroHighlightSurface(&ps, NULL);
- }
+ if(v) {
+ Vertex *pv = FindPoint(v, THEM);
+ if(pv)
+ ZeroHighlightPoint(&pv, NULL);
+ }
+ if(c) {
+ Curve *pc = FindCurve(c, THEM);
+ if(pc)
+ ZeroHighlightCurve(&pc, NULL);
+ }
+ if(s) {
+ Surface *ps = FindSurface(s, THEM);
+ if(ps)
+ ZeroHighlightSurface(&ps, NULL);
+ }
}
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 4ba79899f8..7df60c5e87 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -3,6 +3,8 @@
#include <math.h>
#include <stdio.h>
#include "Numeric.h"
+#include "GmshMatrix.h"
+#include "EigSolve.h"
/*
(X-Xc)^2 = R^2
@@ -11,7 +13,28 @@
*/
+BDS_Vector::BDS_Vector (const BDS_Point &p2,const BDS_Point &p1)
+ : x(p2.X-p1.X),y(p2.Y-p1.Y),z(p2.Z-p1.Z)
+{
+}
+
+BDS_Vector BDS_Point::N() const
+{
+ std::list<BDS_Triangle *> t;
+ getTriangles (t);
+ std::list<BDS_Triangle *>::iterator tit = t.begin();
+ std::list<BDS_Triangle *>::iterator tite = t.end();
+ BDS_Vector SumN;
+ while (tit!=tite)
+ {
+ SumN += (*tit)->N();
+ ++tit;
+ }
+ double NORM = sqrt(SumN*SumN);
+ SumN /= NORM;
+ return SumN;
+}
BDS_Sphere :: BDS_Sphere( BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, BDS_Point *p4)
{
@@ -183,11 +206,11 @@ BDS_Plane::BDS_Plane (BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
}
-void BDS_Point::getTriangles (std::list<BDS_Triangle *> &t)
+void BDS_Point::getTriangles (std::list<BDS_Triangle *> &t) const
{
t.clear ();
- std::list<BDS_Edge*>::iterator it = edges.begin();
- std::list<BDS_Edge*>::iterator ite = edges.end();
+ std::list<BDS_Edge*>::const_iterator it = edges.begin();
+ std::list<BDS_Edge*>::const_iterator ite = edges.end();
while(it!=ite)
{
int NF = (*it)->numfaces();
@@ -342,6 +365,7 @@ void recur_tag ( BDS_Triangle *t , BDS_GeomEntity *g )
if (!t->g)
{
t->g = g;
+ g->t.push_back(t);
if ( ! t->e1->g && t->e1->numfaces() == 2)
{
recur_tag (t->e1->otherFace (t) , g);
@@ -377,18 +401,15 @@ void BDS_Mesh :: reverseEngineerCAD ( )
{
std::set<BDS_Point*,PointLessThan> pts;
- std::list<BDS_Triangle*>::iterator tit = triangles.begin();
- std::list<BDS_Triangle*>::iterator tite = triangles.end();
+ std::list<BDS_Triangle*>::iterator tit = (*it)->t.begin();
+ std::list<BDS_Triangle*>::iterator tite = (*it)->t.end();
while (tit!=tite)
{
- if ((*tit)->g == (*it))
- {
- BDS_Point *nod[3];
- (*tit)->getNodes (nod);
- pts.insert(nod[0]);
- pts.insert(nod[1]);
- pts.insert(nod[2]);
- }
+ BDS_Point *nod[3];
+ (*tit)->getNodes (nod);
+ pts.insert(nod[0]);
+ pts.insert(nod[1]);
+ pts.insert(nod[2]);
tit++;
}
if(pts.size() > 3)
@@ -412,7 +433,14 @@ void BDS_Mesh :: reverseEngineerCAD ( )
if (distmax < 1.e-6 * LC )
{
(*it)->surf = plane;
- printf("plane surface found %g\n",distmax);
+// printf("plane surface found %g\n",distmax);
+ pit = pts.begin();
+ while (pit!=pite)
+ {
+ if ((*pit)->g->classif_degree == 2)
+ (*pit)->radius_of_curvature = 1.e22;
+ pit++;
+ }
}
else
{
@@ -456,12 +484,100 @@ void BDS_Mesh :: reverseEngineerCAD ( )
}
}
+void BDS_Point :: compute_curvature ( )
+{
+ std::list<BDS_Triangle *> t;
+ getTriangles (t);
+
+// printf("curvature using %d triangles\n",t.size());
+
+ if (t.size() > 4)
+ {
+ Double_Matrix M ( t.size() , 4 );
+ Double_Vector Nx ( t.size() );
+ Double_Vector Ny ( t.size() );
+ Double_Vector Nz ( t.size() );
+ Double_Vector Cx ( 4 );
+ Double_Vector Cy ( 4 );
+ Double_Vector Cz ( 4 );
+
+ std::list<BDS_Triangle *>::iterator tit = t.begin();
+ std::list<BDS_Triangle *>::iterator tite = t.end();
+ int k = 0;
+ while (tit!=tite)
+ {
+ const BDS_Vector cog = (*tit)->cog();
+ M(k,0) = X-cog.x;
+ M(k,1) = Y-cog.y;
+ M(k,2) = Z-cog.z;
+ M(k,3) = 1.0;
+ BDS_Vector N = (*tit)->N();
+ Nx(k) = N.x;
+ Ny(k) = N.y;
+ Nz(k) = N.z;
+ k++;
+ ++tit;
+ }
+ M.least_squares (Nx,Cx);
+ M.least_squares (Ny,Cy);
+ M.least_squares (Nz,Cz);
+
+ const double A[9] = { Cx(0),0.5*(Cx(1)+Cy(0)), 0.5*(Cx(2)+Cz(0)) ,
+ 0.5*(Cx(1)+Cy(0)), Cy(1), 0.5*(Cy(2)+Cz(1)) ,
+ 0.5*(Cx(2)+Cz(0)), 0.5*(Cy(2)+Cz(1)) ,Cz(2)};
+
+ double v[9];
+ double wr[3],wi[3];
+
+ EigSolve3x3 (A ,wr,wi,v);
+
+ if(wr[0]!=0.0)radius_of_curvature = fabs(1/wr[0]);
+// printf(" curvature = %g %g %g R = %g\n",wr[0],wr[1],wr[2],radius_of_curvature);
+ }
+}
+
+
void BDS_Mesh :: classify ( double angle )
{
- printf(" classifying \n");
-
+ printf(" classifying \n");
+
static BDS_GeomEntity EDGE_CLAS (0,1);
+
+ // clear everything
+ // pre classified meshes will fail to work, more 2 be done
+ {
+ {
+ std::set<BDS_Point*,PointLessThan>::iterator it = points.begin();
+ std::set<BDS_Point*,PointLessThan>::iterator ite = points.end();
+ while (it != ite){ (*it)->g = 0; ++it;}
+ }
+
+ {
+ std::set<BDS_Edge*, EdgeLessThan>::iterator it = edges.begin();
+ std::set<BDS_Edge*, EdgeLessThan>::iterator ite = edges.end();
+ while (it != ite){ (*it)->g = 0; ++it;}
+ }
+ {
+ std::list<BDS_Triangle*>::iterator it = triangles.begin();
+ std::list<BDS_Triangle*>::iterator ite = triangles.end();
+ while (it!=ite){(*it)->g = 0;++it;}
+ }
+ geom.clear();
+ }
+
+
+ {
+ std::set<BDS_Point*,PointLessThan>::iterator it = points.begin();
+ std::set<BDS_Point*,PointLessThan>::iterator ite = points.end();
+ while (it != ite)
+ {
+ (*it)-> compute_curvature ( );
+ ++it;
+ }
+ }
+
+
{
std::set<BDS_Edge*, EdgeLessThan>::iterator it = edges.begin();
std::set<BDS_Edge*, EdgeLessThan>::iterator ite = edges.end();
@@ -472,19 +588,21 @@ void BDS_Mesh :: classify ( double angle )
{
e.g = &EDGE_CLAS;
}
-
+
else if (e.numfaces() == 2 &&
e.faces(0)->g != e.faces(1)->g )
- {
+ {
e.g = &EDGE_CLAS;
- }
-
+ }
+
else if (e.numfaces() == 2)
{
BDS_Vector N0 = e.faces(0)->N();
BDS_Vector N1 = e.faces(1)->N();
double ag = N0.angle (N1);
- if (fabs(ag) > angle)e.g = &EDGE_CLAS;
+ if (fabs(ag) > angle){
+ e.g = &EDGE_CLAS;
+ }
}
else
{
@@ -493,7 +611,8 @@ void BDS_Mesh :: classify ( double angle )
++it;
}
}
-
+
+/*
{
std::list<BDS_Triangle*>::iterator it = triangles.begin();
std::list<BDS_Triangle*>::iterator ite = triangles.end();
@@ -504,7 +623,10 @@ void BDS_Mesh :: classify ( double angle )
}
geom.clear();
}
+*/
+ printf(" classifying faces\n");
+
{
int tag = 1;
while (1)
@@ -527,7 +649,9 @@ void BDS_Mesh :: classify ( double angle )
recur_tag ( start , g );
tag++;
}
+ printf(" classifying edges (%d face tags found)\n",tag-1);
}
+
int edgetag = 1;
{
std::map< std::pair<int,int> , int > edgetags;
@@ -539,20 +663,31 @@ void BDS_Mesh :: classify ( double angle )
if (e.g)
{
+ int MIN_FAC = 100000;
+ int MAX_FAC = -100000;
std::map< std::pair<int,int> , int >::iterator found;
BDS_GeomEntity *g;
if ( e.numfaces() == 1)
{
found = edgetags.find (std::make_pair ( e.faces(0)->g->classif_tag,-1));
}
- else if (e.numfaces() == 2 && e.faces(0)->g->classif_tag != e.faces(1)->g->classif_tag)
+ else if (e.numfaces() == 2)// && e.faces(0)->g->classif_tag != e.faces(1)->g->classif_tag)
{
if (e.faces(0)->g->classif_tag > e.faces(1)->g->classif_tag)
found = edgetags.find (std::make_pair ( e.faces(1)->g->classif_tag,e.faces(0)->g->classif_tag));
else
found = edgetags.find (std::make_pair ( e.faces(0)->g->classif_tag,e.faces(1)->g->classif_tag));
}
-
+ else
+ {
+ for (int K=0;K<e.numfaces();K++)
+ {
+ if (MIN_FAC > e.faces(K)->g->classif_tag)MIN_FAC = e.faces(K)->g->classif_tag;
+ if (MAX_FAC < e.faces(K)->g->classif_tag)MAX_FAC = e.faces(K)->g->classif_tag;
+ found = edgetags.find (std::make_pair ( MIN_FAC,MAX_FAC ) );
+ }
+ }
+
if (e.g)
{
if (found == edgetags.end())
@@ -572,6 +707,11 @@ void BDS_Mesh :: classify ( double angle )
edgetags[std::make_pair ( e.faces(0)->g->classif_tag,e.faces(1)->g->classif_tag)]
= edgetag;
}
+ else
+ {
+ edgetags[std::make_pair ( MIN_FAC,MAX_FAC) ]
+ = edgetag;
+ }
edgetag++;
}
else
@@ -579,15 +719,20 @@ void BDS_Mesh :: classify ( double angle )
g = get_geom(found->second,1);
}
e.g = g;
+ g->e.push_back(&e);
}
}
else
{
e.g = e.faces(0)->g;
+ e.g->e.push_back(&e);
}
++it;
}
}
+
+ printf(" classifying points\n");
+
int vertextag = 1;
{
std::set<BDS_Point*,PointLessThan>::iterator it = points.begin();
@@ -625,37 +770,30 @@ void BDS_Mesh :: classify ( double angle )
std::set<BDS_Point*,PointLessThan>::iterator ite = points.end();
while (it != ite)
{
- if ((*it)->g->classif_degree > 1)
+ if ((*it)->g->classif_degree > 1)
{
- std::list<BDS_Triangle *> t;
- (*it)->getTriangles (t);
- std::list<BDS_Triangle *>::iterator tit = t.begin();
- std::list<BDS_Triangle *>::iterator tite = t.end();
- BDS_Vector SumN;
- while (tit!=tite)
- {
- SumN += (*tit)->N();
- ++tit;
- }
- double max_angle = 0;
- double NORM = sqrt(SumN*SumN);
- SumN /= NORM;
- tit = t.begin();
- while (tit!=tite)
+ std::list<BDS_Triangle *> t;
+ (*it)->getTriangles (t);
+ BDS_Vector SumN = (*it)->N();
+ std::list<BDS_Triangle *>::iterator tit = t.begin();
+ std::list<BDS_Triangle *>::iterator tite = t.end();
+ double max_angle = 0;
+ while (tit!=tite)
{
- double ag = SumN.angle ((*tit)->N());
- if (fabs(ag) > max_angle) max_angle = fabs(ag);
- ++tit;
+ double ag = SumN.angle ((*tit)->N());
+ if (fabs(ag) > max_angle) max_angle = fabs(ag);
+ ++tit;
}
- if (fabs(max_angle) > angle)
+ if (fabs(max_angle) > angle)
{
- add_geom (vertextag, 0);
+ add_geom (vertextag, 0);
(*it)->g = get_geom(vertextag++,0);
}
- }
- ++it;
+ }
+ ++it;
}
}
+ printf(" reverse engineering surfaces\n");
reverseEngineerCAD ( ) ;
printf(" end classifying %d edgetags %d vertextags\n",edgetag-1,vertextag-1);
@@ -1234,6 +1372,9 @@ bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord)
coord * p1->Y + (1-coord) * p2->Y,
coord * p1->Z + (1-coord) * p2->Z);
+ mid->radius_of_curvature =
+ coord * p1->radius_of_curvature + (1-coord) * p2->radius_of_curvature;
+
// we should project
e->oppositeof (op);
BDS_GeomEntity *g1=0,*g2=0,*ge=e->g;
@@ -1399,7 +1540,7 @@ bool BDS_Mesh ::swap_edge ( BDS_Edge *e)
triangles.push_back(t2);
return true;
}
-bool BDS_Mesh ::collapse_edge ( BDS_Edge *e, BDS_Point *p, const double eps)
+bool BDS_Mesh ::collapse_edge ( BDS_Edge *e, BDS_Point *p, const double eps, const double l)
{
if (e->numfaces() != 2)return false;
@@ -1439,25 +1580,29 @@ bool BDS_Mesh ::collapse_edge ( BDS_Edge *e, BDS_Point *p, const double eps)
// print_face(t);
if (t->e1 != e && t->e2 != e && t->e3 != e)
{
- double n1[3],n2[3];
+ BDS_Vector n1,n2;
BDS_Point *pts[3];
t->getNodes (pts);
p->X = o->X;
p->Y = o->Y;
p->Z = o->Z;
double s_after = surface_triangle (pts[0],pts[1],pts[2]);
- normal_triangle (pts[0],pts[1],pts[2], n1); // normal after collapse
+ n1 = t->N();
p->X = X;
p->Y = Y;
p->Z = Z;
- normal_triangle (pts[0],pts[1],pts[2], n2); // normal before collapse
+ n2 = t->N();
double s_before = surface_triangle (pts[0],pts[1],pts[2]);
// normals should not be opposed or change too dramatically
// this does not concern the triangles with the small edge that
// are collapsed too
- double ps = n1[0] * n2[0] + n1[1] * n2[1] + n1[2] * n2[2];
- if ( ps < eps ) return false;
- if (fabs (s_after - s_before) < 0.01 * fabs (s_after + s_before))return false;
+ double angle = n1.angle(n2);
+ if (fabs(angle) > M_PI/2 )return false;
+ if ( e->length() > 0.1 *l && fabs(angle) > M_PI/4 ) return false;
+ if (s_after < 1.e-2 * s_before)
+ {
+ return false;
+ }
gs[nt] = t->g;
pt[0][nt] = (pts[0] == p) ? o->iD : pts[0]->iD ;
pt[1][nt] = (pts[1] == p) ? o->iD : pts[1]->iD ;
@@ -1517,6 +1662,85 @@ bool BDS_Mesh ::collapse_edge ( BDS_Edge *e, BDS_Point *p, const double eps)
return true;
}
+/*
+
+ t(123)
+
+ p
+
+*/
+
+
+void project_point_on_a_list_of_triangles ( BDS_Point *p ,
+ const std::list<BDS_Triangle*> &t,
+ double &X, double &Y, double &Z)
+{
+ int p_classif = p->g->classif_tag;
+ BDS_Vector N = p->N();
+ std::list<BDS_Triangle *>::const_iterator it = t.begin();
+ std::list<BDS_Triangle *>::const_iterator ite = t.end() ;
+ double dist = 1.e22;
+ double XX = p->X;
+ double YY = p->Y;
+ double ZZ = p->Z;
+
+// printf("looking at %d triangles\n",t.size());
+
+ while (it != ite)
+ {
+ if ((*it)->g->classif_tag == p_classif && fabs(N.angle((*it)->N())) < (M_PI/6))
+ {
+ BDS_Point *pts[3];
+ (*it)->getNodes (pts);
+ double xp,yp,zp;
+ proj_point_plane ( p->X,p->Y,p->Z,pts[0],pts[1],pts[2],xp,yp,zp);
+ BDS_Point ppp (0,xp,yp,zp);
+ double a1 = surface_triangle ( &ppp, pts[0],pts[1] );
+ double a2 = surface_triangle ( &ppp, pts[1],pts[2] );
+ double a3 = surface_triangle ( &ppp, pts[2],pts[0] );
+ double a = surface_triangle ( pts[0], pts[1],pts[2] );
+ if ( fabs (a-a1-a2-a3) < 1.e-2 * a)
+ {
+ double d = sqrt ((xp-X)*(xp-X)+(yp-Y)*(yp-Y)+(zp-Z)*(zp-Z));
+ if (d < dist )
+ {
+// printf("one found %g %g %g %g\n",a,a1,a2,a3);
+ XX = xp; YY = yp; ZZ = zp;
+ dist = d;
+ }
+ }
+ }
+ ++it;
+ }
+ X = XX;
+ Y = YY;
+ Z = ZZ;
+}
+
+bool BDS_Mesh ::smooth_point_b ( BDS_Point *p )
+{
+ if (p->g && p->g->classif_degree <= 1)return false;
+
+ std::list<BDS_Edge*>::iterator eit = p->edges.begin();
+
+ double l_max = 0;
+ BDS_Edge *e_max;
+
+ while (eit!=p->edges.end())
+ {
+ ++eit;
+ }
+
+ if (p->g->surf)
+ {
+// p->g->surf->projection ( X,Y,Z,p->X,p->Y,p->Z);
+// printf("coucou\n");
+ }
+ else
+ {
+ }
+}
+
bool BDS_Mesh ::smooth_point ( BDS_Point *p , BDS_Mesh *geom_mesh )
{
@@ -1548,50 +1772,15 @@ bool BDS_Mesh ::smooth_point ( BDS_Point *p , BDS_Mesh *geom_mesh )
}
else
{
- std::list<BDS_Triangle *> t;
- std::list<BDS_Triangle *>::iterator it ;
- std::list<BDS_Triangle *>::iterator ite;
- if (geom_mesh)
- {
- it = geom_mesh->triangles.begin();
- ite = geom_mesh->triangles.end();
- }
- else
- {
- p->getTriangles (t);
- it = t.begin();
- ite = t.end();
- }
-
- double dist = 1.e22;
- double XX = X,YY=Y,ZZ=Z;
-
- int p_classif = p->g->classif_tag;
-
- while (it != ite)
- {
- if ((*it)->g->classif_tag == p_classif)
- {
- BDS_Point *pts[3];
- (*it)->getNodes (pts);
- double xp,yp,zp;
- proj_point_plane ( X,Y,Z,pts[0],pts[1],pts[2],xp,yp,zp);
- double d = sqrt ((xp-X)*(xp-X)+(yp-Y)*(yp-Y)+(zp-Z)*(zp-Z));
- if (d < dist)
- {
- XX = xp; YY = yp; ZZ = zp;
- dist = d;
- }
- }
- ++it;
- }
- X = XX;
- Y = YY;
- Z = ZZ;
-
+// return false;
p->X = X;
p->Y = Y;
p->Z = Z;
+ std::list<BDS_Triangle *>t;
+ p->getTriangles (t);
+ project_point_on_a_list_of_triangles ( p , t, p->X,p->Y,p->Z);
+ if (geom_mesh)
+ project_point_on_a_list_of_triangles ( p , geom_mesh->triangles, p->X,p->Y,p->Z);
}
return true;
}
@@ -1627,7 +1816,12 @@ int BDS_Mesh :: adapt_mesh ( double l, bool smooth, BDS_Mesh *geom_mesh)
double ll = sqrt((op[0]->X-op[1]->X)*(op[0]->X-op[1]->X)+
(op[0]->Y-op[1]->Y)*(op[0]->Y-op[1]->Y)+
(op[0]->Z-op[1]->Z)*(op[0]->Z-op[1]->Z));
- if (!(*it)->deleted && (*it)->length() > l/0.7 && ll > 0.25 * l){
+ double curv = 0.5 * (*it)->p1->radius_of_curvature +
+ 0.5 * (*it)->p2->radius_of_curvature;
+ double LLL = l;
+/// if (l > .50 * curv ) LLL = .50*curv;
+
+ if (!(*it)->deleted && (*it)->length() > LLL/0.7 && ll > 0.25 * LLL){
split_edge (*it, 0.5 );
nb_modif++;
}
@@ -1654,9 +1848,22 @@ int BDS_Mesh :: adapt_mesh ( double l, bool smooth, BDS_Mesh *geom_mesh)
std::set<BDS_Edge*, EdgeLessThan>::iterator ite = small_to_long.end();
while (it != ite)
{
- if (!(*it)->deleted && (*it)->length() < 0.7 * l ){
- if (collapse_edge (*it, (*it)->p1, 0.1 ))
- nb_modif++;
+ double curv = 0.5 * (*it)->p1->radius_of_curvature +
+ 0.5 * (*it)->p2->radius_of_curvature;
+ double LLL = l;
+// if (l > .50 * curv ) LLL = .50*curv;
+
+ if (!(*it)->deleted && (*it)->length() < 0.7 * LLL ){
+ if (nb_modif %2 )
+ {
+ if (collapse_edge (*it, (*it)->p1, 0.1,LLL ))
+ nb_modif++;
+ }
+ else
+ {
+ if (collapse_edge (*it, (*it)->p2, 0.1,LLL ))
+ nb_modif++;
+ }
}
++it;
}
@@ -1738,6 +1945,7 @@ int BDS_Mesh :: adapt_mesh ( double l, bool smooth, BDS_Mesh *geom_mesh)
BDS_Mesh::BDS_Mesh (const BDS_Mesh &other)
{
+ MAXPOINTNUMBER = other.MAXPOINTNUMBER;
for (std::set<BDS_GeomEntity*,GeomLessThan>::const_iterator it = other.geom.begin();
it != other.geom.end();
@@ -1753,6 +1961,7 @@ BDS_Mesh::BDS_Mesh (const BDS_Mesh &other)
{
BDS_Point *p = add_point((*it)->iD,(*it)->X,(*it)->Y,(*it)->Z);
p->g = ((*it)->g)? get_geom ((*it)->g->classif_tag,(*it)->g->classif_degree) : 0;
+ p->radius_of_curvature = (*it)->radius_of_curvature ;
}
for ( std::set<BDS_Edge*, EdgeLessThan>::const_iterator it = other.edges.begin();
it != other.edges.end();
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 0ea17a2b6d..97c1a988c9 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -30,6 +30,9 @@ public:
int classif_tag;
int classif_degree;
+ std::list<BDS_Triangle *> t;
+ std::list<BDS_Edge *> e;
+
BDS_Surface *surf;
inline bool operator < ( const BDS_GeomEntity & other ) const
@@ -48,39 +51,6 @@ public:
void print_face (BDS_Triangle *t);
-class BDS_Point
-{
-public:
- int iD;
- double X,Y,Z;
- BDS_GeomEntity *g;
- std::list<BDS_Edge*> edges;
- inline bool operator < ( const BDS_Point & other ) const
- {
- return iD < other.iD;
- }
- inline void del (BDS_Edge *e)
- {
- std::list<BDS_Edge*>::iterator it = edges.begin();
- std::list<BDS_Edge*>::iterator ite = edges.end();
- while(it!=ite)
- {
- if (*it == e)
- {
- edges.erase(it);
- break;
- }
- ++it;
- }
- }
- void getTriangles (std::list<BDS_Triangle *> &t);
-
- BDS_Point ( int id, double x=0, double y=0, double z=0 )
- : iD(id),X(x),Y(y),Z(z),g(0)
- {
- }
-};
-
class BDS_Vector
{
public:
@@ -141,11 +111,8 @@ public:
{
return (x*v.x+y*v.y+z*v.z);
}
- BDS_Vector (const BDS_Point &p2,const BDS_Point &p1)
- : x(p2.X-p1.X),y(p2.Y-p1.Y),z(p2.Z-p1.Z)
- {
-
- }
+ BDS_Vector (const BDS_Point &p2,const BDS_Point &p1);
+
BDS_Vector (double ix=0, double iy=0, double iz=0)
// : x(ix),(iy),z(iz)
{
@@ -156,6 +123,44 @@ public:
};
+class BDS_Point
+{
+public:
+ int iD;
+ double X,Y,Z,radius_of_curvature;
+ BDS_GeomEntity *g;
+ std::list<BDS_Edge*> edges;
+
+ BDS_Vector N() const;
+
+ inline bool operator < ( const BDS_Point & other ) const
+ {
+ return iD < other.iD;
+ }
+ inline void del (BDS_Edge *e)
+ {
+ std::list<BDS_Edge*>::iterator it = edges.begin();
+ std::list<BDS_Edge*>::iterator ite = edges.end();
+ while(it!=ite)
+ {
+ if (*it == e)
+ {
+ edges.erase(it);
+ break;
+ }
+ ++it;
+ }
+ }
+ void getTriangles (std::list<BDS_Triangle *> &t) const;
+
+ void compute_curvature ( );
+
+ BDS_Point ( int id, double x=0, double y=0, double z=0 )
+ : iD(id),X(x),Y(y),Z(z),radius_of_curvature(1.e22),g(0)
+ {
+ }
+};
+
class BDS_Edge
{
std::vector <BDS_Triangle *> _faces;
@@ -242,6 +247,15 @@ public:
BDS_Vector N() const ;
BDS_GeomEntity *g;
+ inline BDS_Vector cog() const
+ {
+ BDS_Point *n[3];
+ getNodes (n);
+ return BDS_Vector ((n[0]->X+n[1]->X+n[2]->X)/3.,
+ (n[0]->Y+n[1]->Y+n[2]->Y)/3.,
+ (n[0]->Z+n[1]->Z+n[2]->Z)/3.);
+ }
+
inline void getNodes (BDS_Point *n[3]) const
{
n[0] = e1->commonvertex (e3);
@@ -362,8 +376,9 @@ class BDS_Mesh
BDS_Edge *find_edge (BDS_Point *p1, BDS_Point *p2, BDS_Triangle *t)const;
BDS_GeomEntity *get_geom (int p1, int p2);
bool swap_edge ( BDS_Edge *);
- bool collapse_edge ( BDS_Edge *, BDS_Point*, const double eps);
- bool smooth_point ( BDS_Point* , BDS_Mesh *geom = 0);
+ bool collapse_edge ( BDS_Edge *, BDS_Point*, const double eps, const double l);
+ bool smooth_point ( BDS_Point* , BDS_Mesh *geom = 0);
+ bool smooth_point_b ( BDS_Point* );
bool split_edge ( BDS_Edge *, double coord);
void classify ( double angle);
void reverseEngineerCAD ( ) ;
@@ -378,3 +393,5 @@ class BDS_Mesh
void save_gmsh_format (const char *filename);
};
void normal_triangle (BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3]);
+void project_point_on_a_list_of_triangles ( BDS_Point *p , const std::list<BDS_Triangle*> &t,
+ double &X, double &Y, double &Z);
diff --git a/Mesh/DiscreteSurface.cpp b/Mesh/DiscreteSurface.cpp
index 80f2468eab..82d320df74 100644
--- a/Mesh/DiscreteSurface.cpp
+++ b/Mesh/DiscreteSurface.cpp
@@ -1,4 +1,4 @@
-// $Id: DiscreteSurface.cpp,v 1.15 2005-06-03 17:32:29 geuzaine Exp $
+// $Id: DiscreteSurface.cpp,v 1.16 2005-06-27 15:03:46 remacle Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -63,6 +63,13 @@ void Mesh_To_BDS(Surface *s, BDS_Mesh *m)
void BDS_To_Mesh(Mesh *m)
{
+ Tree_Action(m->Surfaces, Free_Surface);
+ Tree_Action(m->Curves, Free_Curve);
+ Tree_Delete(m->Surfaces);
+ Tree_Delete(m->Curves);
+ m->Curves = Tree_Create(sizeof(Curve *), compareCurve);
+ m->Surfaces = Tree_Create(sizeof(Surface *), compareSurface);
+
std::set<BDS_GeomEntity*,GeomLessThan>::iterator it = m->bds->geom.begin();
std::set<BDS_GeomEntity*,GeomLessThan>::iterator ite = m->bds->geom.end();
@@ -70,21 +77,30 @@ void BDS_To_Mesh(Mesh *m)
{
if ((*it)->classif_degree ==2 )
{
- Surface *_Surf = Create_Surface((*it)->classif_tag, MSH_SURF_DISCRETE);
+ Surface *_Surf = 0;
+ _Surf = FindSurface((*it)->classif_tag, m);
+ if (!_Surf)
+ _Surf = Create_Surface((*it)->classif_tag, MSH_SURF_DISCRETE);
_Surf->bds = m->bds;
End_Surface(_Surf);
Tree_Add(m->Surfaces, &_Surf);
}
else if ((*it)->classif_degree ==1 )
{
- Curve *_c = Create_Curve((*it)->classif_tag, MSH_SEGM_DISCRETE, 1, NULL, NULL, -1, -1, 0., 1.);
+ Curve *_c = 0;
+ _c = FindCurve((*it)->classif_tag, m);
+ if (!_c)
+ _c = Create_Curve((*it)->classif_tag, MSH_SEGM_DISCRETE, 1, NULL, NULL, -1, -1, 0., 1.);
_c->bds = m->bds;
End_Curve(_c);
Tree_Add(m->Curves, &_c);
}
++it;
}
- printf ("%d surfaces\n",Tree_Nbr(m->Surfaces));
+
+ CTX.mesh.changed = 1;
+
+ printf ("%d surfaces %d curves\n",Tree_Nbr(m->Surfaces), Tree_Nbr(m->Curves) );
}
diff --git a/Parser/OpenFile.cpp b/Parser/OpenFile.cpp
index 7d7856dfb8..9eced090e5 100644
--- a/Parser/OpenFile.cpp
+++ b/Parser/OpenFile.cpp
@@ -1,4 +1,4 @@
-// $Id: OpenFile.cpp,v 1.79 2005-05-16 00:50:11 geuzaine Exp $
+// $Id: OpenFile.cpp,v 1.80 2005-06-27 15:03:46 remacle Exp $
//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
@@ -293,6 +293,9 @@ int MergeProblem(char *name, int warn_if_missing)
}
THEM->bds->save_gmsh_format ( "1.msh" );
THEM->bds->classify ( M_PI / 8 );
+#if defined(HAVE_FLTK)
+ WID->create_surface_mesh_wizard();
+#endif
THEM->bds->save_gmsh_format ( "2.msh" );
BDS_To_Mesh (THEM);
SetBoundingBox();
--
GitLab