function [Elem]=ElemStiffElasMatP1Bb(q1,q2,q3,area,lambda,mu)
% function [Elem]=ElemStiffElasMatP1Bb(q1,q2,q3,area,lambda,mu)
%  Computation of the element stiffness elasticity matrix for
%  P1-Lagrange method. 
%  The method for numbering the degrees of freedom is local block numbering.
%
%  Hooke's matrix :
%  C=[L + 2*M       L       0]
%    [      L L + 2*M       0]
%    [      0       0       M]
%  Numbering of local points in reference element is :
%    P=[(0, 0), (1, 0), (0, 1)]
%
% Parameters:
%  q1    	: array of coordinates of the first point of the triangle
%  q2    	: array of coordinates of the second point of the triangle
%  q3    	: array of coordinates of the third point of the triangle
%  area 	: triangle area
%  lambda 	: first  Lame coefficient in Hooke's law
%  mu     	: second Lame coefficient in Hooke's law
%
% Return values:
%  Elem 	: element stiffness elasticity matrix, 6-by-6 matrix
%
% Example: 
%    q1=[0;0];q2=[1;0];q3=[0;1];
%    area=1/2.;
%    lambda=1.; mu=1.;
%    KElem=ElemStiffElasMatP1Bb(q1,q2,q3,area,lambda,mu);

%
% Copyright (C) 2013  CJS (LAGA)
%   see README for details

u=q2-q3;
v=q3-q1; 
w=q1-q2;
% Matrice de Hooke
C=[lambda+2*mu,lambda,0;lambda,lambda+2*mu,0;0,0,mu];
% Matrice des déformations (x par 2*area)
B=[u(2),v(2),w(2),0,0,0; ...
   0,0,0,-u(1),-v(1),-w(1); ...
   -u(1),-v(1),-w(1),u(2),v(2),w(2)];
Elem=B'*C*B/(4*area);
