1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
|
#include "stdio.h"
#include "mesh.h"
#define SMOOTH_ITERATIONS 20
/*********
this routine will smooth out the mesh
it will only move points whose label is 0
*************/
/* relaxing step */
#define STEP(CEN,WEST,EST,SOUTH,NORTH) (((NORTH)+(SOUTH)+(EST)+(WEST)) / 4.0)
void smooth_mesh(MeshT *mesh,
int preserve
/* if true, a point cannot enter a neighbouring cell
but the algo is far from perfect in this respect*/
)
{
int xi, yi;
double x,y;
int lp;
for(lp =SMOOTH_ITERATIONS ; lp ; lp--) {
for(xi=0; xi < mesh->nx ; xi++) {
for(yi=0; yi<mesh->ny ; yi++) {
if( 0 == meshGetLabel(mesh,xi,yi)) {
if ( xi == 0 || xi == mesh->nx -1)
x=meshGetx(mesh, xi,yi);
else
{
x=STEP(meshGetx(mesh, xi,yi),
meshGetxRefl(mesh, xi+1,yi) ,
meshGetxRefl(mesh, xi,yi+1) ,
meshGetxRefl(mesh, xi-1,yi) ,
meshGetxRefl(mesh, xi,yi-1)) ;
if(preserve)
{
if ( x< meshGetxRefl(mesh, xi-1,yi))
x= meshGetxRefl(mesh, xi-1,yi);
else
if(x> meshGetxRefl(mesh, xi+1,yi))
x= meshGetxRefl(mesh, xi+1,yi);
else
if ( x< meshGetxRefl(mesh, xi-1,yi+1))
x= meshGetxRefl(mesh, xi-1,yi+1);
else
if(x> meshGetxRefl(mesh, xi+1,yi-1))
x= meshGetxRefl(mesh, xi+1,yi-1);
}
}
if ( yi == 0 || yi == mesh->ny -1)
y=meshGety(mesh, xi,yi);
else
{
y=STEP(meshGetx(mesh, xi,yi),
meshGetyRefl(mesh, xi+1,yi),
meshGetyRefl(mesh, xi,yi+1),
meshGetyRefl(mesh, xi-1,yi),
meshGetyRefl(mesh, xi,yi-1)) ;
if(preserve)
{
if ( y< meshGetyRefl(mesh, xi,yi-1))
y= meshGetyRefl(mesh, xi,yi-1);
else
if(y> meshGetyRefl(mesh, xi,yi+1))
y= meshGetyRefl(mesh, xi,yi+1);
else
if ( y< meshGetyRefl(mesh, xi-1,yi-1))
y= meshGetyRefl(mesh, xi-1,yi-1);
else
if(y> meshGetyRefl(mesh, xi+1,yi+1))
y= meshGetyRefl(mesh, xi+1,yi+1);
}
}
meshSetNoundo(mesh,xi,yi,x,y);
}
}
}
}
}
|
