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

*** empty log message ***

parent 5494d7ff
No related branches found
No related tags found
No related merge requests found
// $Id: CommandLine.cpp,v 1.79 2006-09-03 07:44:09 geuzaine Exp $
// $Id: CommandLine.cpp,v 1.80 2006-09-07 16:03:32 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
......@@ -428,7 +428,9 @@ void Get_Options(int argc, char *argv[])
else if(!strncmp(argv[i], "netgen", 6))
CTX.mesh.algo3d = FRONTAL_NETGEN;
else if(!strncmp(argv[i], "tetgen", 6))
{
CTX.mesh.algo3d = DELAUNAY_TETGEN;
}
else {
fprintf(stderr, ERROR_STR "Unknown mesh algorithm\n");
exit(1);
......
// $Id: meshGFaceTransfinite.cpp,v 1.2 2006-09-06 10:25:24 remacle Exp $
// $Id: meshGFaceTransfinite.cpp,v 1.3 2006-09-07 16:03:32 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
......@@ -86,6 +86,15 @@ int MeshTransfiniteSurface( GFace *gf)
int L = N2-N1;
int H = N3-N2;
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;
}
std::vector<double> lengths_i;
std::vector<double> lengths_j;
......@@ -155,6 +164,8 @@ int MeshTransfiniteSurface( GFace *gf)
double VC3 = V[N3];
double VC4 = V[N4];
//create points using transfinite interpolation
for(int i = 1; i<L; i++)
{
double u = lengths_i[i]/L_i;
......@@ -176,8 +187,11 @@ int MeshTransfiniteSurface( GFace *gf)
}
}
if (gf->geomType() == GEntity::Plane) // some work 2 be done to smooth the stuff in parametric coords
// parametric coords are available (MFaceVertex) !
{
// elliptic smoother
for (int IT = 0;IT< 15;IT++)
for (int IT = 0;IT< 10;IT++)
{
for(int i = 1; i<L; i++)
{
......@@ -194,13 +208,20 @@ int MeshTransfiniteSurface( GFace *gf)
MVertex *v33 = tab[std::make_pair(i+1,j+1)];
double alpha = 0.25 * (DSQR(v23->x() - v21->x()) +
DSQR(v23->y() - v21->y()));
DSQR(v23->y() - v21->y()) +
DSQR(v23->z() - v21->z())
);
double gamma = 0.25 * (DSQR(v32->x() - v12->x()) +
DSQR(v32->y() - v12->y()));
DSQR(v32->y() - v12->y()) +
DSQR(v32->z() - v12->z())
);
double beta = 0.0625 * ((v32->x() - v12->x()) *
(v23->x() - v21->x()) +
(v32->y() - v12->y()) *
(v23->y() - v21->y()));
(v23->y() - v21->y()) +
(v32->z() - v12->z()) *
(v23->z() - v21->z())
);
v22->x() = 0.5 * (alpha * (v32->x() + v12->x()) +
gamma * (v23->x() + v21->x()) -
......@@ -220,7 +241,8 @@ int MeshTransfiniteSurface( GFace *gf)
}
}
}
}
// create elements
for(int i = 0; i < L ; i++)
{
for(int j = 0; j < H; j++)
......
......@@ -5,21 +5,10 @@
#include "Message.h"
#include <vector>
#if defined(HAVE_NETGEN)
namespace nglib {
#include "nglib.h"
#include "nglib_addon.h"
}
using namespace nglib;
Ng_Mesh * buildNetgenStructure (GRegion *gr, int importVolumeMesh, std::vector<MVertex*> & numberedV)
void getAllBoundingVertices (GRegion *gr, int importVolumeMesh, std::set<MVertex*> &allBoundingVertices )
{
NgAddOn_Init();
Ng_Mesh *_ngmesh = Ng_NewMesh();
std::list<GFace*> faces = gr->faces();
std::list<GFace*>::iterator it = faces.begin();
std::set<MVertex*> allBoundingVertices;
if (importVolumeMesh)
allBoundingVertices.insert ( gr->mesh_vertices.begin() , gr->mesh_vertices.end() );
......@@ -36,6 +25,121 @@ Ng_Mesh * buildNetgenStructure (GRegion *gr, int importVolumeMesh, std::vector<M
}
++it;
}
}
#if defined(HAVE_TETGEN)
#include "tetgen.h"
void buildTetgenStructure ( GRegion *gr, tetgenio &in, std::vector<MVertex*> & numberedV)
{
std::set<MVertex*> allBoundingVertices;
getAllBoundingVertices (gr, 0, allBoundingVertices );
in.mesh_dim = 3;
in.firstnumber = 1;
in.numberofpoints = allBoundingVertices.size();
in.pointlist = new REAL[in.numberofpoints * 3];
in.pointmarkerlist = NULL;
std::set<MVertex*>::iterator itv = allBoundingVertices.begin();
int I=1;
while (itv != allBoundingVertices.end())
{
in.pointlist[(I-1)*3 + 0] = (*itv)->x();
in.pointlist[(I-1)*3 + 1] = (*itv)->y();
in.pointlist[(I-1)*3 + 2] = (*itv)->z();
(*itv)->setNum(I++);
numberedV.push_back(*itv);
++itv;
}
int nbFace = 0;
std::list<GFace*> faces = gr->faces();
std::list<GFace*>::iterator it = faces.begin();
while (it != faces.end())
{
GFace *gf = (*it);
nbFace += gf->triangles.size();
++it;
}
Msg(INFO,"%d model faces -- NumFaces = %d",nbFace,faces.size());
in.numberoffacets = nbFace;
in.facetlist = new tetgenio::facet[in.numberoffacets];
in.facetmarkerlist = new int[in.numberoffacets];
it = faces.begin();
I= 0;
while (it != faces.end())
{
GFace *gf = (*it);
for (int i = 0; i< gf->triangles.size(); i++)
{
MTriangle *t = gf->triangles[i];
tetgenio::facet *f = &in.facetlist[I];
tetgenio::init(f);
f->numberofholes = 0;
f->numberofpolygons = 1;
tetgenio::polygon *p = f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
tetgenio::init(p);
p->numberofvertices = 3;
p->vertexlist = new int[p->numberofvertices];
p->vertexlist[0] = t->getVertex(0)->getNum();
p->vertexlist[1] = t->getVertex(1)->getNum();
p->vertexlist[2] = t->getVertex(2)->getNum();
// Msg(INFO,"Faces %d = %d %d %d",I,p->vertexlist[0] ,p->vertexlist[1] ,p->vertexlist[2] );
in.facetmarkerlist[I] = gf->tag();
++I;
}
++it;
}
}
void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out, std::vector<MVertex*> & numberedV)
{
int I = numberedV.size() + 1;
for (int i = 0; i < out.numberofpoints; i++)
{
MVertex *v = new MVertex (out.pointlist[i * 3 + 0],out.pointlist[i * 3 + 1],out.pointlist[i * 3 + 2]);
gr->mesh_vertices.push_back(v);
numberedV.push_back(v);
}
for (int i = 0; i < out.numberoftetrahedra; i++) {
MVertex *v1 = numberedV[out.tetrahedronlist[i * 4 + 0] -1];
MVertex *v2 = numberedV[out.tetrahedronlist[i * 4 + 1] -1];
MVertex *v3 = numberedV[out.tetrahedronlist[i * 4 + 2] -1];
MVertex *v4 = numberedV[out.tetrahedronlist[i * 4 + 3] -1];
MTetrahedron *t = new MTetrahedron(v1,v2,v3,v4);
gr->tetrahedra.push_back(t);
}
}
#endif
#if defined(HAVE_NETGEN)
namespace nglib {
#include "nglib.h"
#include "nglib_addon.h"
}
using namespace nglib;
Ng_Mesh * buildNetgenStructure (GRegion *gr, int importVolumeMesh, std::vector<MVertex*> & numberedV)
{
NgAddOn_Init();
Ng_Mesh *_ngmesh = Ng_NewMesh();
std::set<MVertex*> allBoundingVertices;
getAllBoundingVertices (gr, importVolumeMesh, allBoundingVertices );
std::set<MVertex*>::iterator itv = allBoundingVertices.begin();
int I=1;
while (itv != allBoundingVertices.end())
......@@ -50,7 +154,8 @@ Ng_Mesh * buildNetgenStructure (GRegion *gr, int importVolumeMesh, std::vector<M
++itv;
}
it = faces.begin();
std::list<GFace*> faces = gr->faces();
std::list<GFace*>::iterator it = faces.begin();
while (it != faces.end())
{
GFace *gf = (*it);
......@@ -261,6 +366,23 @@ void meshGRegion :: operator() (GRegion *gr)
// orient the triangles of with respect to this region
meshNormalsPointOutOfTheRegion (gr);
if(CTX.mesh.algo3d == DELAUNAY_TETGEN)
{
#if !defined(HAVE_TETGEN)
Msg(GERROR, "Tetgen is not compiled in this version of Gmsh");
#else
tetgenio in, out;
std::vector<MVertex*> numberedV;
char opts[128];
buildTetgenStructure ( gr, in, numberedV);
sprintf(opts, "pq1.4Ya%f%c", (float)CTX.mesh.quality,
(CTX.verbosity < 3)? 'Q': (CTX.verbosity > 6)? 'V': '\0');
Msg(STATUS2, "Meshing with volume constraint %f", (float)CTX.mesh.quality);
tetrahedralize(opts, &in, &out);
TransferTetgenMesh(gr, in, out, numberedV);
#endif
}
if(CTX.mesh.algo3d == FRONTAL_NETGEN)
{
#if !defined(HAVE_NETGEN)
......
......@@ -10,4 +10,4 @@ Line(4) = {1,4};
Line Loop(5) = {2,3,4,1};
Plane Surface(6) = {5};
Extrude Surface { 6, {0,0.0,1} };
Mesh.Algorithm3D = 4;
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment