src/tool/Utilitaire/utqtet.F

   1       subroutine utqtet ( letetr, qualit, qualij, volume,
   2      >                    coonoe, somare, aretri,
   3      >                    tritet, cotrte, aretet )
   4 c ______________________________________________________________________
   5 c
   6 c                             H O M A R D
   7 c
   8 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
   9 c
  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
  15 c
  16 c    HOMARD est une marque deposee d'Electricite de France
  17 c
  18 c Copyright EDF 1996
  19 c Copyright EDF 1998
  20 c Copyright EDF 2002
  21 c Copyright EDF 2020
  22 c ______________________________________________________________________
  23 c
  24 c     UTilitaire : Qualite d'un TETraedre
  25 c     --           -            ---
  26 c ______________________________________________________________________
  27 c
  28 c    on utilise le critere decrit dans
  29 c     'Maillages, applications aux elements finis'
  30 c     Pascal Jean Frey, Paul-Louis George
  31 c     Hermes, 1999
  32 c     Chapitre 18.2, page 606
  33 c                                           h
  34 c    le critere de qualite, q, vaut alpha * -
  35 c                                           r
  36 c    h est le diametre du tetraedre, i.e. son plus grand cote
  37 c    r est le rayon de la sphere inscrite
  38 c    alpha est un coefficient de normalisation pour que le critere q
  39 c    vaille 1 pour un tetraedre regulier ==> alpha = 1/racine(24)
  40 c
  41 c    pour tout autre tetraedre, le critere est donc superieur a 1
  42 c
  43 c                              max(ak) * somme des si
  44 c    tous calculs faits q vaut ----------------------
  45 c                               3 * racine(24) * vol
  46 c
  47 c    ou si est la surface de la i-eme face,
  48 c       ak est la longueur du k-eme cote
  49 c       vol le volume du tetraedre.
  50 c
  51 c Un tetraedre regulier :
  52 c noeud 1
  53 c     x = 0.5d0
  54 c     y = 0.5d0*sqrt(3.0d0)/3.d0
  55 c     z = sqrt(2.0d0/3.0d0)
  56 c noeud 2
  57 c     x = 0.0d0
  58 c     y = 0.0d0
  59 c     z = 0.0d0
  60 c noeud 3
  61 c     x = 1.0d0
  62 c     y = 0.0d0
  63 c     z = 0.0d0
  64 c noeud 4
  65 c     x = 0.5d0
  66 c     y = 0.5d0*sqrt(3.0d0)
  67 c     z = 0.0d0
  68 c
  69 c     . Jacobien normalise
  70 c ______________________________________________________________________
  71 c .        .     .        .                                            .
  72 c .  nom   . e/s . taille .           description                      .
  73 c .____________________________________________________________________.
  74 c . letetr . e   .  1     . numero du tetraedre a examiner             .
  75 c . qualit .  s  .  1     . qualite                                    .
  76 c . qualij .  s  .  1     . qualite par le jacobien normalise          .
  77 c . volume .  s  .  1     . volume                                     .
  78 c . coonoe . e   . nbnoto . coordonnees des noeuds                     .
  79 c .        .     . * sdim .                                            .
  80 c . somare . e   .2*nbarto. numeros des extremites d'arete             .
  81 c . aretri . e   .nbtrto*3. numeros des 3 aretes des triangles         .
  82 c . tritet . e   .nbtecf*4. numeros des 4 triangles des tetraedres     .
  83 c . cotrte . e   .nbtecf*4. code des 4 triangles des tetraedres        .
  84 c . aretet . e   .nbteca*6. numeros des 6 aretes des tetraedres        .
  85 c .____________________________________________________________________.
  86 c
  87 c====
  88 c 0. declarations et dimensionnement
  89 c====
  90 c
  91 c 0.1. ==> generalites
  92 c
  93       implicit none
  94       save
  95 c
  96 #include "fractl.h"
  97 #include "fracte.h"
  98 c
  99 c 0.2. ==> communs
 100 c
 101 #include "nombno.h"
 102 #include "nombar.h"
 103 #include "nombtr.h"
 104 #include "nombte.h"
 105 c
 106 c 0.3. ==> arguments
 107 c
 108       double precision qualit, qualij, volume, coonoe(nbnoto,3)
 109 c
 110       integer letetr
 111       integer somare(2,nbarto)
 112       integer aretri(nbtrto,3)
 113       integer tritet(nbtecf,4), cotrte(nbtecf,4), aretet(nbteca,6)
 114 c
 115 c 0.4. ==> variables locales
 116 c
 117       integer iaux, jaux
 118       integer aresom(0:3,4)
 119       integer listar(6), listso(4)
 120 c
 121       double precision coosom(3,4)
 122       double precision sf1, sf2, sf3, sf4, sixvol
 123       double precision ar1, ar2, ar3, ar4, ar5, ar6
 124       double precision v1(3), v3(3), v4(3), v6(3), vn(3)
 125       double precision daux
 126 c
 127 c 0.5. ==> initialisations
 128 c ______________________________________________________________________
 129 c
 130 c====
 131 c 1. les aretes et sommets de ce tetraedre
 132 c====
 133 c
 134       call utaste ( letetr,
 135      >              nbtrto, nbtecf, nbteca,
 136      >              somare, aretri,
 137      >              tritet, cotrte, aretet,
 138      >              listar, listso )
 139 c
 140       do 11 , jaux = 1, 4
 141         do 111 , iaux = 1, 3
 142           coosom(iaux,jaux) = coonoe(listso(jaux),iaux)
 143   111   continue
 144    11 continue
 145 c
 146 c====
 147 c 2. longueurs des aretes 2 (n1-n3) et 5 (n2-n4)
 148 c====
 149 c
 150       vn(1) = coosom(1,3) - coosom(1,1)
 151       vn(2) = coosom(2,3) - coosom(2,1)
 152       vn(3) = coosom(3,3) - coosom(3,1)
 153       ar2 = sqrt ( vn(1)*vn(1) + vn(2)*vn(2) + vn(3)*vn(3) )
 154 c
 155       vn(1) = coosom(1,4) - coosom(1,2)
 156       vn(2) = coosom(2,4) - coosom(2,2)
 157       vn(3) = coosom(3,4) - coosom(3,2)
 158       ar5 = sqrt ( vn(1)*vn(1) + vn(2)*vn(2) + vn(3)*vn(3) )
 159 c
 160 c====
 161 c 3. memorisation des vecteurs liees aux aretes 1 (n1-n2),
 162 c          3 (n1-n4), 4 (n2-n3) et 6 (n3-n4)
 163 c        et calcul des longueurs de ces aretes
 164 c====
 165 c
 166       v1(1) = coosom(1,2) - coosom(1,1)
 167       v1(2) = coosom(2,2) - coosom(2,1)
 168       v1(3) = coosom(3,2) - coosom(3,1)
 169       ar1 = sqrt ( v1(1)*v1(1) + v1(2)*v1(2) + v1(3)*v1(3) )
 170 c
 171       v3(1) = coosom(1,4) - coosom(1,1)
 172       v3(2) = coosom(2,4) - coosom(2,1)
 173       v3(3) = coosom(3,4) - coosom(3,1)
 174       ar3 = sqrt ( v3(1)*v3(1) + v3(2)*v3(2) + v3(3)*v3(3) )
 175 c
 176       v4(1) = coosom(1,3) - coosom(1,2)
 177       v4(2) = coosom(2,3) - coosom(2,2)
 178       v4(3) = coosom(3,3) - coosom(3,2)
 179       ar4 = sqrt ( v4(1)*v4(1) + v4(2)*v4(2) + v4(3)*v4(3) )
 180 c
 181       v6(1) = coosom(1,4) - coosom(1,3)
 182       v6(2) = coosom(2,4) - coosom(2,3)
 183       v6(3) = coosom(3,4) - coosom(3,3)
 184       ar6 = sqrt ( v6(1)*v6(1) + v6(2)*v6(2) + v6(3)*v6(3) )
 185 cgn      write(1,*) ar1,ar2,ar3,ar4,ar5,ar6
 186 c
 187 c====
 188 c 4. calcul des 4 surfaces (plutot 2 fois les surfaces)
 189 c           on rappelle que la surface d'un triangle est egale
 190 c           a la moitie de la norme du produit vectoriel de deux
 191 c           des vecteurs representant les aretes.
 192 c====
 193 c
 194       vn(1) = v6(2)*v4(3) - v6(3)*v4(2)
 195       vn(2) = v6(3)*v4(1) - v6(1)*v4(3)
 196       vn(3) = v6(1)*v4(2) - v6(2)*v4(1)
 197       sf1 = vn(1)*vn(1) + vn(2)*vn(2) + vn(3)*vn(3)
 198 c
 199       vn(1) = v6(2)*v3(3) - v6(3)*v3(2)
 200       vn(2) = v6(3)*v3(1) - v6(1)*v3(3)
 201       vn(3) = v6(1)*v3(2) - v6(2)*v3(1)
 202       sf2 = vn(1)*vn(1) + vn(2)*vn(2) + vn(3)*vn(3)
 203 c
 204       vn(1) = v1(2)*v3(3) - v1(3)*v3(2)
 205       vn(2) = v1(3)*v3(1) - v1(1)*v3(3)
 206       vn(3) = v1(1)*v3(2) - v1(2)*v3(1)
 207       sf3 = vn(1)*vn(1) + vn(2)*vn(2) + vn(3)*vn(3)
 208 c
 209       vn(1) = v1(2)*v4(3) - v1(3)*v4(2)
 210       vn(2) = v1(3)*v4(1) - v1(1)*v4(3)
 211       vn(3) = v1(1)*v4(2) - v1(2)*v4(1)
 212       sf4 = vn(1)*vn(1) + vn(2)*vn(2) + vn(3)*vn(3)
 213 cgn      write(1,*) sf1,sf2,sf3,sf4
 214 c
 215 c====
 216 c 5. calcul du volume du tetraedre
 217 c====
 218 c
 219       call utvot0 ( coosom,  volume )
 220 c
 221 c====
 222 c 6. volume et qualite
 223 c====
 224 c
 225       sixvol = 6.d0*volume
 226 c
 227       qualit = sqrt(uns24) * max(ar1,ar2,ar3,ar4,ar5,ar6)
 228      >     * (sqrt(sf1)+sqrt(sf2)+sqrt(sf3)+sqrt(sf4)) / sixvol
 229 c
 230 c====
 231 c 7. qualite par le jacobien normalise
 232 c====
 233 c 7.1. ==> Liens sommet/aretes
 234 c
 235       aresom(0,1) = 1
 236       aresom(1,1) = 1
 237       aresom(2,1) = 3
 238       aresom(3,1) = 2
 239 c
 240       aresom(0,2) = 2
 241       aresom(1,2) = 5
 242       aresom(2,2) = 1
 243       aresom(3,2) = 4
 244 c
 245       aresom(0,3) = 3
 246       aresom(1,3) = 4
 247       aresom(2,3) = 2
 248       aresom(3,3) = 6
 249 c
 250       aresom(0,4) = 4
 251       aresom(1,4) = 6
 252       aresom(2,4) = 3
 253       aresom(3,4) = 5
 254 c
 255 c 7.2. ==> fonction generique
 256 c
 257       iaux = 4
 258       daux = sqrt(2.d0)/2.d0
 259       call utqjno (   iaux, aresom,   daux,
 260      >              listar, listso, somare, coonoe,
 261      >              qualij )
 262 cgn      write(1,*) '==> qualij : ', qualij
 263 c
 264       end