titr 'Thermique' ;
****************** thermique3.dgibi
opti dime 2 elem tri3 ;
saut ligne ;
******** paramètres
l = 10.0 ;
r = 3.0 ;
n = 10 ;
n2 = enti (n/2.0) ;
n4 = enti (n/3.0) ;
mf = 1.414*n ;
m = enti mf ;
******** géométrie extérieure
p1 = 0.0 0.0 ;
p2 = l 0.0 ;
p3 = l (2.0*l) ;
p4 = (2.0*l) (2.0*l) ;
p5 = (2.0*l) (3.0*l) ;
p6 = 0.0 (3.0*l) ;
p7 = 0.0 (2.0*l) ;
p8 = ((-1)*0.50*l) (1.50*l) ;
p9 = 0.0 l ;
c1 = 0.0 (1.50*l) ;
l1 = droi n p1 p2 ;
l2 = droi (2*n) p2 p3 ;
l3 = droi n p3 p4 ;
l4 = droi n p4 p5 ;
l5 = droi (2*n) p5 p6 ;
l6 = droi n p6 p7 ;
l7 = cerc n2 p7 c1 p8 ;
l8 = cerc n2 p8 c1 p9 ;
l9 = droi n p9 p1 ;
li1 = droi m p3 p6 ;
cont1 = l1 et l2 et li1 et l6 et l7 et l8 et l9 ;
cont2 = l3 et l4 et l5 et li1 ;
********* géométrie intérieure
c2 = (0.50*l) (1.50*l) ;
pi1 = (0.50*l) (1.50*l+r) ;
pi2 = (0.50*l-r) (1.50*l) ;
pi3 = (0.50*l) (1.50*l-r) ;
pi4 = (0.50*l+r) (1.50*l) ;
li1 = cerc n4 pi1 c2 pi2 ; 
li2 = cerc n4 pi2 c2 pi3 ; 
li3 = cerc n4 pi3 c2 pi4 ; 
li4 = cerc n4 pi4 c2 pi1 ; 
conti1 = li1 et li2 et li3 et li4 ;
conti1 = inve conti1 ;
c2 = (1.50*l) (2.50*l) ;
r2 = r / 2.0 ;
pi1 = (1.50*l) (2.50*l+r2) ;
pi2 = (1.50*l-r2) (2.50*l) ;
pi3 = (1.50*l) (2.50*l-r2) ;
pi4 = (1.50*l+r2) (2.50*l) ;
li1 = cerc n4 pi1 c2 pi2 ; 
li2 = cerc n4 pi2 c2 pi3 ; 
li3 = cerc n4 pi3 c2 pi4 ; 
li4 = cerc n4 pi4 c2 pi1 ; 
conti2 = li1 et li2 et li3 et li4 ;
conti2 = inve conti2 ;
trac qual (cont1 et cont2 et conti1 et conti2) ;
******** maillage
mail1 = surf (cont1 et conti1) plan ;
mail2 = surf (cont2 et conti2) plan ;
mail0 = mail1 et mail2 ;
*list mail0 ;
trac qual mail0 ;
******** modèle
k1 = 50.0 ;
mod1 = mode mail1 thermique isotrope ;
mat1 = mate mod1 k k1 ;
k2 = 50.0 ;
mod2 = mode mail2 thermique isotrope ;
mat2 = mate mod2 k k2 ;
******** matrice de rigidité
cond1 = cond mod1 mat1 ;
cond2 = cond mod2 mat2 ;
condk = cond1 et cond2 ;
******** conditions limites
cl1 = bloq t l1 ;
dcl1 = depi cl1 100.0 ;
******** chargement
mod3 = mode l4 convection ;
mat3 = mate mod3 h 20.0 ;
flu3 = conv mod3 mat3 t 10.0 ;
******** modif' de la matrice de conductivité 
cond3 = cond mod3 mat3 ;
******** matrice de rigidité modifiée par conditions limites
condk0 = condk et cond3 et cl1 ;
******** second membre
flut0 = flu3 et dcl1 ;
******** résolution
temp0 = reso condk0 flut0 ;
*list temp0 ;
******** post traitement
titre 'Champ de temperatures' ;
trac temp0 mail0 ;
******** grad T sur les éléments des 2 domaines
grad1 = grad mod1 temp0 ;
grad2 = grad mod2 temp0 ;
******** vecteur flux sur les éléments des 2 domaines
flux1 = (-1)*k1*grad1 ;
flux2 = (-1)*k2*grad2 ;
flux0 = flux1 et flux2 ;
*list flux0 ;
mod0 = mod1 et mod2 ;
******** vecteur flux (q) moyenné aux noeuds
fluxno0 = chan chpo flux0 mod0 ;
amp0 = 0.010 ;
vflux0 = vect fluxno0 amp0 t,x t,y rouge ;
titre 'Vecteur flux' ;
trac mail0 vflux0 ;
******** module du flux
fluxx0 = exco t,x fluxno0 ;
fluxy0 = exco t,y fluxno0 ;
mflux0 = ( fluxx0**2 + fluxy0**2 )**0.50 ;
titre 'Module du vecteur flux' ;
trac mail0 mflux0 ;
******** q.-y sur l1
fluxyl1 = evol rouge chpo fluxy0 l1 ;
*list fluxyl1 ;
dess fluxyl1 titx 'ligne 1' tity 'q.n (W/m2)' 
             titr 'flux suivant y sur l1' ;
puis1 = somm  fluxyl1 ;
mess "Puissance traversant ligne 1 (W) :" ;
list puis1 ;
******** - q.x sur l4
fluxxl4 = evol rouge chpo fluxx0 l4 ;
fluxxl4 = fluxxl4*(-1) ;
dess fluxxl4 titx 'ligne 4' tity 'q.n (W/m2)' 
             titr 'flux suivant x sur l4' ;
puis4 = somm  fluxxl4 ;
mess "Puissance traversant ligne 4 (W) :" ;
list puis4 ;
********
fin ;