x y
dx d d
x y
d X = dy = d = J d = J d
T T
T
d X = (d JJ T d )1/2
clc
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% basics on how to compute an integral over a curved line
% compute the vector dX= (dx,dy)’ given dKSI= (dksi,deta)’
% and its norm dL= ||dX||
% example of Q9 element
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms ksi eta dksi deta dKSI dETA real;
syms x1 x2 x3 x4 x5 x6 x7 x8 x9 real;
syms y1 y2 y3 y4 y5 y6 y7 y8 y9 real ;
coor(1,1)= x1;
coor(1,2)= y1;
coor(2,:)= [x2,y2];
coor(3,:)= [x3,y3];
coor(4,:)= [x4,y4];
coor(5,:)= [x5,y5];
coor(6,:)= [x6,y6];
coor(7,:)= [x7,y7];
coor(8,:)= [x8,y8];
coor(9,:)= [x9,y9];
coor;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[N ,B] = Q9Shapefunctions(ksi,eta );
% this function can be found with COURS6
% Recall that x(ksi,eta) = Sum (Ni xi) and y(ksi,eta)= sum(Ni yi)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
J= B*coor;
%choose the edge number (1,2,3 or 4)
edge=1;
switch (edge)
case 1
eta= -1;
eval(J);
deta=0;
dKSI(1,1)= dksi;
dKSI(2,1)= deta;
dX= J'*dKSI;
'Vector dX'
eval(dX)
'norm'
norm( dX)
case 2
ksi= 1;
eval(J);
dksi=0;
dKSI(1,1)= dksi;
dKSI(2,1)= deta;
dX= J'*dKSI;
'Vector dX'
eval(dX)
'norm'
norm( dX)
case 3
eta= 1;
eval(J);
deta=0;
dKSI(1,1)= dksi;
dKSI(2,1)= deta;
dX= J'*dKSI;
'Vector dX'
eval(dX)
'norm'
norm( dX)
case 4
ksi= -1;
eval(J);
dksi=0;
dKSI(1,1)= dksi;
dKSI(2,1)= deta;
dX= J'*dKSI;
'Vector dX'
eval(dX)
'norm'
norm( dX)
end
RESULTS:
Edge=1
dx= dksi*(x1*(ksi - 1/2) - 2*ksi*x2 + x3*(ksi + 1/2))
dy=dksi*(y1*(ksi - 1/2) - 2*ksi*y2 + y3*(ksi + 1/2))
Edge=2
dx=deta*(x3*(eta - 1/2) - 2*eta*x4 + x5*(eta + 1/2))
dy=deta*(y3*(eta - 1/2) - 2*eta*y4 + y5*(eta + 1/2))
Edge= 3
dx=dksi*(x5*(ksi + 1/2) - 2*ksi*x6 + x7*(ksi - 1/2))
dy=dksi*(y5*(ksi + 1/2) - 2*ksi*y6 + y7*(ksi - 1/2))
pose s= -ksi -1<= s <=1 (on this edge: node 5 is the first node and node 7 is the last one!)
dx=-ds*(x5*(-s + 1/2) + 2*s*x6 + x7*(-s - 1/2)) = ds (x5*(s - 1/2) -2*s*x6 + x7*(s + 1/2))
dy= -ds*(y5*(-s+ 1/2) +2*ksi*y6 + y7*(-s- 1/2))= ds (y5*(s- 1/2) -2*ksi*y6 + y7*(s+1/2)
Edge= 4
dx=eta*(x1*(eta - 1/2) - 2*eta*x8 + x7*(eta + 1/2))
dy=deta*(y1*(eta - 1/2) - 2*eta*y8 + y7*(eta + 1/2))
pose s= -eta -1<= s <=1 (on this edge: node 7 is the first node and node 1 is the last one!)
dx=-ds*(x1*(-s - 1/2) + 2*s*x8 + x7*(-s + 1/2))= ds*( x7*(s - 1/2)- 2*s*x8 + x1*(s + 1/2) )
dy=-ds*(y1*(-s - 1/2) +2*s*y8 + y7*(-s+ 1/2)) = ds*(y7*(s - 1/2) - 2*s*y8 + y1*(s+ 1/2))
Conclusion: in general we can write
dx= ds*(x1*(s- 1/2) - 2*s*x2 + x3*(s+ 1/2))= a(s) ds
dy= ds*(y1*(s- 1/2) - 2*s*y2 + y3*(s+ 1/2))= b(s) ds
where (x1,y1) is the first node, (x2,y2) is the second node and (x3,y3) is the third node of
the ‘1D contour element’, and s is the parametric curvilinear coordinate going from -1 to
1.
Note that :
(s-1/2)= N1(s)’ since N1(s) = s*(s-1)/2 , N1, N2, N3 are the shape functions in 1D!
(s+1/2)= N3(s)’ N3(s)= s*(s+1)/2
-2s= N2(s)’ N2(s)= 1-s2
and
x(s)= N1(s)*x1+ N2(s)*x2+x3(s)*x3
x(s)= N1(s)*y1+ N2(s)*y2+x3(s)*y3
Therefore: the element length dL = sqrt(dx2+dy2) = sqrt(a2+ b2) ds
Any integral over the contour (A (x1,y1)-->B(x3,y3)) can be computed using integration
over the 1D reference element as:
𝐵 1
∮ 𝑓(𝑥, 𝑦)𝑑𝐿 = ∫ 𝑓(𝑥(𝑠), 𝑦(𝑠)) sqrt(a^2 + b^2) ds
𝐴 −1
The factor sqrt(a2+ b2) is the determinant of the geometric mapping from the curvilinear
real element to the 1D reference element.