Select Git revision
DefaultOptions.h
Forked from
gmsh / gmsh
Source project has a limited visibility.
2D_Mesh.cpp 23.78 KiB
// $Id: 2D_Mesh.cpp,v 1.80 2006-01-06 00:34:25 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA.
//
// Please report all bugs and problems to <gmsh@geuz.org>.
#include "Gmsh.h"
#include "Numeric.h"
#include "Geo.h"
#include "CAD.h"
#include "Mesh.h"
#include "Utils.h"
#include "Create.h"
#include "2D_Mesh.h"
#include "Context.h"
#include "Interpolation.h"
extern Mesh *THEM;
extern Context_T CTX;
PointRecord *gPointArray;
DocRecord *BGMESH, *FGMESH;
double LC2D;
static int is_3D = 0;
static Surface *THESURFACE, *THESUPPORT;
void ProjetteSurface(void *a, void *b)
{
Vertex *v;
v = *(Vertex **) a;
if(!v->ListCurves)
ProjectPointOnSurface(THESUPPORT, *v);
}
void Projette_Plan_Moyen(void *a, void *b)
{
Vertex *v = *(Vertex **) a;
Projette(v, THESURFACE->plan);
}
void Projette_Inverse(void *a, void *b)
{
Vertex *v = *(Vertex **) a;
Projette(v, THESURFACE->invplan);
}
void Plan_Moyen(void *data, void *dum)
{
int i, j;
Vertex *v;
Curve *pC;
static List_T *points;
static int deb = 1;
Surface *s = *(Surface **) data;
if(deb) {
points = List_Create(10, 10, sizeof(Vertex *));
deb = 0;
}
else
List_Reset(points);
switch (s->Typ) {
case MSH_SURF_PLAN:
case MSH_SURF_TRIC:
case MSH_SURF_REGL:
case MSH_SURF_NURBS:
case MSH_SURF_TRIMMED:
for(i = 0; i < List_Nbr(s->Generatrices); i++) {
List_Read(s->Generatrices, i, &pC);
for(j = 0; j < List_Nbr(pC->Vertices); j++) {
List_Read(pC->Vertices, j, &v);
List_Add(points, &v);
Tree_Insert(s->Vertices, List_Pointer(pC->Vertices, j));
}
}
break;
}
MeanPlane(points, s);
}
int Calcule_Contours(Surface * s)
{
int i, j, ori, ORI = 0;
Vertex *v, *first, *last;
List_T *l, *linv;
Curve *c;
double n[] = { 0., 0., 1. };
l = List_Create(10, 10, sizeof(Vertex *));
for(i = 0; i < List_Nbr(s->Generatrices); i++) {
List_Read(s->Generatrices, i, &c);
if(!List_Nbr(l)) {
List_Read(c->Vertices, 0, &first);
}
for(j = 1; j < List_Nbr(c->Vertices); j++) {
List_Read(c->Vertices, j, &v);
List_Add(l, &v);
if(j == List_Nbr(c->Vertices) - 1) {
last = v;
}
}
if(!compareVertex(&last, &first)) {
ori = Oriente(l, n);
/* Le premier contour est oriente aire a droite
Les autes sont des trous orientes aire a gauche */
if(!List_Nbr(s->Contours))
ORI = ori;
if((!List_Nbr(s->Contours) && !ori) || (List_Nbr(s->Contours) && ori)) {
linv = List_Create(10, 10, sizeof(Vertex *));
List_Invert(l, linv);
List_Delete(l);
l = linv;
}
List_Add(s->Contours, &l);
l = List_Create(10, 10, sizeof(Vertex *));
} }
return ORI;
}
void Calcule_Z(void *data, void *dum)
{
Vertex **pv, *v;
double Z, U, V;
pv = (Vertex **) data;
v = *pv;
if(v->Frozen || v->Num < 0)
return;
XYtoUV(THESUPPORT, &v->Pos.X, &v->Pos.Y, &U, &V, &Z, 1.0);
v->Pos.Z = Z;
}
void Calcule_Z_Plan(void *data, void *dum)
{
Vertex **pv, *v;
Vertex V;
V.Pos.X = THESURFACE->a;
V.Pos.Y = THESURFACE->b;
V.Pos.Z = THESURFACE->c;
Projette(&V, THESURFACE->plan);
pv = (Vertex **) data;
v = *pv;
if(V.Pos.Z != 0.0)
v->Pos.Z = (THESURFACE->d - V.Pos.X * v->Pos.X - V.Pos.Y * v->Pos.Y)
/ V.Pos.Z;
else
v->Pos.Z = 0.0;
}
void Add_In_Mesh(void *a, void *b)
{
if(!Tree_Search(THEM->Vertices, a))
Tree_Add(THEM->Vertices, a);
}
void Freeze_Vertex(void *a, void *b)
{
Vertex *pv;
pv = *(Vertex **) a;
pv->Freeze.X = pv->Pos.X;
pv->Freeze.Y = pv->Pos.Y;
pv->Freeze.Z = pv->Pos.Z;
pv->Frozen = 1;
}
void deFreeze_Vertex(void *a, void *b)
{
Vertex *pv;
pv = *(Vertex **) a;
if(!pv->Frozen)
return;
pv->Pos.X = pv->Freeze.X;
pv->Pos.Y = pv->Freeze.Y;
pv->Pos.Z = pv->Freeze.Z;
pv->Frozen = 0;
}
void PutVertex_OnSurf(void *a, void *b)
{ Vertex *pv;
pv = *(Vertex **) a;
if(!pv->ListSurf)
pv->ListSurf = List_Create(1, 1, sizeof(Surface *));
List_Add(pv->ListSurf, &THESURFACE);
}
/* remplis la structure Delaunay d'un triangle cree */
Delaunay *testconv(avlptr *root, int *conv, DocRecord * ptr)
{
double qual;
Delaunay *pdel = NULL;
*conv = 0;
if(*root == NULL)
*conv = 1;
else {
pdel = findrightest(*root, ptr);
qual = pdel->t.quality_value;
if(qual < CONV_VALUE)
*conv = 1;
}
return pdel;
}
int mesh_domain(ContourPeek * ListContours, int numcontours,
maillage * mai, int *numpoints, int OnlyTheInitialMesh)
{
List_T *kill_L, *del_L;
MPoint pt;
DocRecord docm, *doc;
int conv, i, j, nump, numact, numlink, numkil, numaloc, numtri, numtrr,
numtrwait;
Delaunay *del, *del_P, *deladd;
PointNumero aa, bb, cc, a, b, c;
DListPeek list, p, q;
avlptr *root, *root_w;
double volume_old, volume_new;
root = (avlptr *) Malloc(sizeof(avlptr));
root_w = (avlptr *) Malloc(sizeof(avlptr));
nump = 0;
for(i = 0; i < numcontours; i++)
nump += ListContours[i]->numpoints;
/* creation du maillage initial grace au puissant algorithme "divide and conquer" */
doc = &docm;
FGMESH = &docm;
doc->points = (PointRecord *) Malloc((nump + 100) * sizeof(PointRecord));
gPointArray = doc->points;
numaloc = nump;
doc->delaunay = 0;
doc->numPoints = nump;
doc->numTriangles = 0;
mai->numtriangles = 0;
mai->numpoints = 0;
if(!numcontours) {
Msg(GERROR, "No contour");
return -1;
}
numact = 0;
for(i = 0; i < numcontours; i++) { for(j = 0; j < ListContours[i]->numpoints; j++) {
doc->points[numact + j].where.h =
ListContours[i]->oriented_points[j].where.h
+ ListContours[i]->perturbations[j].h;
doc->points[numact + j].where.v =
ListContours[i]->oriented_points[j].where.v
+ ListContours[i]->perturbations[j].v;
doc->points[numact + j].quality =
ListContours[i]->oriented_points[j].quality;
doc->points[numact + j].initial =
ListContours[i]->oriented_points[j].initial;
ListContours[i]->oriented_points[j].permu = numact + j;
doc->points[numact + j].permu = numact + j;
doc->points[numact + j].numcontour = ListContours[i]->numerocontour;
doc->points[numact + j].adjacent = NULL;
}
numact += ListContours[i]->numpoints;
}
DelaunayAndVoronoi(doc);
Conversion(doc);
remove_all_dlist(doc->numPoints, doc->points);
if(!is_3D || CTX.mesh.constrained_bgmesh)
BGMESH = doc;
else
BGMESH = NULL;
InitBricks(BGMESH);
/* elimination des triangles exterieurs + verification de l'existence
des edges (coherence) */
makepermut(nump);
del_L = List_Create(doc->numTriangles, 1000, sizeof(Delaunay *));
for(i = 0; i < doc->numTriangles; i++) {
del_P = &doc->delaunay[i];
if((del_P->t.a == del_P->t.b) && (del_P->t.a == del_P->t.c)) {
Msg(GERROR, "Initial mesh is wrong. Try the new isotropic algorithm.");
}
else
List_Add(del_L, &del_P);
}
doc->numTriangles = List_Nbr(del_L);
verify_edges(del_L, ListContours, numcontours);
verify_inside(doc->delaunay, doc->numTriangles);
/* creation des liens ( triangles voisins ) */
CreateLinks(del_L, doc->numTriangles, ListContours, numcontours);
/* initialisation de l'arbre */
(*root) = (*root_w) = NULL;
int onlyinit = OnlyTheInitialMesh;
if(doc->numTriangles == 1) {
doc->delaunay[0].t.position = INTERN;
Msg(INFO, "Only one triangle in initial mesh (skipping refinement)");
onlyinit = 1;
}
for(i = 0; i < doc->numTriangles; i++) {
if(doc->delaunay[i].t.position != EXTERN) {
del = &doc->delaunay[i];
Insert_Triangle(root, del);
}
}
/* maillage proprement dit :
1) Les triangles sont tries dans l'arbre suivant leur qualite
2) on ajoute un noeud au c.g du plus mauvais
3) on reconstruit le delaunay par Bowyer-Watson
4) on ajoute les triangles dans l'arbre et on recommence
jusque quand tous les triangles sont OK */
del = testconv(root, &conv, doc);
PushgPointArray(doc->points);
List_Reset(del_L);
kill_L = List_Create(1000, 1000, sizeof(Delaunay *));
if(onlyinit)
conv = 1;
while(conv != 1) {
pt = Localize(del, doc);
numlink = 0;
numkil = 0;
list = NULL;
if(!PE_Del_Triangle(del, pt, &list, kill_L, del_L, &numlink, &numkil)) {
Msg(WARNING, "Triangle non deleted");
Delete_Triangle(root, del);
Delete_Triangle(root_w, del);
del->t.quality_value /= 10.;
Insert_Triangle(root, del);
numlink = numkil = 0;
if(list != NULL) {
p = list;
do {
q = p;
p = Pred(p);
Free(q);
}
while(p != list);
}
list = NULL;
del = testconv(root, &conv, doc);
continue;
}
/* Test du Volume */
i = 0;
p = list;
volume_new = 0.0;
do {
q = p->next;
bb = p->point_num;
cc = q->point_num;
volume_new += fabs((doc->points[bb].where.h - pt.h) *
(doc->points[cc].where.v - doc->points[bb].where.v) -
(doc->points[cc].where.h - doc->points[bb].where.h) *
(doc->points[bb].where.v - pt.v));
p = q;
} while(p != list);
volume_old = 0.0;
for(i = 0; i < numkil; i++) {
del_P = *(Delaunay **) List_Pointer(kill_L, i);
aa = del_P->t.a;
bb = del_P->t.b;
cc = del_P->t.c;
volume_old += fabs((doc->points[bb].where.h - doc->points[aa].where.h) * (doc->points[cc].where.v - doc->points[bb].where.v) -
(doc->points[cc].where.h - doc->points[bb].where.h) *
(doc->points[bb].where.v - doc->points[aa].where.v));
}
if((volume_old - volume_new) / (volume_old + volume_new) > 1.e-6) {
Msg(WARNING, "Volume has changed (%g->%g)", volume_old, volume_new);
Delete_Triangle(root, del);
Delete_Triangle(root_w, del);
del->t.quality_value /= 10.;
Insert_Triangle(root, del);
numlink = numkil = 0;
if(list != NULL) {
p = list;
do {
q = p;
p = Pred(p);
Free(q);
} while(p != list);
}
list = NULL;
del = testconv(root, &conv, doc);
continue;
}
/* Fin test du volume */
for(i = 0; i < numkil; i++) {
del_P = *(Delaunay **) List_Pointer(kill_L, i);
Delete_Triangle(root, del_P);
// MEMORY_LEAK -JF
// Free(del_P);
}
*numpoints = doc->numPoints;
Insert_Point(pt, numpoints, &numaloc, doc, BGMESH, is_3D);
doc->points = gPointArray;
doc->numPoints = *numpoints;
i = 0;
p = list;
do {
q = p->next;
aa = doc->numPoints - 1;
bb = p->point_num;
cc = q->point_num;
deladd = (Delaunay *) Malloc(sizeof(Delaunay));
filldel(deladd, aa, bb, cc, doc->points, BGMESH);
List_Put(del_L, numlink + i, &deladd);
i++;
p = q;
} while(p != list);
p = list;
do {
q = p;
p = Pred(p);
Free(q);
} while(p != list);
CreateLinks(del_L, i + numlink, ListContours, numcontours);
for(j = 0; j < i; j++) { del_P = *(Delaunay **) List_Pointer(del_L, j + numlink);
Insert_Triangle(root, del_P);
}
del = testconv(root, &conv, doc);
}
List_Delete(kill_L);
numtrwait = 0;
if(*root_w != NULL)
avltree_count(*root_w, &numtrwait);
numtrr = 0;
avltree_count(*root, &numtrr);
mai->listdel = (delpeek *) Malloc((numtrr + numtrwait + 1) * sizeof(delpeek));
numtri = 0;
avltree_print(*root, mai->listdel, &numtri);
if(numtrwait != 0) {
numtri = 0;
avltree_print(*root_w, (Delaunay **) del_L->array, &numtri);
for(i = 0; i < numtrwait; i++) {
mai->listdel[i + numtrr] = *(Delaunay **) List_Pointer(del_L, i);
}
avltree_remove(root_w);
}
avltree_remove(root);
List_Delete(del_L);
mai->numtriangles = numtrr + numtrwait;
mai->numpoints = doc->numPoints;
mai->points = doc->points;
/* Tous Les Triangles doivent etre orientes */
if(!mai->numtriangles) {
Msg(GERROR, "No triangles in surface mesh");
return 0;
//mai->numtriangles = 1;
}
for(i = 0; i < mai->numtriangles; i++) {
a = mai->listdel[i]->t.a;
b = mai->listdel[i]->t.b;
c = mai->listdel[i]->t.c;
mai->listdel[i]->v.voisin1 = NULL;
if(((doc->points[b].where.h - doc->points[a].where.h) *
(doc->points[c].where.v - doc->points[b].where.v) -
(doc->points[c].where.h - doc->points[b].where.h) *
(doc->points[b].where.v - doc->points[a].where.v)) > 0.0) {
mai->listdel[i]->t.a = b;
mai->listdel[i]->t.b = a;
}
}
// MEMORY LEAK (JF)
// Free(doc->delaunay);
// doc->delaunay = 0;
return 1;
}
void Maillage_Automatique_VieuxCode(Surface * pS, Mesh * m, int ori)
{
ContourRecord *cp, **liste;
List_T *c;
Vertex *v, V[3], *ver[3], **pp[3]; maillage M;
int err, i, j, k, N, a, b, d;
Simplex *s;
double Xmin = 0., Xmax = 0., Ymin = 0., Ymax = 0.;
if(m->BGM.Typ == WITHPOINTS) {
is_3D = 0;
}
else {
is_3D = 1;
}
liste =
(ContourPeek *) Malloc(List_Nbr(pS->Contours) * sizeof(ContourPeek));
k = 0;
for(i = 0; i < List_Nbr(pS->Contours); i++) {
cp = (ContourRecord *) Malloc(sizeof(ContourRecord));
List_Read(pS->Contours, i, &c);
cp->oriented_points =
(PointRecord *) Malloc(List_Nbr(c) * sizeof(PointRecord));
cp->perturbations = (MPoint *) Malloc(List_Nbr(c) * sizeof(MPoint));
cp->numerocontour = i;
for(j = 0; j < List_Nbr(c); j++) {
List_Read(c, j, &v);
if(!j) {
Xmin = Xmax = v->Pos.X;
Ymin = Ymax = v->Pos.Y;
}
else {
Xmin = DMIN(Xmin, v->Pos.X);
Xmax = DMAX(Xmax, v->Pos.X);
Ymin = DMIN(Ymin, v->Pos.Y);
Ymax = DMAX(Ymax, v->Pos.Y);
}
}
LC2D = sqrt(DSQR(Xmax - Xmin) + DSQR(Ymax - Ymin));
for(j = 0; j < List_Nbr(c); j++) {
List_Read(c, j, &v);
cp->oriented_points[j].where.h = v->Pos.X;
cp->oriented_points[j].where.v = v->Pos.Y;
// On pourrait imaginer diviser les perturbations par un facteur
// dependant de v->lc, mais ca n'a pas l'air d'ameliorer le bignou
cp->perturbations[j].h = CTX.mesh.rand_factor * LC2D *
(double)rand() / (double)RAND_MAX;
cp->perturbations[j].v = CTX.mesh.rand_factor * LC2D *
(double)rand() / (double)RAND_MAX;
cp->oriented_points[j].numcontour = i;
cp->oriented_points[j].quality = v->lc;
cp->oriented_points[j].permu = k++;
cp->oriented_points[j].initial = v->Num;
}
cp->numpoints = List_Nbr(c);
liste[i] = cp;
}
int res_mesh_domain = 1;
if(pS->Method)
res_mesh_domain = mesh_domain(liste, List_Nbr(pS->Contours), &M, &N,
(CTX.mesh.initial_only == 2));
for(i = 0; i < M.numpoints; i++) { if(gPointArray[i].initial < 0) {
gPointArray[i].initial = ++THEM->MaxPointNum;
v = Create_Vertex(gPointArray[i].initial, gPointArray[i].where.h,
gPointArray[i].where.v, 0.0, gPointArray[i].quality,
0.0);
if(!Tree_Search(pS->Vertices, &v))
Tree_Add(pS->Vertices, &v);
}
}
for(i = 0; i < 3; i++)
ver[i] = &V[i];
for(i = 0; i < M.numtriangles; i++) {
a = M.listdel[i]->t.a;
b = M.listdel[i]->t.b;
d = M.listdel[i]->t.c;
ver[0]->Num = gPointArray[a].initial;
ver[1]->Num = gPointArray[b].initial;
ver[2]->Num = gPointArray[d].initial;
err = 0;
for(j = 0; j < 3; j++) {
if(!(pp[j] = (Vertex **) Tree_PQuery(pS->Vertices, &ver[j]))) {
err = 1;
Msg(GERROR, "Unknown vertex %d", ver[j]->Num);
}
}
/*
Je n'ai pas l'impression que ceci fonctionne comme voulu (essaie
e.g. avec demos/machine.geo)...
if (ori && !err)
s = Create_Simplex (*pp[0], *pp[1], *pp[2], NULL);
else if (!err)
s = Create_Simplex (*pp[0], *pp[2], *pp[1], NULL);
if (!err){
s->iEnt = pS->Num;
Tree_Insert (pS->Simplexes, &s);
}
*/
if(!err) {
s = Create_Simplex(*pp[0], *pp[2], *pp[1], NULL);
s->iEnt = pS->Num;
Tree_Insert(pS->Simplexes, &s);
}
// MEMORY LEAK (JF)
//printf("%d %d %d\n", M.listdel[i]->t.a, M.listdel[i]->t.b, M.listdel[i]->t.c);
//xxxxFree(M.listdel[i]);
}
if(res_mesh_domain >= 0)
Free(M.listdel);
Free(gPointArray);
}
void Make_Mesh_With_Points(DocRecord * ptr, PointRecord * Liste,
int Numpoints)
{
ptr->numTriangles = 0;
ptr->points = Liste;
ptr->numPoints = Numpoints;
ptr->delaunay = 0;
DelaunayAndVoronoi(ptr);
Conversion(ptr);
remove_all_dlist(ptr->numPoints, ptr->points);
}
void filldel(Delaunay * deladd, int aa, int bb, int cc,
PointRecord * points, DocRecord * mesh)
{
double qual, newqual;
MPoint pt2;
Vertex *v, *dum;
deladd->t.a = aa;
deladd->t.b = bb;
deladd->t.c = cc;
deladd->t.info = TOLINK;
deladd->t.info2 = 0;
deladd->v.voisin1 = NULL;
deladd->v.voisin2 = NULL;
deladd->v.voisin3 = NULL;
CircumCircle(points[aa].where.h, points[aa].where.v,
points[bb].where.h, points[bb].where.v,
points[cc].where.h, points[cc].where.v,
&deladd->t.xc, &deladd->t.yc);
pt2.h = deladd->t.xc;
pt2.v = deladd->t.yc;
if(!is_3D) {
if(mesh) {
newqual = find_quality(pt2, mesh);
}
else {
newqual =
(points[aa].quality + points[bb].quality + points[cc].quality) / 3.;
}
v = Create_Vertex(-1, pt2.h, pt2.v, 0.0, 0.0, 0.0);
Calcule_Z_Plan(&v, &dum);
Projette_Inverse(&v, &dum);
Free_Vertex(&v, 0);
}
else {
v = Create_Vertex(-1, pt2.h, pt2.v, 0.0, 0.0, 0.0);
Calcule_Z_Plan(&v, &dum);
Projette_Inverse(&v, &dum);
qual = Lc_XYZ(v->Pos.X, v->Pos.Y, v->Pos.Z, THEM);
if(CTX.mesh.constrained_bgmesh) {
if(mesh) {
newqual = MIN(qual, find_quality(pt2, mesh));
}
else {
newqual =
MIN(qual,
(points[aa].quality + points[bb].quality +
points[cc].quality) / 3.);
}
}
else
newqual = qual;
Free_Vertex(&v, 0);
}
deladd->t.quality_value =
sqrt((deladd->t.xc - points[cc].where.h) * (deladd->t.xc -
points[cc].where.h) +
(deladd->t.yc - points[cc].where.v) * (deladd->t.yc -
points[cc].where.v)
) / newqual;
deladd->t.position = INTERN;
}
void ActionEndTheCurve(void *a, void *b)
{
Curve *c = *(Curve **) a; End_Curve(c);
}
void ActionInvertTri(void *a, void *b)
{
Simplex *s = *(Simplex **) a;
Vertex *tmp = s->V[1];
s->V[1] = s->V[2];
s->V[2] = tmp;
}
void ActionInvertQua(void *a, void *b)
{
Quadrangle *q = *(Quadrangle **) a;
Vertex *tmp = q->V[1];
q->V[1] = q->V[3];
q->V[3] = tmp;
}
// A little routine that re-orients the surface meshes so that they
// match the geometrical orientations. It would be better to generate
// the meshes correctly in the first place (!), but it's actually
// pretty nice to do it here since we use many algorithms (some
// closed) and since the mean plane computation doesn't take the
// original orientation into account.
void ReOrientSurfaceMesh(Surface *s)
{
if(!List_Nbr(s->Generatrices))
return;
Curve *c;
List_Read(s->Generatrices, 0, &c);
if(List_Nbr(c->Vertices) < 2)
return;
// take the first edge (even if there are holes, we know that this
// edge belongs to the exterior loop, due to the ordering of
// s->Generatrices)
Vertex *v1, *v2;
List_Read(c->Vertices, 0, &v1);
List_Read(c->Vertices, 1, &v2);
Edge edge;
edge.V[0] = v1;
edge.V[1] = v2;
if(Tree_Nbr(s->Simplexes)){
EdgesContainer edges;
edges.AddSimplexTree(s->Simplexes);
Edge *e = (Edge*)Tree_PQuery(edges.AllEdges, &edge);
if(e && e->V[0]->Num == v2->Num && e->V[1]->Num == v1->Num){
Tree_Action(s->Simplexes, ActionInvertTri);
if(Tree_Nbr(s->Quadrangles)) // for partially recombined meshes
Tree_Action(s->Quadrangles, ActionInvertQua);
}
}
if(Tree_Nbr(s->Quadrangles)){
EdgesContainer edges;
edges.AddQuadrangleTree(s->Quadrangles);
Edge *e = (Edge*)Tree_PQuery(edges.AllEdges, &edge);
if(e && e->V[0]->Num == v2->Num && e->V[1]->Num == v1->Num){
Tree_Action(s->Quadrangles, ActionInvertQua);
if(Tree_Nbr(s->Simplexes)) // for partially recombined meshes
Tree_Action(s->Simplexes, ActionInvertTri);
}
}
}
void Maillage_Surface(void *data, void *dum)
{
Tree_T *tnxe;
int ori;
Surface *s = *(Surface **) data;
if(!s->Support)
return;
Msg(STATUS3, "Meshing surface %d", s->Num);
int tag = MeshDiscreteSurface(s);
if(tag == 1){
Tree_Action(THEM->Points, PutVertex_OnSurf);
Tree_Action(s->Vertices, PutVertex_OnSurf);
Tree_Action(s->Vertices, Add_In_Mesh);
return;
}
else if (tag == 2)
{
return;
}
THESUPPORT = s->Support;
THESURFACE = s;
if(Tree_Nbr(s->Simplexes))
Tree_Delete(s->Simplexes);
s->Simplexes = Tree_Create(sizeof(Simplex *), compareQuality);
if(Tree_Nbr(s->Vertices))
Tree_Delete(s->Vertices);
s->Vertices = Tree_Create(sizeof(Vertex *), compareVertex);
if(MeshTransfiniteSurface(s) ||
MeshEllipticSurface(s) ||
MeshCylindricalSurface(s) ||
MeshParametricSurface(s) ||
Extrude_Mesh(s)) {
Tree_Action(THEM->Points, PutVertex_OnSurf);
Tree_Action(s->Vertices, PutVertex_OnSurf);
Tree_Action(s->Vertices, Add_In_Mesh);
}
else{
if(!List_Nbr(s->Generatrices)){
Msg(GERROR, "Cannot mesh surface with no boundary");
return;
}
int TypSurface = s->Typ;
Plan_Moyen(&s, dum);
Tree_Action(THEM->Points, Freeze_Vertex);
Tree_Action(s->Vertices, Freeze_Vertex);
Tree_Action(THEM->Points, Projette_Plan_Moyen);
Tree_Action(s->Vertices, Projette_Plan_Moyen);
Tree_Action(THEM->Curves, ActionEndTheCurve);
s->Typ = MSH_SURF_PLAN;
End_Surface(s, 0);
ori = Calcule_Contours(s);
if(CTX.mesh.algo2d == DELAUNAY_ISO)
Maillage_Automatique_VieuxCode(s, THEM, ori);
else if(CTX.mesh.algo2d == DELAUNAY_ANISO)
AlgorithmeMaillage2DAnisotropeModeJF(s);
else
Mesh_Triangle(s);
if(CTX.mesh.nb_smoothing) {
Msg(STATUS3, "Smoothing surface %d", s->Num);
tnxe = Tree_Create(sizeof(NXE), compareNXE);
create_NXE(s->Vertices, s->Simplexes, tnxe);
for(int i = 0; i < CTX.mesh.nb_smoothing; i++)
Tree_Action(tnxe, ActionLiss);
delete_NXE(tnxe);
}
if(s->Recombine){
Msg(STATUS3, "Recombining surface %d", s->Num);
Recombine(s->Vertices, s->Simplexes, s->Quadrangles, s->RecombineAngle);
}
s->Typ = TypSurface;
if(s->Typ != MSH_SURF_PLAN) {
if(s->Extrude)
s->Extrude->Rotate(s->plan);
Tree_Action(s->Vertices, Calcule_Z);
if(s->Extrude)
s->Extrude->Rotate(s->invplan);
}
else
Tree_Action(s->Vertices, Calcule_Z_Plan);
Tree_Action(s->Vertices, Projette_Inverse);
Tree_Action(THEM->Points, Projette_Inverse);
Tree_Action(THEM->Points, deFreeze_Vertex);
Tree_Action(s->Vertices, deFreeze_Vertex);
Tree_Action(THEM->Points, PutVertex_OnSurf);
Tree_Action(s->Vertices, PutVertex_OnSurf);
Tree_Action(s->Vertices, Add_In_Mesh);
Tree_Action(THEM->Curves, ActionEndTheCurve);
End_Surface(s->Support, 0);
End_Surface(s, 0);
ReOrientSurfaceMesh(s);
}
}