1 subroutine sfseno ( cono, numlig, unst2x, epsid2,
2 > seglig, somseg, geocoo, abscur,
4 c ______________________________________________________________________
8 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
10 c Version originale enregistree le 18 juin 1996 sous le numero 96036
11 c aupres des huissiers de justice Simart et Lavoir a Clamart
12 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
13 c aupres des huissiers de justice
14 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
16 c HOMARD est une marque deposee d'Electricite de France
22 c ______________________________________________________________________
24 c Suivi de Frontiere - SEgment - NOeud
26 c ______________________________________________________________________
28 c but : recherche du segment de la ligne auquel appartient le noeud
29 c ______________________________________________________________________
31 c . nom . e/s . taille . description .
32 c .____________________________________________________________________.
33 c . cono . es . sdim . coordonnees du noeud .
34 c . numlig . e . 1 . numero de ligne sur laquelle on cherche .
35 c . unst2x . e . 1 . inverse de la taille maximale au carre .
36 c . epsid2 . e . 1 . precision relative pour carre de distance .
37 c . seglig . e .0:sfnbli. pointeur dans le tableau somseg : les .
38 c . . . . segments de la ligne i sont aux places de .
39 c . . . . seglig(i-1)+1 a seglig(i)-1 inclus .
40 c . somseg . e . sfnbse . liste des sommets des lignes separees par .
42 c . geocoo . e .sfnbso**. coordonnees des sommets de la frontiere .
43 c . seg . s . 1 . numero de segment trouve .
44 c - : noeud sur extremite de ligne .
45 c + : noeud hors extremites .
46 c . acnoeu . s . 1 . abscisse curviligne du noeud sur la ligne .
47 c ______________________________________________________________________
50 c 0. declarations et dimensionnement
53 c 0.1. ==> generalites
66 double precision unst2x, epsid2
67 double precision geocoo(sfnbso,*), cono(sdim)
68 integer seglig(0:sfnbli), somseg(sfnbse)
70 double precision acnoeu
71 double precision abscur(sfnbse)
73 c 0.4. ==> variables locales
75 integer iaux, jaux, kaux
76 integer seg1, seg2, segmin
78 double precision cooa(3), coob(3)
79 double precision daux, daux1, daux2
88 c 0.5. ==> initialisations
89 c ______________________________________________________________________
99 write(ulsort,*) 'Entree dans SFSENO'
100 write(ulsort,90004) 'unst2x',unst2x
101 write(ulsort,90004) 'epsid2',epsid2
102 write(ulsort,90004) 'cono', (cono(iaux),iaux = 1 , sdim)
103 write(ulsort,90002) 'numlig', numlig
107 c 1.1. ==> segment pas encore trouve
111 c 1.2. ==> bornes de la ligne dans la numerotation des segments
113 seg1 = seglig(numlig-1) + 1
114 seg2 = seglig(numlig ) - 2
115 #ifdef _DEBUG_HOMARD_
116 if ( glop.ne.0 ) then
117 write(ulsort,90006) 'Pointeurs seglig pour la ligne', numlig,
118 >' de',seg1,' a ',seg2+2
119 write(ulsort,90002) '==> Segments extremites',seg1,seg2
120 write(ulsort,90002) '==> Sommets extremites ',
121 > somseg(seg1),somseg(seg2+1)
126 c 2. Le noeud est-il une extremite de la ligne ?
127 c Remarque : on commence par la fin de la ligne pour pouvoir
128 c initialiser correctement l'etape suivante en cas de
130 c Remarque : comme ce sont les coordonnees du debut du segment qui
131 c sont stockees, il faut examiner le segment fictif
132 c (dernier+1) pour trouver les coordonnees de sa fin
137 if ( jaux.eq.1 ) then
144 do 21 , iaux = 1 , sdim
145 coob(iaux) = geocoo(somseg(kaux),iaux)
146 daux = daux + (coob(iaux)-cono(iaux))**2
149 #ifdef _DEBUG_HOMARD_
150 if ( glop.ne.0 ) then
151 if ( jaux.eq.1 ) then
152 write(ulsort,90002) 'segment (jaux=1)', seg2
154 write(ulsort,90002) 'segment (jaux=2)', seg1
156 write(ulsort,90024) 'Noeud', somseg(kaux),
157 > (geocoo(somseg(kaux),iaux),iaux=1,sdim)
158 write(ulsort,90024) '==> carre distance a l''extremite',
162 if ( daux*unst2x.le.epsid2 ) then
164 acnoeu = abscur(kaux)
165 #ifdef _DEBUG_HOMARD_
166 if ( glop.ne.0 ) then
167 write(ulsort,90015) '==> Le noeud est l''extremite',
168 > 1+mod(jaux,2),' de la ligne', numlig
178 c 3. Le noeud n'est pas une extremite de la ligne
179 c On va chercher le segment dont la premiere extremite est
180 c la plus proche du noeud courant.
181 c Si une de ces extremites est le noeud courant, on le note
185 #ifdef _DEBUG_HOMARD_
186 if ( glop.ne.0 ) then
187 write(ulsort,90002) 'Le noeud n''est pas extremite de la ligne',
194 do 30 , kaux = seg1 , seg2
196 #ifdef _DEBUG_HOMARD_
197 if ( glop.ne.0 ) then
198 if ( kaux.eq.seg1 .or. kaux.le.-1 ) then
199 write(ulsort,90006) ' Segment [',
200 > somseg(kaux),' ',somseg(kaux+1),'] '
205 c 3.1. ==> Calcul de la distance au point A, debut du segment
208 do 31 , iaux = 1 , sdim
209 daux = daux + (geocoo(somseg(kaux),iaux)-cono(iaux))**2
211 cgn if ( kaux.eq.seg1 .or. kaux.le.-1 ) then
212 cgn write(ulsort,90004) ' Carre de la distance', daux
215 if ( daux*unst2x.le.epsid2 ) then
216 #ifdef _DEBUG_HOMARD_
217 if ( glop.ne.0 ) then
218 write(ulsort,90004) 'Carre de la distance', daux
220 > '==> Le noeud est au debut du segment', kaux
224 acnoeu = abscur(kaux)
229 c 3.2. ==> Memorisation du minimum
231 cgn write(ulsort,90004) ' Distance ', daux
232 if ( daux.le.daux1 ) then
236 cgn write(ulsort,*) ' Minimum pour le segment ',segmin,
237 cgn > ' : [',somseg(segmin),',',somseg(segmin+1),'] '
243 #ifdef _DEBUG_HOMARD_
244 if ( glop.ne.0 ) then
245 write(ulsort,90004) 'Distance minimale', daux1
246 write(ulsort,90006) ' atteinte sur le segment ',segmin,
247 > '[',somseg(segmin),' ',somseg(segmin+1),'] '
252 c 4. Le noeud n'est pas une extremite d'un segment
253 c On precise les points A et B
256 #ifdef _DEBUG_HOMARD_
257 if ( glop.ne.0 ) then
258 write(ulsort,90002) 'Le premier segment est seg1', seg1
259 write(ulsort,90002) 'Le dernier segment est seg2', seg2
261 > 'Le segment le plus proche du noeud est segmin', segmin
265 c 4.1. ==> On a trouve le segment dont le debut, A, est le sommet le
267 c Si c'est le premier segment, il n'y a pas d'equivoque :
268 c N sera place entre A et B
269 c Sinon, il faut preciser entre les deux segments concernes,
270 c segmin-1 et segmin. On postule que le segment le plus proche
271 c est celui dont la deuxieme extremite est la plus proche de N.
276 c segmin-1 . . segmin
280 if ( segmin.ne.seg1 ) then
284 do 41 , iaux = 1 , sdim
286 > (geocoo(somseg(segmin-1),iaux)-cono(iaux))**2
288 > (geocoo(somseg(segmin+1),iaux)-cono(iaux))**2
291 #ifdef _DEBUG_HOMARD_
292 if ( glop.ne.0 ) then
293 write(ulsort,90024) '(C-A) Carre de la distance au debut de',
295 write(ulsort,90024) '(A-B) Carre de la distance au debut de',
297 write(ulsort,90024) '(B- ) Carre de la distance au debut de',
299 write(ulsort,90004) ' N',(cono(iaux),iaux=1,sdim)
300 write(ulsort,90004) ' C',
301 > (geocoo(somseg(segmin-1),iaux),iaux=1,sdim), abscur(segmin-1)
302 write(ulsort,90004) ' A',
303 > (geocoo(somseg(segmin),iaux),iaux=1,sdim), abscur(segmin)
304 write(ulsort,90004) ' B',
305 > (geocoo(somseg(segmin+1),iaux),iaux=1,sdim), abscur(segmin+1)
309 if ( daux.lt.daux2 ) then
321 c 4.2. ==> Definition des points A et B
323 #ifdef _DEBUG_HOMARD_
324 if ( glop.ne.0 ) then
325 write(ulsort,90002) '... ==> Le noeud est lie au segment', seg
329 do 42 , iaux = 1 , sdim
330 cooa(iaux) = geocoo(somseg(seg ),iaux)
331 coob(iaux) = geocoo(somseg(seg+1),iaux)
334 #ifdef _DEBUG_HOMARD_
335 if ( glop.ne.0 ) then
336 write(ulsort,90004) ' de A', (cooa(iaux),iaux=1,sdim),
338 write(ulsort,90004) ' a B', (coob(iaux),iaux=1,sdim),
343 c 4.3. ==> Il faut recoller N sur la frontiere le cas echeant,
344 c puis redefinir son abscisse curviligne.
346 c 4.3.1. ==> Calcul du produit scalaire AB.AN
349 do 431 , iaux = 1 , sdim
351 > ( coob(iaux) - cooa(iaux) ) *
352 > ( cono(iaux) - cooa(iaux) )
354 #ifdef _DEBUG_HOMARD_
355 if ( glop.ne.0 ) then
356 write(ulsort,90004) 'Produit scalaire AB.AN', daux
360 c 4.3.2. ==> Positionnement
361 c 4.3.2.1. ==> Si le produit scalaire est negatif, c'est que P est
362 c "en arriere" de A. On ramene P sur A.
364 if ( daux.le.0 ) then
366 do 4321 , iaux = 1 , sdim
367 cono(iaux) = cooa(iaux)
370 #ifdef _DEBUG_HOMARD_
371 if ( glop.ne.0 ) then
372 write(ulsort,90004) '.. ==> Le noeud a ete replace sur A'
378 c daux1 : distance AB
379 daux1 = abscur(seg+1)-abscur(seg)
380 daux = daux / (daux1**2)
381 #ifdef _DEBUG_HOMARD_
382 if ( glop.ne.0 ) then
383 write(ulsort,90004) 'Distance AB', daux1
384 write(ulsort,90004) 'Carre de la distance AB ', daux1**2
388 c 4.3.2.2. ==> Si le produit scalaire est superieur au carre de la
389 c distance AB, c'est que P est "en avant" de B.
390 c On ramenera P sur B.
392 if ( daux.ge.1.d0 ) then
394 do 4322 , iaux = 1 , sdim
395 cono(iaux) = coob(iaux)
397 acnoeu = abscur(seg+1)
398 #ifdef _DEBUG_HOMARD_
399 if ( glop.ne.0 ) then
400 write(ulsort,90004) '.. ==> Le noeud a ete replace sur B'
404 c 4.3.2.3. ==> N est "entre" A et B.
405 c On decompose le vecteur AN en une partie le long du
406 c segment, alpha.AB, et une partie orthogonale.
407 c Cela revient a projeter N orthogonalement au segment.
411 do 4323 , iaux = 1 , sdim
412 cono(iaux) = cooa(iaux) + daux*(coob(iaux)-cooa(iaux))
414 acnoeu = abscur(seg) + daux*daux1
415 #ifdef _DEBUG_HOMARD_
416 if ( glop.ne.0 ) then
417 write(ulsort,90004) '.. ==> Le noeud a ete place entre A et B'
424 #ifdef _DEBUG_HOMARD_
425 if ( glop.ne.0 ) then
426 write(ulsort,90004) '. Nouveau noeud', (cono(iaux),iaux=1,sdim)
434 #ifdef _DEBUG_HOMARD_
435 write(ulsort,*) '5. sortie'
440 #ifdef _DEBUG_HOMARD_
441 if ( glop.ne.0 ) then
442 write(ulsort,*) 'Bilan :'
443 write(ulsort,90002) '==> Segment ', seg
444 write(ulsort,90004) '==> Abcisse curviligne', acnoeu
445 write(ulsort,*) 'Sortie de SFSENO'