0
$\begingroup$

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

bK=zeros(6,1); % elements load vector

for q=1:length(qwgts)                                  % quadrature loop   
r=rspts(q,1);                                         % quadrature r-coordinate
s=rspts(q,2);                                         % quadrature s-coordinate
[S,dSdx,dSdy,detJ]=Isopmap(x,y,r,s,@P2shapes);        % map

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

bK=bK+S*f(mean(r),mean(s))*wxarea; % elements load vector
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;
...
$\endgroup$
2
  • 1
    $\begingroup$ You haven't included a question mark (?) in your question; this makes it harder to know what you are asking. $\endgroup$
    – Richard
    Feb 8 at 20:30
  • $\begingroup$ @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/… $\endgroup$
    – A. AchbK
    Feb 8 at 21:47

1 Answer 1

0
$\begingroup$

I see where the thing fails. This problem has already been asked and answered before, see in Computation of stiffness matrix with variable coefficient.

I edited my question, problem solved. Thank you, it was a bit stupid, but as we say, the devil is in the details.

$\endgroup$
1
  • $\begingroup$ As it’s currently written, your answer is unclear. Please edit to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center. $\endgroup$
    – Community Bot
    Feb 9 at 16:46

Not the answer you're looking for? Browse other questions tagged or ask your own question.