1 |
|
1 |
/*****************************************************************
|
2 |
|
|
|
3 |
|
|
m1d_lissage.c Type:Func
|
4 |
|
|
|
5 |
|
|
Lissage du maillage lineique
|
6 |
|
|
|
7 |
|
|
Date de creation : Thu Dec 5 10:23:03 1996
|
8 |
|
|
|
9 |
|
|
Derniere version : Fri May 9 16:05:11 1997
|
10 |
|
|
|
11 |
|
|
|
12 |
|
|
|
13 |
|
|
|
14 |
|
|
|
15 |
|
|
|
16 |
|
|
|
17 |
|
|
|
18 |
|
|
|
19 |
|
|
Vincent FRANCOIS
|
20 |
|
|
|
21 |
|
|
*****************************************************************/
|
22 |
|
|
|
23 |
|
|
|
24 |
|
|
|
25 |
|
|
|
26 |
|
|
|
27 |
|
|
/**************************/
|
28 |
|
|
/* include */
|
29 |
|
|
#include <stdio.h>
|
30 |
|
|
#include <string.h>
|
31 |
|
|
#include <math.h>
|
32 |
|
|
#include <stdlib.h>
|
33 |
|
|
#include "const.h"
|
34 |
|
|
#include "struct.h"
|
35 |
|
|
#include "memoire.h"
|
36 |
|
|
#include "prototype.h"
|
37 |
|
|
|
38 |
|
|
|
39 |
|
|
/**************************/
|
40 |
|
|
/* variables globales */
|
41 |
|
|
extern struct s_param *para;
|
42 |
|
|
extern struct environnement env;
|
43 |
|
|
extern struct s_acis *acis;
|
44 |
|
|
extern struct s_mesh *mesh;
|
45 |
|
|
|
46 |
|
|
|
47 |
|
|
#define FDN2(s) (dnm*(s-s1)*(s-s2)/((sm-s1)*(sm-s2))+dn1*(s-sm)*(s-s2)/((s1-sm)*(s1-s2))+dn2*(s-s1)*(s-sm)/((s2-s1)*(s2-sm)))
|
48 |
|
|
|
49 |
|
|
/**************************/
|
50 |
|
|
/* programme principal */
|
51 |
|
|
|
52 |
|
|
void m1d_lissage(int segment1,int segment2)
|
53 |
|
|
{
|
54 |
|
|
struct s_edge *edge;
|
55 |
|
|
struct s_segment *seg,*sega;
|
56 |
|
|
struct s_face *face;
|
57 |
|
|
struct s_noeud *no,*no1,*no2,*no1a,*no2a;
|
58 |
|
|
int i,j,num,degre,nfin;
|
59 |
|
|
float t,a,dens,l,vec[4],coord[4],vec1[4],vec2[4];
|
60 |
|
|
float ti,tii,aprec,tprec,limit,dn1,dn2,dnm,sm,s1,s2,s,dis,t1,t2;
|
61 |
|
|
int n1,n2,n,noe,numo,raffine;
|
62 |
|
|
float *densite;
|
63 |
|
|
int *iordre;
|
64 |
|
|
|
65 |
|
|
|
66 |
|
|
/* moyenne des longueurs */
|
67 |
|
|
/*moy=0.;
|
68 |
|
|
nb=0;
|
69 |
|
|
for (i=0;i<mesh->nb_noeud;i++)
|
70 |
|
|
{
|
71 |
|
|
no=ADRESSE(i,noeud,mesh->);
|
72 |
|
|
for (j=0;j<no->nb_segment;j++)
|
73 |
|
|
{
|
74 |
|
|
valeur=no->segment[j]->longueur;
|
75 |
|
|
if (valeur>env.dens) valeur=env.dens;
|
76 |
|
|
moy=moy+valeur;
|
77 |
|
|
nb++;
|
78 |
|
|
}
|
79 |
|
|
}
|
80 |
|
|
moy=moy/nb; */
|
81 |
|
|
/* ecart type */
|
82 |
|
|
/*ec=0.;
|
83 |
|
|
nb=0;
|
84 |
|
|
for (i=0;i<mesh->nb_noeud;i++)
|
85 |
|
|
{
|
86 |
|
|
no=ADRESSE(i,noeud,mesh->);
|
87 |
|
|
for (j=0;j<no->nb_segment;j++)
|
88 |
|
|
{
|
89 |
|
|
valeur=no->segment[j]->longueur;
|
90 |
|
|
if (valeur>env.dens) valeur=env.dens;
|
91 |
|
|
ec=ec+(valeur-moy)*(valeur-moy);
|
92 |
|
|
nb++;
|
93 |
|
|
}
|
94 |
|
|
}
|
95 |
|
|
ec=(float)sqrt((double)(ec/(nb-1)));*/
|
96 |
|
|
/* calcul de densite */
|
97 |
|
|
densite=(float *)calloc(mesh->nb_noeud+1,sizeof(float));
|
98 |
|
|
ERREUR_ALLOC(densite);
|
99 |
|
|
iordre=(int *)calloc(mesh->nb_noeud+1,sizeof(int));
|
100 |
|
|
ERREUR_ALLOC(iordre);
|
101 |
|
|
for (i=0;i<mesh->nb_noeud;i++)
|
102 |
|
|
{
|
103 |
|
|
no=ADRESSE(i,noeud,mesh->);
|
104 |
|
|
dens=env.dens;
|
105 |
|
|
no->dens=dens;
|
106 |
|
|
for (j=0;j<no->nb_segment;j++)
|
107 |
|
|
// if (no->segment[j]->longueur>moy-ec)
|
108 |
|
|
MINI(no->dens,no->dens,no->segment[j]->longueur)
|
109 |
|
|
}
|
110 |
|
|
for (i=segment1;i<segment2;i++)
|
111 |
|
|
{
|
112 |
|
|
voir();
|
113 |
|
|
seg=ADRESSE(i,segment,mesh->);
|
114 |
|
|
edge=(struct s_edge *)acis->entity[seg->num_ent];
|
115 |
|
|
no1=ADRESSE(seg->n1,noeud,mesh->);
|
116 |
|
|
no2=ADRESSE(seg->n2,noeud,mesh->);
|
117 |
|
|
dn1=no1->dens;
|
118 |
|
|
dn2=no2->dens;
|
119 |
|
|
dis=seg->longueur;
|
120 |
|
|
s1=0.;
|
121 |
|
|
s2=dis;
|
122 |
|
|
if (no1->type==EDGE) t1=no1->pos; else t1=edge->t1;
|
123 |
|
|
if (no2->type==EDGE) t2=no2->pos-t1; else t2=edge->t1+edge->t2-t1;
|
124 |
|
|
if (edge->nb_segment!=0)
|
125 |
|
|
{
|
126 |
|
|
sega=edge->segment[0];
|
127 |
|
|
no1a=ADRESSE(sega->n1,noeud,mesh->);
|
128 |
|
|
no2a=ADRESSE(sega->n2,noeud,mesh->);
|
129 |
|
|
VEC(vec1,no1a,no2a);
|
130 |
|
|
t=geo_ind_point_edge(edge,no1a);
|
131 |
|
|
eval_edge(edge,t,DERIVE,vec2);
|
132 |
|
|
if (PSCA(vec1,vec2)<0.)
|
133 |
|
|
{
|
134 |
|
|
if (no2->type==EDGE) t1=no2->pos; else t1=edge->t1;
|
135 |
|
|
if (no1->type==EDGE) t2=no1->pos-t1; else t2=edge->t1+edge->t2-t1;
|
136 |
|
|
}
|
137 |
|
|
}
|
138 |
|
|
raffine=0;
|
139 |
|
|
face=edge->owning_coedge->owning_loop->owning_face;
|
140 |
|
|
if (strcmp(acis->type_entite[face->surface],"plane-surface")!=0) raffine=1;
|
141 |
|
|
if (edge->owning_coedge->partner_coedge!=NULL)
|
142 |
|
|
{
|
143 |
|
|
face=edge->owning_coedge->partner_coedge->owning_loop->owning_face;
|
144 |
|
|
if (strcmp(acis->type_entite[face->surface],"plane-surface")!=0) raffine=1;
|
145 |
|
|
}
|
146 |
|
|
if (raffine==1)
|
147 |
|
|
{
|
148 |
|
|
dn1=edge->start_vertex->noeud->dens;
|
149 |
|
|
dn2=edge->end_vertex->noeud->dens;
|
150 |
|
|
}
|
151 |
|
|
if ( (dis<1.5*env.dens) || (raffine==1) )
|
152 |
|
|
if (dn1==dn2) degre=0;
|
153 |
|
|
else degre=1;
|
154 |
|
|
else
|
155 |
|
|
if ((dn1==dn2) && (dn1==env.dens)) degre=0;
|
156 |
|
|
else if ( (dn1!=env.dens) && (dn2==env.dens)) degre=1;
|
157 |
|
|
else if ( (dn1==env.dens) && (dn2!=env.dens)) degre=1;
|
158 |
|
|
else degre=2;
|
159 |
|
|
if (degre==2)
|
160 |
|
|
{
|
161 |
|
|
dnm=env.dens;
|
162 |
|
|
sm=(dn1*s1+dn2*s2)/(dn1+dn2);
|
163 |
|
|
}
|
164 |
|
|
else
|
165 |
|
|
{
|
166 |
|
|
dnm=0.5*(dn1+dn2);
|
167 |
|
|
sm=dis/2;
|
168 |
|
|
}
|
169 |
|
|
MINI(l,dn1,dn2);
|
170 |
|
|
num=(int)floor((double)10*dis/l)+1;
|
171 |
|
|
a=0.;
|
172 |
|
|
for (j=0;j<num;j++)
|
173 |
|
|
{
|
174 |
|
|
ti=(t1+t2*j/num);
|
175 |
|
|
tii=(t1+t2*(j+1)/num);
|
176 |
|
|
t=(0.7886751345*ti+0.2113248654*tii);
|
177 |
|
|
eval_edge(edge,t,DERIVE,coord);
|
178 |
|
|
s=(t-t1)/t2*dis;
|
179 |
|
|
a=a+(float)sqrt(coord[0]*coord[0]+coord[1]*coord[1]+coord[2]*coord[2])/FDN2(s);
|
180 |
|
|
t=(0.7886751345*tii+0.2113248654*ti);
|
181 |
|
|
eval_edge(edge,t,DERIVE,coord);
|
182 |
|
|
s=(t-t1)/t2*dis;
|
183 |
|
|
a=a+(float)sqrt(coord[0]*coord[0]+coord[1]*coord[1]+coord[2]*coord[2])/FDN2(s);
|
184 |
|
|
}
|
185 |
|
|
a=(float)fabs((double)(a*t2/num/2));
|
186 |
|
|
n=(int)floor((double)a);
|
187 |
|
|
if (a-floor((double)a)>0.5) n++;
|
188 |
|
|
if (n<1) n=1;
|
189 |
|
|
limit=a/(float)n;
|
190 |
|
|
n1=seg->n1;
|
191 |
|
|
nfin=seg->n2;
|
192 |
|
|
j=0;
|
193 |
|
|
a=0.;
|
194 |
|
|
if (n!=1)
|
195 |
|
|
{
|
196 |
|
|
for (noe=0;noe<n;noe++)
|
197 |
|
|
{
|
198 |
|
|
if (noe!=n-1)
|
199 |
|
|
{
|
200 |
|
|
while (a<limit)
|
201 |
|
|
{
|
202 |
|
|
aprec=a;
|
203 |
|
|
tprec=tii;
|
204 |
|
|
ti=(t1+t2*j/num);
|
205 |
|
|
tii=(t1+t2*(j+1)/num);
|
206 |
|
|
t=(0.7886751345*ti+0.2113248654*tii);
|
207 |
|
|
eval_edge(edge,t,DERIVE,coord);
|
208 |
|
|
s=(t-t1)/t2*dis;
|
209 |
|
|
a=a+(float)fabs((double)t2)/num/2*(float)sqrt(coord[0]*coord[0]+coord[1]*coord[1]+coord[2]*coord[2])/FDN2(s);
|
210 |
|
|
t=(0.7886751345*tii+0.2113248654*ti);
|
211 |
|
|
eval_edge(edge,t,DERIVE,coord);
|
212 |
|
|
s=(t-t1)/t2*dis;
|
213 |
|
|
a=a+(float)fabs((double)t2)/num/2*(float)sqrt(coord[0]*coord[0]+coord[1]*coord[1]+coord[2]*coord[2])/FDN2(s);
|
214 |
|
|
j++;
|
215 |
|
|
}
|
216 |
|
|
t=tprec+(tii-tprec)*(limit-aprec)/(a-aprec);
|
217 |
|
|
NEW_ENTITE(no,noeud,mesh->);
|
218 |
|
|
no->num=mesh->nb_noeud-1;
|
219 |
|
|
eval_edge(edge,t,FONCTION,coord);
|
220 |
|
|
no->x=coord[0];
|
221 |
|
|
no->y=coord[1];
|
222 |
|
|
no->z=coord[2];
|
223 |
|
|
no->dens=env.dens;
|
224 |
|
|
no->type=EDGE;
|
225 |
|
|
no->num_ent=edge->num;
|
226 |
|
|
no->pos=t;
|
227 |
|
|
n2=no->num;
|
228 |
|
|
}
|
229 |
|
|
else n2=nfin;
|
230 |
|
|
if (noe==0) seg=ADRESSE(i,segment,mesh->);
|
231 |
|
|
else
|
232 |
|
|
{
|
233 |
|
|
NEW_ENTITE(seg,segment,mesh->);
|
234 |
|
|
seg->num=mesh->nb_segment-1;
|
235 |
|
|
NEW_POINTEUR(numo,segment,edge->);
|
236 |
|
|
edge->segment[numo]=seg;
|
237 |
|
|
}
|
238 |
|
|
seg->n1=n1;
|
239 |
|
|
seg->n2=n2;
|
240 |
|
|
seg->type=EDGE;
|
241 |
|
|
seg->num_ent=edge->num;
|
242 |
|
|
n1=n2;
|
243 |
|
|
a=a-limit;
|
244 |
|
|
}
|
245 |
|
|
}
|
246 |
|
|
}
|
247 |
|
|
/* mise a jour des liens noeud segment et de la densite aux noeuds */
|
248 |
|
|
|
249 |
|
|
for (i=0;i<mesh->nb_noeud;i++)
|
250 |
|
|
{
|
251 |
|
|
no=ADRESSE(i,noeud,mesh->);
|
252 |
|
|
if (no->nb_segment!=0)
|
253 |
|
|
{
|
254 |
|
|
free(no->segment);
|
255 |
|
|
no->segment=NULL;
|
256 |
|
|
no->nb_segment=0;
|
257 |
|
|
}
|
258 |
|
|
}
|
259 |
|
|
for (i=0;i<mesh->nb_segment;i++)
|
260 |
|
|
{
|
261 |
|
|
seg=ADRESSE(i,segment,mesh->);
|
262 |
|
|
no1=ADRESSE(seg->n1,noeud,mesh->);
|
263 |
|
|
no2=ADRESSE(seg->n2,noeud,mesh->);
|
264 |
|
|
VEC(vec,no1,no2);
|
265 |
|
|
seg->longueur=(float)sqrt((double)PSCA(vec,vec));
|
266 |
|
|
NEW_POINTEUR(num,segment,no1->);
|
267 |
|
|
no1->segment[num]=seg;
|
268 |
|
|
NEW_POINTEUR(num,segment,no2->);
|
269 |
|
|
no2->segment[num]=seg;
|
270 |
|
|
}
|
271 |
|
|
for (i=0;i<mesh->nb_noeud;i++)
|
272 |
|
|
{
|
273 |
|
|
no=ADRESSE(i,noeud,mesh->);
|
274 |
|
|
dens=env.dens;
|
275 |
|
|
no->dens=dens;
|
276 |
|
|
for (j=0;j<no->nb_segment;j++)
|
277 |
|
|
MINI(no->dens,no->dens,no->segment[j]->longueur)
|
278 |
|
|
}
|
279 |
|
|
}
|