From 37922a573f5bc3e31d23edccd439eba7d05cc0b5 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Mon, 27 Nov 2006 01:33:28 +0000
Subject: [PATCH] extrude edge+face (only "recombined" for now; we'll implement
subdivision into simplices "globally" later)
---
Common/Context.h | 2 +-
Common/DefaultOptions.h | 2 -
Common/Options.cpp | 9 +---
Common/Options.h | 1 -
Geo/gmshEdge.cpp | 3 +-
Geo/gmshRegion.cpp | 3 +-
Mesh/Makefile | 3 +-
Mesh/meshGEdge.cpp | 58 ++++++++++-----------
Mesh/meshGEdge.h | 2 +
Mesh/meshGFaceExtruded.cpp | 102 +++++++++++++++++++++++++++----------
10 files changed, 114 insertions(+), 71 deletions(-)
diff --git a/Common/Context.h b/Common/Context.h
index 915c2c33be..981c9c3324 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -185,7 +185,7 @@ public :
return val;
}
int oldxtrude, oldxtrude_recombine;
- int allow_degenerated_extrude, save_all, stl_binary, msh_binary, bdf_field_format;
+ int save_all, stl_binary, msh_binary, bdf_field_format;
char *triangle_options;
int smooth_normals, reverse_all_normals;
double angle_smooth_normals;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index b290b3ebcc..7dee45baaf 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -841,8 +841,6 @@ StringXNumber MeshOptions_Number[] = {
"2D mesh algorithm (1=isotropic, 2=anisotropic, 3=triangle)" },
{ F|O, "Algorithm3D" , opt_mesh_algo3d , DELAUNAY_ISO ,
"3D mesh algorithm (1=isotropic, 4=netgen, 5=tetgen)" },
- { F, "AllowDegeneratedExtrude" , opt_mesh_allow_degenerated_extrude , 0. ,
- "Allow the generation of degenerated hexahedra or prisms during extrusion" },
{ F|O, "AngleSmoothNormals" , opt_mesh_angle_smooth_normals , 30.0 ,
"Threshold angle below which normals are not smoothed" },
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 12aca03032..b4126c403b 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.317 2006-11-25 20:08:39 geuzaine Exp $
+// $Id: Options.cpp,v 1.318 2006-11-27 01:33:27 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -4976,13 +4976,6 @@ double opt_mesh_cut_planed(OPT_ARGS_NUM)
return CTX.mesh.cut_planed;
}
-double opt_mesh_allow_degenerated_extrude(OPT_ARGS_NUM)
-{
- if(action & GMSH_SET)
- CTX.mesh.allow_degenerated_extrude = (int)val;
- return CTX.mesh.allow_degenerated_extrude;
-}
-
double opt_mesh_save_all(OPT_ARGS_NUM)
{
if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index 78f5b802c0..dc6eb0d607 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -478,7 +478,6 @@ double opt_mesh_cut_planea(OPT_ARGS_NUM);
double opt_mesh_cut_planeb(OPT_ARGS_NUM);
double opt_mesh_cut_planec(OPT_ARGS_NUM);
double opt_mesh_cut_planed(OPT_ARGS_NUM);
-double opt_mesh_allow_degenerated_extrude(OPT_ARGS_NUM);
double opt_mesh_save_all(OPT_ARGS_NUM);
double opt_mesh_color_carousel(OPT_ARGS_NUM);
double opt_mesh_nb_nodes(OPT_ARGS_NUM);
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index caed9ae6ef..4fd610c710 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshEdge.cpp,v 1.19 2006-11-25 16:52:43 geuzaine Exp $
+// $Id: gmshEdge.cpp,v 1.20 2006-11-27 01:33:27 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -34,6 +34,7 @@ gmshEdge::gmshEdge(GModel *model, Curve *edge, GVertex *v1, GVertex *v2)
meshAttributes.nbPointsTransfinite = c->ipar[0];
meshAttributes.coeffTransfinite = c->dpar[0];
meshAttributes.typeTransfinite = c->ipar[1];
+ meshAttributes.extrude = c->Extrude;
}
gmshEdge::gmshEdge(GModel *model, int num)
diff --git a/Geo/gmshRegion.cpp b/Geo/gmshRegion.cpp
index 280721bf22..d29ab45453 100644
--- a/Geo/gmshRegion.cpp
+++ b/Geo/gmshRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshRegion.cpp,v 1.8 2006-11-25 16:52:43 geuzaine Exp $
+// $Id: gmshRegion.cpp,v 1.9 2006-11-27 01:33:27 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -39,6 +39,7 @@ gmshRegion::gmshRegion(GModel *m, ::Volume * volume)
l_faces.push_back(f);
l_dirs.push_back(ori);
}
+ meshAttributes.extrude = volume->Extrude;
}
gmshRegion::gmshRegion(GModel *m, int num)
diff --git a/Mesh/Makefile b/Mesh/Makefile
index b290d039e8..5cef9c9845 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.147 2006-11-26 04:36:46 geuzaine Exp $
+# $Id: Makefile,v 1.148 2006-11-27 01:33:28 geuzaine Exp $
#
# Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
#
@@ -32,6 +32,7 @@ CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = Generator.cpp \
meshGEdge.cpp \
+ meshGEdgeExtruded.cpp \
meshGFace.cpp \
meshGFaceTransfinite.cpp \
meshGFaceExtruded.cpp \
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index c8e70dc1b1..cc2ed927de 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdge.cpp,v 1.19 2006-11-25 18:03:49 geuzaine Exp $
+// $Id: meshGEdge.cpp,v 1.20 2006-11-27 01:33:28 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -38,20 +38,20 @@ double F_LC_ANALY (double xx, double yy, double zz)
// return 0.005 + 0.05*fabs (sin(5*xx) + sin(15*yy) + sin(15*zz));
// return 0.02;
// return 0.002 + 0.04*fabs (sin(6*xx) + sin(6*yy) + sin(6*zz));
- return 0.003 + 0.05*fabs (sin(8*xx) + sin(8*yy) + sin(8*zz));
- return 0.02 + 0.1*fabs (sin(3*xx) + sin(3*yy) + sin(3*zz));
- return 0.01 + 0.1*fabs(sin((xx*xx+(zz-0.7)*(zz-0.7)-.25))) ;
- return 0.05 + .1*fabs(xx*yy) ;
+ return 0.003 + 0.05*fabs(sin(8*xx) + sin(8*yy) + sin(8*zz));
+ return 0.02 + 0.1*fabs(sin(3*xx) + sin(3*yy) + sin(3*zz));
+ return 0.01 + 0.1*fabs(sin((xx*xx+(zz-0.7)*(zz-0.7)-.25)));
+ return 0.05 + 0.1*fabs(xx*yy);
}
-double F_Lc_bis(double t)
+double F_Lc(double t)
{
// const double nb_points_per_radius_of_curv = 2;
- GPoint point = _myGEdge -> point(t) ;
- const double fact = (t-t_begin)/(t_end-t_begin);
- double lc_here = lc_begin + fact * (lc_end-lc_begin);
- SVector3 der = _myGEdge -> firstDer(t) ;
- const double d = norm(der);
+ GPoint point = _myGEdge->point(t) ;
+ const double fact = (t - t_begin) / (t_end - t_begin);
+ double lc_here = lc_begin + fact * (lc_end - lc_begin);
+ SVector3 der = _myGEdge->firstDer(t) ;
+ const double d = norm(der);
if(CTX.mesh.bgmesh_type == ONFILE) {
const double Lc = BGMXYZ(point.x(), point.y(), point.z());
@@ -64,16 +64,16 @@ double F_Lc_bis(double t)
return d/lc_here;
}
-double F_Transfini_bis(double t)
+double F_Transfinite(double t)
{
double val, r;
- SVector3 der = _myGEdge -> firstDer(t) ;
+ SVector3 der = _myGEdge->firstDer(t) ;
double d = norm(der);
double coef = _myGEdge->meshAttributes.coeffTransfinite;
- int type = _myGEdge->meshAttributes.typeTransfinite;
- int nbpt = _myGEdge->meshAttributes.nbPointsTransfinite;
+ int type = _myGEdge->meshAttributes.typeTransfinite;
+ int nbpt = _myGEdge->meshAttributes.nbPointsTransfinite;
if(coef <= 0.0 || coef == 1.0) {
// coef < 0 should never happen
@@ -121,9 +121,9 @@ double F_Transfini_bis(double t)
return val;
}
-double F_One_bis(double t)
+double F_One(double t)
{
- SVector3 der = _myGEdge -> firstDer(t) ;
+ SVector3 der = _myGEdge->firstDer(t) ;
return norm(der);
}
@@ -132,7 +132,7 @@ typedef struct{
double t, lc, p;
}IntPoint;
-double trapeze(IntPoint * P1, IntPoint * P2)
+double trapezoidal(IntPoint * P1, IntPoint * P2)
{
return (0.5 * (P1->lc + P2->lc) * (P2->t - P1->t));
}
@@ -149,9 +149,9 @@ void RecursiveIntegration(IntPoint * from, IntPoint * to,
P.t = 0.5 * (from->t + to->t);
P.lc = f(P.t);
- val1 = trapeze(from, to);
- val2 = trapeze(from, &P);
- val3 = trapeze(&P, to);
+ val1 = trapezoidal(from, to);
+ val2 = trapezoidal(from, &P);
+ val3 = trapezoidal(&P, to);
err = fabs(val1 - val2 - val3);
// Msg(INFO,"Int %22.15 E %22.15 E %22.15 E\n", val1,val2,val3);
@@ -210,6 +210,8 @@ void meshGEdge :: operator() (GEdge *ge)
deMeshGEdge dem;
dem(ge);
+ if(MeshExtrudedCurve(ge)) return;
+
// Create a list of integration points
List_T *Points = List_Create(10, 10, sizeof(IntPoint));
@@ -218,31 +220,30 @@ void meshGEdge :: operator() (GEdge *ge)
// to pass an extra argument...
_myGEdge = ge;
-
// compute bounds
_myGEdgeBounds = ge->parBounds(0) ;
t_begin = _myGEdgeBounds.low();
- t_end = _myGEdgeBounds.high();
+ t_end = _myGEdgeBounds.high();
// first compute the length of the curve by integrating one
_myGEdgeLength = Integration(_myGEdgeBounds.low(), _myGEdgeBounds.high(),
- F_One_bis, Points, 1.e-4);
+ F_One, Points, 1.e-4);
List_Reset(Points);
- lc_begin = _myGEdge->getBeginVertex()->prescribedMeshSizeAtVertex();
- lc_end = _myGEdge->getEndVertex()->prescribedMeshSizeAtVertex();
+ lc_begin = _myGEdge->getBeginVertex()->prescribedMeshSizeAtVertex();
+ lc_end = _myGEdge->getEndVertex()->prescribedMeshSizeAtVertex();
// Integrate detJ/lc du
double a;
int N;
if(ge->meshAttributes.Method == TRANSFINI){
a = Integration(_myGEdgeBounds.low(), _myGEdgeBounds.high(),
- F_Transfini_bis, Points, 1.e-7);
+ F_Transfinite, Points, 1.e-7);
N = ge->meshAttributes.nbPointsTransfinite;
}
else{
a = Integration(_myGEdgeBounds.low(), _myGEdgeBounds.high(),
- F_Lc_bis, Points, 1.e-4);
+ F_Lc, Points, 1.e-4);
N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.));
}
const double b = a / (double)(N - 1);
@@ -282,4 +283,3 @@ void meshGEdge :: operator() (GEdge *ge)
ge->lines.push_back(new MLine(v0, v1));
}
}
-
diff --git a/Mesh/meshGEdge.h b/Mesh/meshGEdge.h
index 40d22be02b..03bce0c693 100644
--- a/Mesh/meshGEdge.h
+++ b/Mesh/meshGEdge.h
@@ -35,4 +35,6 @@ class deMeshGEdge
void operator () ( GEdge * );
};
+int MeshExtrudedCurve(GEdge *ge);
+
#endif
diff --git a/Mesh/meshGFaceExtruded.cpp b/Mesh/meshGFaceExtruded.cpp
index 7264c806bf..68e9554a84 100644
--- a/Mesh/meshGFaceExtruded.cpp
+++ b/Mesh/meshGFaceExtruded.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceExtruded.cpp,v 1.4 2006-11-26 21:56:26 geuzaine Exp $
+// $Id: meshGFaceExtruded.cpp,v 1.5 2006-11-27 01:33:28 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -24,7 +24,8 @@
#include "GModel.h"
#include "Message.h"
-int extrudeMesh(GEdge *from, GFace *to)
+void extrudeMesh(GEdge *from, GFace *to,
+ std::set<MVertex*, MVertexLessThanLexicographic> &pos)
{
ExtrudeParams *ep = to->meshAttributes.extrude;
@@ -32,35 +33,69 @@ int extrudeMesh(GEdge *from, GFace *to)
for(unsigned int i = 0; i < from->mesh_vertices.size(); i++){
MVertex *v = from->mesh_vertices[i];
for(int j = 0; j < ep->mesh.NbLayer; j++) {
- for(int k = 1; k < ep->mesh.NbElmLayer[j]; k++) {
+ for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
double x = v->x(), y = v->y(), z = v->z();
- ep->Extrude(j, k, x, y, z);
- to->mesh_vertices.push_back(new MVertex(x, y, z, to));
+ ep->Extrude(j, k + 1, x, y, z);
+ if(j != ep->mesh.NbLayer - 1 || k != ep->mesh.NbElmLayer[j] - 1){
+ MVertex *newv = new MVertex(x, y, z, to);
+ to->mesh_vertices.push_back(newv);
+ pos.insert(newv);
+ }
}
}
}
// create elements
+ std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp;
+ for(unsigned int i = 0; i < from->lines.size(); i++){
+ MVertex *v0 = from->lines[i]->getVertex(0);
+ MVertex *v1 = from->lines[i]->getVertex(1);
+ for(int j = 0; j < ep->mesh.NbLayer; j++) {
+ for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
+ std::vector<MVertex*> verts;
+ double x[4], y[4], z[4];
+ x[0] = v0->x(); y[0] = v0->y(); z[0] = v0->z();
+ ep->Extrude(j, k, x[0], y[0], z[0]);
+ x[1] = v1->x(); y[1] = v1->y(); z[1] = v1->z();
+ ep->Extrude(j, k, x[1], y[1], z[1]);
+ x[2] = v0->x(); y[2] = v0->y(); z[2] = v0->z();
+ ep->Extrude(j, k + 1, x[2], y[2], z[2]);
+ x[3] = v1->x(); y[3] = v1->y(); z[3] = v1->z();
+ ep->Extrude(j, k + 1, x[3], y[3], z[3]);
+ for(int p = 0; p < 4; p++){
+ MVertex tmp(x[p], y[p], z[p], 0, -1);
+ itp = pos.find(&tmp);
+ if(itp == pos.end()) {
+ Msg(GERROR, "Could not find extruded vertex in surface %d", to->tag());
+ return;
+ }
+ verts.push_back(*itp);
+ }
+ if(verts[0] == verts[1] || verts[1] == verts[3]){
+ to->triangles.push_back(new MTriangle(verts[0], verts[3], verts[2]));
+ }
+ else if(verts[0] == verts[2] || verts[2] == verts[3]){
+ to->triangles.push_back(new MTriangle(verts[0], verts[1], verts[3]));
+ }
+ else if(verts[0] == verts[3] || verts[1] == verts[2]){
+ Msg(GERROR, "Uncoherent quadrangle in extrusion");
+ return;
+ }
+ else{
+ to->quadrangles.push_back(new MQuadrangle(verts[0], verts[1],
+ verts[3], verts[2]));
+ }
+ }
+ }
+ }
}
-int copyMesh(GFace *from, GFace *to)
+void copyMesh(GFace *from, GFace *to,
+ std::set<MVertex*, MVertexLessThanLexicographic> &pos)
{
ExtrudeParams *ep = to->meshAttributes.extrude;
- // build a set with all the vertices on the boundary of to
- std::set<MVertex*, MVertexLessThanLexicographic> pos;
- std::list<GEdge*> edges = to->edges();
- std::list<GEdge*>::iterator it = edges.begin();
- while(it != edges.end()){
- pos.insert((*it)->mesh_vertices.begin(), (*it)->mesh_vertices.end());
- pos.insert((*it)->getBeginVertex()->mesh_vertices.begin(),
- (*it)->getBeginVertex()->mesh_vertices.end());
- pos.insert((*it)->getEndVertex()->mesh_vertices.begin(),
- (*it)->getEndVertex()->mesh_vertices.end());
- ++it;
- }
-
- // create new vertices
+ // create vertices
for(unsigned int i = 0; i < from->mesh_vertices.size(); i++){
MVertex *v = from->mesh_vertices[i];
double x = v->x(), y = v->y(), z = v->z();
@@ -71,19 +106,19 @@ int copyMesh(GFace *from, GFace *to)
pos.insert(newv);
}
- // create new elements
+ // create elements
std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp;
for(unsigned int i = 0; i < from->triangles.size(); i++){
std::vector<MVertex*> verts;
for(int j = 0; j < 3; j++){
MVertex *v = from->triangles[i]->getVertex(j);
- MVertex tmp(v->x(), v->y(), v->z(), 0, 0);
+ MVertex tmp(v->x(), v->y(), v->z(), 0, -1);
ep->Extrude(ep->mesh.NbLayer - 1, ep->mesh.NbElmLayer[ep->mesh.NbLayer - 1],
tmp.x(), tmp.y(), tmp.z());
itp = pos.find(&tmp);
if(itp == pos.end()) {
Msg(GERROR, "Could not find extruded vertex in surface %d", to->tag());
- return 0;
+ return;
}
verts.push_back(*itp);
}
@@ -93,13 +128,13 @@ int copyMesh(GFace *from, GFace *to)
std::vector<MVertex*> verts;
for(int j = 0; j < 4; j++){
MVertex *v = from->quadrangles[i]->getVertex(j);
- MVertex tmp(v->x(), v->y(), v->z(), 0, 0);
+ MVertex tmp(v->x(), v->y(), v->z(), 0, -1);
ep->Extrude(ep->mesh.NbLayer - 1, ep->mesh.NbElmLayer[ep->mesh.NbLayer - 1],
tmp.x(), tmp.y(), tmp.z());
itp = pos.find(&tmp);
if(itp == pos.end()) {
Msg(GERROR, "Could not find extruded vertex in surface %d", to->tag());
- return 0;
+ return;
}
verts.push_back(*itp);
}
@@ -114,17 +149,30 @@ int MeshExtrudedSurface(GFace *gf)
if(!ep || !ep->mesh.ExtrudeMesh)
return 0;
+ // build a set with all the vertices on the boundary of gf
+ std::set<MVertex*, MVertexLessThanLexicographic> pos;
+ std::list<GEdge*> edges = gf->edges();
+ std::list<GEdge*>::iterator it = edges.begin();
+ while(it != edges.end()){
+ pos.insert((*it)->mesh_vertices.begin(), (*it)->mesh_vertices.end());
+ pos.insert((*it)->getBeginVertex()->mesh_vertices.begin(),
+ (*it)->getBeginVertex()->mesh_vertices.end());
+ pos.insert((*it)->getEndVertex()->mesh_vertices.begin(),
+ (*it)->getEndVertex()->mesh_vertices.end());
+ ++it;
+ }
+
if(ep->geo.Mode == EXTRUDED_ENTITY) {
// surface is extruded from a curve
GEdge *from = gf->model()->edgeByTag(std::abs(ep->geo.Source));
if(!from) return 0;
- extrudeMesh(from, gf);
+ extrudeMesh(from, gf, pos);
}
else {
// surface is a copy of another surface ("chapeau")
GFace *from = gf->model()->faceByTag(std::abs(ep->geo.Source));
if(!from) return 0;
- copyMesh(from, gf);
+ copyMesh(from, gf, pos);
}
return 1;
--
GitLab