ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/acismesh/m1d_lissage.cpp
Revision: 1
Committed: Mon Jun 11 22:53:07 2007 UTC (17 years, 11 months ago)
File size: 7782 byte(s)
Log Message:

File Contents

# User Rev Content
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     }