====================================================================

  programmes Maple  -  24 mars 2006 - 2 décembre 2009

  Méthode des éléments finis : élasticité plane

     Yves DEBARD

     Institut Universitaire de Technologie du MANS
     Département Génie Mécanique et Productique
     Avenue Olivier Messiaen
     72085 LE MANS CEDEX 9

====================================================================

  Eléments iso-paramétriques

  Calcul des fonctions d'interpolation et de leurs dérivées


tri3_int  : triangle à 3 noeuds

restart:with(linalg):

n:=3;

# base polynomiale

P:=(xi,eta)->[1,xi,eta];

# coordonnées nodales

xi_nod:=[0,1,0]:
eta_nod:=[0,0,1]:

C:=matrix([seq(P(xi_nod[i],eta_nod[i]),i=1..n)]);

# fonctions d'interpolation

N:=multiply(P(xi,eta),inverse(C));

dNxi:=map(diff,N,xi);
dNeta:=map(diff,N,eta);

-----------------------------------

tri6_int : triangle à 6 noeuds

restart:with(linalg):

n:=6;

# base polynomiale

P:=(xi,eta)->[1,xi,eta,xi*xi,eta*eta,xi*eta];

# coordonnées nodales

xi_nod:=[0,1/2,1,1/2,0,0]:
eta_nod:=[0,0,0,1/2,1,1/2]:

C:=matrix([seq(P(xi_nod[i],eta_nod[i]),i=1..n)]);

# fonctions d'interpolation

N:=multiply(P(xi,eta),inverse(C));

N:=map(factor,N);

dNxi:=map(diff,N,xi);
dNeta:=map(diff,N,eta);

------------------------------------

quad4_int : quadrangle à 4 noeuds

restart:with(linalg):

n:=4;

# base polynomiale

P:=(xi,eta)->[1,xi,eta,xi*eta];

# coordonnées nodales

xi_nod:=[-1,1,1,-1]:
eta_nod:=[-1,-1,1,1]:

C:=matrix([seq(P(xi_nod[i],eta_nod[i]),i=1..n)]);

# fonctions d'interpolation

N:=multiply(P(xi,eta),inverse(C));

N:=map(factor,N);

dNxi:=map(diff,N,xi);
dNeta:=map(diff,N,eta);

# représentations graphiques

for i from 1 to n do
  plot3d(N[i],xi=-1..1,eta=-1..1,axes=boxed);
end do;

------------------------------------

quad8_int : quadrangle à 8 noeuds

restart:with(linalg):

n:=8;

# base polynomiale

P:=(xi,eta)->[1,xi,eta,xi*eta,xi^2,eta^2,xi^2*eta,xi*eta^2];

# coordonnées nodales

xi_nod:=[-1,0,1,1,1,0,-1,-1]:
eta_nod:=[-1,-1,-1,0,1,1,1,0]:

C:=matrix([seq(P(xi_nod[i],eta_nod[i]),i=1..n)]);

# fonctions d'interpolation

N:=multiply(P(xi,eta),inverse(C));

N:=map(factor,N);

dNxi:=map(diff,N,xi);
dNeta:=map(diff,N,eta);

# représentations graphiques

for i from 1 to n do
  plot3d(N[i],xi=-1..1,eta=-1..1,axes=boxed);
end do;

----------------------------------------

quad9_int : quadrangle à 9 noeuds

restart:with(linalg):

n:=9;

# base polynomiale

P:=(xi,eta)->[1,xi,eta,xi*eta,xi^2,eta^2,xi^2*eta,xi*eta^2,xi^2*eta^2];

# coordonnées nodales

xi_nod:=[-1,0,1,1,1,0,-1,-1,0]:
eta_nod:=[-1,-1,-1,0,1,1,1,0,0]:

C:=matrix([seq(P(xi_nod[i],eta_nod[i]),i=1..n)]);

# fonctions d'interpolation

N:=multiply(P(xi,eta),inverse(C));

N:=map(factor,N);

dNxi:=map(diff,N,xi);
dNeta:=map(diff,N,eta);

# représentations graphiques

for i from 1 to n do
  plot3d(N[i],xi=-1..1,eta=-1..1,axes=boxed);
end do;


------------------------------------

qualité du jacobien : quadrangle à 4 noeuds

restart:with(linalg):

n:=4;

# base polynomiale

P:=(xi,eta)->[1,xi,eta,xi*eta];

# coordonnées nodales

xi_nod:=[-1,1,1,-1]:
eta_nod:=[-1,-1,1,1]:

C:=matrix([seq(P(xi_nod[i],eta_nod[i]),i=1..n)]);

# fonctions d'interpolation

N:=multiply(P(xi,eta),inverse(C));

N:=map(factor,N);

dNxi:=map(diff,N,xi);
dNeta:=map(diff,N,eta);

dN:=matrix([dNxi,dNeta]);


coord:=matrix([[d,d],[L,0],[L,L],[0,L]]);

J:=multiply(dN,coord):J:=simplify(%);

detJ:=det(J);

-------------------------------------

qualité du jacobien : triangle à 6 noeuds

restart:with(linalg):

n:=6;

# base polynomiale

P:=(xi,eta)->[1,xi,eta,xi*xi,eta*eta,xi*eta];

# coordonnées nodales

xi_nod:=[0,1/2,1,1/2,0,0]:
eta_nod:=[0,0,0,1/2,1,1/2]:

C:=matrix([seq(P(xi_nod[i],eta_nod[i]),i=1..n)]);

# fonctions d'interpolation

N:=multiply(P(xi,eta),inverse(C));

N:=map(factor,N);

dNxi:=map(diff,N,xi);
dNeta:=map(diff,N,eta);

dN:=matrix([dNxi,dNeta]);


coord:=matrix([[0,0],[L/2,h],[L,0],[L/2,L/2],[0,L],[0,L/2]]);

J:=multiply(dN,coord):J:=simplify(%);

detJ:=det(J);