1 subroutine vcms22 ( maextr,
3 > nbse2d, nbtr2d, nbqu2d, nbele2,
5 > fameel, typele, noeele,
6 > fame2d, type2d, noee2d,
8 > ulsort, langue, codret )
9 c ______________________________________________________________________
13 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
15 c Version originale enregistree le 18 juin 1996 sous le numero 96036
16 c aupres des huissiers de justice Simart et Lavoir a Clamart
17 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
18 c aupres des huissiers de justice
19 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
21 c HOMARD est une marque deposee d'Electricite de France
27 c ______________________________________________________________________
29 c aVant adaptation - Conversion de Maillage -
31 c Saturne 2D - phase 2 - Neptune 2D
33 c ______________________________________________________________________
35 c . nom . e/s . taille . description .
36 c .____________________________________________________________________.
37 c . maextr . e . 1 . maillage extrude .
38 c . . . . 0 : non (defaut) .
39 c . . . . 1 : selon X .
40 c . . . . 2 : selon Y .
41 c . . . . 3 : selon Z (cas de Saturne ou Neptune) .
42 c . nbnoto . e . 1 . nombre de noeuds du maillage externe .
43 c . nbse2d . e . 1 . nombre de segments du maillage 2d .
44 c . nbtr2d . e . 1 . nombre de triangles du maillage 2d .
45 c . nbqu2d . e . 1 . nombre de quadrangles du maillage 2d .
46 c . nbelem . e . 1 . nombre d'elements du maillage externe .
47 c . nu2dno . e . nbnoto . numero du calcul des noeuds saturne/neptune.
48 c . coonoe . e . nbnoto . coordonnees des noeuds .
50 c . fameel . e . nbelem . famille med des elements .
51 c . typele . e . nbelem . type des elements pour le code de calcul .
52 c . noeele . e . nbelem . noeuds des elements .
54 c . fame2d . s . nbele2 . famille med des elements du maillage 2d .
55 c . type2d . s . nbele2 . type des elements du maillage 2d .
56 c . noee2d . s . nbele2 . noeuds des elements du maillage 2d .
58 c . faminf . s . 1 . famille med des quad de la face inferieure .
59 c . famsup . s . 1 . famille med des quad de la face superieure .
60 c . ulsort . e . 1 . numero d'unite logique de la liste standard.
61 c . langue . e . 1 . langue des messages .
62 c . . . . 1 : francais, 2 : anglais .
63 c . codret . es . 1 . code de retour des modules .
64 c . . . . 0 : pas de probleme .
65 c ______________________________________________________________________
68 c 0. declarations et dimensionnement
71 c 0.1. ==> generalites
77 parameter ( nompro = 'VCMS22' )
97 integer nbse2d, nbtr2d, nbqu2d, nbele2
99 integer faminf, famsup
100 integer nu2dno(nbnoto)
101 integer fameel(nbelem), typele(nbelem), noeele(nbelem,*)
102 integer fame2d(nbele2), type2d(nbele2), noee2d(nbele2,*)
104 double precision coonoe(nbnoto,*)
106 integer ulsort, langue, codret
108 c 0.4. ==> variables locales
110 integer iaux, jaux, kaux, laux
118 double precision xymil(2)
119 double precision v1(2), vn(2)
120 double precision daux1, daux2
121 double precision epsilo
122 double precision daux(0:4)
125 parameter ( nbmess = 20 )
126 character*80 texte(nblang,nbmess)
128 c 0.5. ==> initialisations
129 c ______________________________________________________________________
137 #ifdef _DEBUG_HOMARD_
138 write (ulsort,texte(langue,1)) 'Entree', nompro
142 texte(1,4) = '(''Maille numero :'',i10,'', de noeuds '',8i10)'
143 texte(1,5) = '(i1,'' noeud(s) sont dans le plan cooinf.'')'
144 texte(1,6) = '(''Pour un '',a,'', il en faudrait '',a)'
145 texte(1,7) = '(''Famille de la face '',a,'' : '',i6)'
146 texte(1,8) = '(''Famille du '',a,i10,'' : '',i6)'
148 >'(''Nombre de '',a,'' attendus pour le maillage 2D :'',i10)'
150 >'(''Nombre de '',a,'' trouves pour le maillage 2D :'',i10)'
151 texte(1,11) = '(''Element '',i10,'' ('',a,''), numloc = '',4i10)'
153 texte(2,4) = '(''Mesh # :'',i10,'', with nodes '',8i10)'
154 texte(2,5) = '(i1,'' node(s) are in cooinf plane.'')'
155 texte(2,6) = '(''For '',a,'', '',a,'' were expected.'')'
156 texte(2,7) = '(''Family for '',a,'' face : '',i6)'
157 texte(2,8) = '(''Family for '',a,'' #'',i10,'' : '',i6)'
159 > '(''Expected number of '',a,'' for the 2D mesh :'',i10)'
161 > '(''Found number of '',a,'' for the 2D mesh :'',i10)'
162 texte(2,11) = '(''Element '',i10,'' ('',a,''), numloc = '',4i10)'
168 if ( maextr.eq.1 ) then
172 elseif ( maextr.eq.2 ) then
176 elseif ( maextr.eq.3 ) then
185 #ifdef _DEBUG_HOMARD_
186 write (ulsort,90002) 'maextr', maextr
187 write (ulsort,90002) 'nbelem', nbelem
188 write (ulsort,90002) 'nbele2', nbele2
189 write (ulsort,90002) 'nbse2d', nbse2d
190 write (ulsort,90002) 'nbqu2d', nbqu2d
191 write (ulsort,90002) 'nbtr2d', nbtr2d
197 c 2. transformations des quadrangles en segments
198 c Attention a les mettre avant les quadrangles pour respecter
199 c les conventions ...
201 #ifdef _DEBUG_HOMARD_
202 write (ulsort,90002) '2. quad -> segm ; codret', codret
208 do 21 , el = 1 , nbelem
210 if ( codret.eq.0 ) then
212 if ( typele(el).eq.edqua4 ) then
214 c 2.1. ==> recherche des noeuds dans le plan cooinf
217 do 211 , iaux = 1 , 4
218 if ( nu2dno(noeele(el,iaux)).ne.0 ) then
220 numloc(jaux) = noeele(el,iaux)
221 if ( jaux.eq.1 ) then
229 #ifdef _DEBUG_HOMARD_
230 if ( jaux.gt.0 .and. el.lt.1 ) then
231 write (ulsort,texte(langue,11)) el, 'quad',
232 > (numloc(iaux),iaux=1,jaux)
236 c 2.2. ==> si exactement 2 noeuds sont dans le plan, on cree le segment
237 c Le segment est cree avec les noeuds dans l'ordre
238 c d'apparition dans la connectivite du quadrangle. Cela
239 c permettra a la reconstitution du quadrangle apres adaptation
240 c de retrouver la meme orientation de la face.
241 c attention au retournemnt eventuel ...
243 if ( jaux.eq.2 ) then
246 if ( kaux.eq.1 .and. laux.eq.4 ) then
248 numloc(2) = numloc(1)
251 noee2d(nuel2d,1) = nu2dno(numloc(1))
252 noee2d(nuel2d,2) = nu2dno(numloc(2))
253 fame2d(nuel2d) = fameel(el)
254 type2d(nuel2d) = edseg2
256 c 2.3. ==> si exactement 4 noeuds sont dans ce plan, c'est la face
257 c inferieure. On memorise le numero de famille
259 elseif ( jaux.eq.4 ) then
261 if ( faminf.eq.1 ) then
264 if ( fameel(el).ne.faminf ) then
265 write (ulsort,texte(langue,7)) 'inf', faminf
266 write (ulsort,texte(langue,8))
267 > mess14(langue,1,4), el, fameel(el)
272 c 2.4. ==> si exactement 4 noeuds sont dans l'autre plan, c'est la face
273 c superieure. On memorise le numero de famille
275 elseif ( jaux.eq.0 ) then
277 if ( famsup.eq.1 ) then
280 if ( fameel(el).ne.famsup ) then
281 write (ulsort,texte(langue,7)) 'sup', famsup
282 write (ulsort,texte(langue,8))
283 > mess14(langue,1,4), el, fameel(el)
288 c 2.5. ==> sinon, c'est louche ...
292 write (ulsort,texte(langue,4)) el,(noeele(el,iaux),iaux=1,4)
293 write (ulsort,texte(langue,5)) jaux
294 write (ulsort,texte(langue,6)) mess14(langue,1,4),
305 if ( codret.eq.0 ) then
307 cgn write (ulsort,90002) 'nuel2d', nuel2d
308 cgn write (ulsort,90002) 'nbse2d', nbse2d
309 if ( nuel2d.ne.nbse2d ) then
310 write (ulsort,texte(langue,9)) mess14(langue,3,1), nbse2d
311 write (ulsort,texte(langue,10))
312 > mess14(langue,3,1), nuel2d-nbse2d
318 #ifdef _DEBUG_HOMARD_
319 write (ulsort,texte(langue,7)) 'inf', faminf
320 write (ulsort,texte(langue,7)) 'sup', famsup
324 c 3. transformations des hexaedres en quadrangles
326 #ifdef _DEBUG_HOMARD_
327 write (ulsort,90002) '3. hexa -> quad ; codret', codret
330 if ( nbqu2d.ne.0 ) then
332 do 31 , el = 1 , nbelem
334 if ( codret.eq.0 ) then
336 if ( typele(el).eq.edhex8 ) then
338 c 3.1. ==> recherche des noeuds dans le plan cooinf
340 cgn write (ulsort,90012) 'noeele',el,(noeele(el,iaux),iaux=1,8)
342 do 311 , iaux = 1 , 8
343 if ( nu2dno(noeele(el,iaux)).ne.0 ) then
345 numloc(jaux) = noeele(el,iaux)
349 #ifdef _DEBUG_HOMARD_
350 write (ulsort,90002) 'jaux', jaux
351 if ( nuel2d.eq.nbse2d.or. nuel2d.eq.988.or.nuel2d.eq.987 ) then
352 write (ulsort,texte(langue,11)) el, 'hexa',
353 > (numloc(iaux),iaux=1,jaux)
354 write (ulsort,90004) 'X',(coonoe(numloc(iaux),1),iaux = 1,jaux)
355 write (ulsort,90004) 'Y',(coonoe(numloc(iaux),2),iaux = 1,jaux)
356 write (ulsort,90004) 'Z',(coonoe(numloc(iaux),3),iaux = 1,jaux)
360 c 3.2. ==> si exactement 4 noeuds sont dans ce plan, creation du
361 c quadrangle asssocie
362 c attention a bien ranger les noeuds pour ne pas croiser !
363 c pour cela, on les positionne en fonction de leur milieu :
372 c ........M...........> X
379 if ( jaux.eq.4 ) then
381 c Le milieu du quadrangle
385 do 313 , iaux = 1 , 4
386 xymil(1) = xymil(1) + coonoe(numloc(iaux),iaux1)
387 xymil(2) = xymil(2) + coonoe(numloc(iaux),iaux2)
389 xymil(1) = unsqu * xymil(1)
390 xymil(2) = unsqu * xymil(2)
391 cgn write (ulsort,90004) 'xymil', xymil
393 c Le vecteur entre le milieu et le premier noeud
395 v1(1) = coonoe(numloc(1),iaux1) - xymil(1)
396 v1(2) = coonoe(numloc(1),iaux2) - xymil(2)
397 daux1 = sqrt(v1(1)**2+v1(2)**2)
400 cgn write (ulsort,90004) 'v1', v1
402 c Les angles entre le segment (milieu,noeud 1) et
403 c (milieu,noeud suivant)
405 do 314 , iaux = 2 , 4
406 vn(1) = coonoe(numloc(iaux),iaux1) - xymil(1)
407 vn(2) = coonoe(numloc(iaux),iaux2) - xymil(2)
408 daux1 = sqrt(vn(1)**2+vn(2)**2)
411 cgn write (ulsort,*) ' '
412 cgn write (ulsort,90114) 'vn', iaux,vn
413 daux1 = v1(1)*vn(1) + v1(2)*vn(2)
414 daux2 = v1(1)*vn(2) - v1(2)*vn(1)
415 cgn write (ulsort,90114) 'p scal', iaux, daux1
416 cgn write (ulsort,90114) 'p vect', iaux, daux2
418 if ( (daux1+1.d0).le.zeroma ) then
421 daux(iaux) = acos(daux1)
423 if ( daux2.le.0.d0 ) then
424 daux(iaux) = deuxpi - daux(iaux)
427 #ifdef _DEBUG_HOMARD_
428 if ( nuel2d.lt.1 ) then
429 write (ulsort,90004) 'angles', daux(2), daux(3), daux(4)
433 c Classement des angles
435 daux1 = min(daux(2), daux(3), daux(4))
436 daux2 = max(daux(2), daux(3), daux(4))
437 cgn write (ulsort,90004) 'angles min', daux1
438 cgn write (ulsort,90004) 'angles max', daux2
440 do 315 , iaux = 2 , 4
441 if ( abs(daux(iaux)-daux1).le.epsilo ) then
444 if ( abs(daux(iaux)-daux2).le.epsilo ) then
449 tbiaux(3) = fp0123(tbiaux(2)-1, tbiaux(4)-1) + 1
450 cgn write (ulsort,90002) 'tbiaux final ', tbiaux
452 c Transfert de la connectivite
455 do 316 , iaux = 1 , 4
456 noee2d(nuel2d,iaux) = nu2dno(numloc(tbiaux(iaux)))
458 fame2d(nuel2d) = fameel(el)
459 type2d(nuel2d) = edqua4
460 #ifdef _DEBUG_HOMARD_
461 if ( nuel2d.lt.1 ) then
462 write (ulsort,90015) 'typele(',nuel2d,') = ',type2d(nuel2d)
463 do 32221 , iaux = 1 , 4
465 > ' noeele',nuel2d,iaux,noee2d(nuel2d,iaux)
471 write (ulsort,texte(langue,4)) el,(noeele(el,iaux),iaux=1,8)
472 write (ulsort,texte(langue,5)) jaux
473 write (ulsort,texte(langue,6)) mess14(langue,1,9), '4'
484 if ( codret.eq.0 ) then
486 cgn write (ulsort,90002) 'nuel2d', nuel2d
487 cgn write (ulsort,90002) 'nbqu2d', nbqu2d
488 cgn write (ulsort,90002) 'nuel2d-nbse2d', nuel2d-nbse2d
489 if ( (nuel2d-nbse2d).ne.nbqu2d ) then
490 write (ulsort,texte(langue,9)) mess14(langue,3,4), nbqu2d
491 write (ulsort,texte(langue,10)) mess14(langue,3,4),
501 c 4. transformations des pentaedres en triangles
503 #ifdef _DEBUG_HOMARD_
504 write (ulsort,90002) '4. pent -> tria ; codret', codret
507 if ( nbtr2d.ne.0 ) then
509 if ( codret.eq.0 ) then
511 do 41 , el = 1 , nbelem
513 if ( typele(el).eq.edpen6 ) then
515 c 4.1. ==> recherche des noeuds dans le plan cooinf
518 do 411 , iaux = 1 , 6
519 if ( nu2dno(noeele(el,iaux)).ne.0 ) then
521 numloc(jaux) = noeele(el,iaux)
525 #ifdef _DEBUG_HOMARD_
526 if ( jaux.gt.0 ) then
527 write (ulsort,texte(langue,11)) el, 'pent',
528 > (numloc(iaux),iaux=1,jaux)
532 c 4.2. ==> si exactement 3 noeuds sont dans ce plan, creation du
535 if ( jaux.eq.3 ) then
538 noee2d(nuel2d,1) = nu2dno(numloc(3))
539 noee2d(nuel2d,2) = nu2dno(numloc(2))
540 noee2d(nuel2d,3) = nu2dno(numloc(1))
541 fame2d(nuel2d) = fameel(el)
542 type2d(nuel2d) = edtri3
545 write (ulsort,texte(langue,4)) el,(noeele(el,iaux),iaux=1,8)
546 write (ulsort,texte(langue,5)) jaux
547 write (ulsort,texte(langue,6)) mess14(langue,1,9), '4'
556 if ( codret.eq.0 ) then
558 cgn write (ulsort,90002) 'nuel2d', nuel2d
559 cgn write (ulsort,90002) 'nbtr2d', nbtr2d
560 cgn write (ulsort,90002) 'nuel2d-nbse2d-nbqu2d', nuel2d-nbse2d-nbqu2d
561 if ( (nuel2d-nbse2d-nbqu2d).ne.nbtr2d ) then
562 write (ulsort,texte(langue,9)) mess14(langue,3,2), nbtr2d
563 write (ulsort,texte(langue,10)) mess14(langue,3,2),
564 > nuel2d-nbse2d-nbqu2d
578 #ifdef _DEBUG_HOMARD_
579 write (ulsort,90002) '5. fin ; codret', codret
582 if ( codret.ne.0 ) then
586 write (ulsort,texte(langue,1)) 'Sortie', nompro
587 write (ulsort,texte(langue,2)) codret
590 #ifdef _DEBUG_HOMARD_
591 write (ulsort,texte(langue,1)) 'Sortie', nompro