I work to solve PDE using FEM in the case P2 on Matlab. I try to correctly assemble load vector using quadratic Lagrange shape functions $$b_i =(f,\phi_i)=\sum_{q=1}^{nq}f(r_q,s_q)*\phi_{i}(r_q,s_q)*w_q*det(J(r_q,s_q)).$$ Here is what I did

[rspts,qwgts] = Gausspoints(precision); % quadrature rule
np      = size(p,2);                   % number of nodes
nt      = size(t,2);                   % number of elements

b = sparse(np,1);

for i=1:nt % loop over elements

nodes = t(1:6,i); % node numbers
x = p(1,nodes);   % node x-coordinates
y = p(2,nodes);   % node y-coordinates

[S,dSdx,dSdy,detJ]=Isopmap(x,y,r,s,@P2shapes);        % map

wxarea=qwgts(q)*detJ/2;                               % weight times area

end
b(nodes)       = b(nodes) +bK;
end



My problem... When I take PDE with source term f=0, the Plot of Numerical solution and exact solution gives the same result. Howerver if I put for example f(x,y)=x+y, I end with a plot difference between the two solutions.

Brief explanation I have found helpful information: Computation of stiffness matrix with variable coefficient

SOLUTION: It comes that I didn't understand properly the logic of the mapping (x,y) domain to (r,s) domain. The integral can be written like this: (exactly in the line of elements load vector )

 ...

x_physical= dot(x,S);
y_physical= dot(y,S);
bK=bK+S*f(x_physical,y_physical)'*wxarea;
...

• You haven't included a question mark (?) in your question; this makes it harder to know what you are asking. Feb 8 at 20:30
• @Richard, my question is : How to compute of load vector in Matlab for solving PDE using FEM in the case P2 ? There is the same topic for stiffness matrix in Python, pls see scicomp.stackexchange.com/questions/27420/… Feb 8 at 21:47