Hi All, I would like to solve Poisson equation using FEM.

$-\Delta u=f,\quad(x,y)\in\Omega=(0,1)\times(0,1)$

$u=0,\quad \Gamma=\partial\Omega$

where $ f(x,y)=2?^2sin(?x)sin(?y)$ using FEM and grid as belove. (Exact solution is $u(x,y)=sin(?x)sin(?y)$)

Here is my code

ClearAll["Global`*"]

n = 8;

(*vertex coordinates of triangles*)
K = 1/n Flatten[
    Drop[Table[{{{i + 1, j + 1}, {i, j + 1}, {i, j}}, {{i + 1, 
         j + 1}, {i, j}, {i + 1, j}}, {{i + 1, j + 1}, {i + 1, 
         j}, {i + 2, j + 1}}, {{i + 1, j + 1}, {i + 2, j + 1}, {i + 2,
          j + 2}}, {{i + 1, j + 1}, {i + 2, j + 2}, {i + 1, 
         j + 2}}, {{i + 1, j + 1}, {i + 1, j + 2}, {i, j + 1}}}, {j, 
       0, n}, {i, 0, n}], -2, -2], 1];

(*area of triangles*)
A = 1/(2 n^2);

a[i_, j_] := 
 K[[i, j, 2, 1]] K[[i, j, 3, 2]] - K[[i, j, 3, 1]] K[[i, j, 2, 2]]
b[i_, j_] := K[[i, j, 2, 2]] - K[[i, j, 3, 2]]
c[i_, j_] := K[[i, j, 3, 1]] - K[[i, j, 2, 1]]

\[Phi][i_, j_, x_, y_] := 1/(2 A) (a[i, j] + b[i, j] x + c[i, j] y)

(*basis function*)
\[Phi] = Table[\[Phi][i, j, x, y], {i, Length@K}, {j, 6}] // Expand;

(*main diagonal of stiffness matrix*)
d1 = Total /@ 
   Table[A (Grad[\[Phi][[i, j]], {x, y}].Grad[\[Phi][[i, j]], {x, 
         y}]), {i, (n - 1)^2}, {j, 6}];

(*second diagonal of stiffness matrix*)
d2 = Riffle[
   Join @@ Table[
     A (Grad[\[Phi][[i + j - 1, 3]], {x, y}].Grad[\[Phi][[i + j, 
           1]], {x, y}]) + 
      A (Grad[\[Phi][[i + j - 1, 4]], {x, y}].Grad[\[Phi][[i + j, 
           6]], {x, y}]), {j, 1, 3 (n - 1), 3}, {i, n - 2}], 0, n - 1];

(*fourth diagonal of stiffness matrix*)
d4 = Table[
   A (Grad[\[Phi][[i, 6]], {x, y}].Grad[\[Phi][[i + 3, 2]], {x, y}]) +
     A (Grad[\[Phi][[i, 5]], {x, y}].Grad[\[Phi][[i + 3, 3]], {x, 
         y}]), {i, (n - 1)^2 - 3}];

(*fifth diagonal of stiffness matrix*)
d5 = Riffle[
   Join @@ Table[
     A (Grad[\[Phi][[i + j - 1, 5]], {x, y}].Grad[\[Phi][[i + j + 3, 
           1]], {x, y}]) + 
      A (Grad[\[Phi][[i + j - 1, 4]], {x, y}].Grad[\[Phi][[i + j + 3, 
           2]], {x, y}]), {j, 1, 3, 2}, {i, n - 2}], 0, n - 1];

(* stiffness matrix*)
stiffMatrix = 
  SparseArray[{Band[{1, 1}] -> 
     d1, {Band[{1, 2}], Band[{2, 1}]} -> {d2, d2}, {Band[{1, 4}], 
      Band[{4, 1}]} -> {d4, d4}, {Band[{1, 5}], Band[{5, 1}]} -> {d5, 
      d5}}];
MatrixForm@stiffMatrix;

(*right hand side of poisson equation*)
f[x_, y_] := 2 \[Pi]^2 Sin[\[Pi] x] Sin[\[Pi] y]

(*exact solution*)
exact[x_, y_] := Sin[\[Pi] x] Sin[\[Pi] y]

(*load vector version 1*)
loadVector = 
  A Total /@ (Table[
        f @@ K[[i, j, 1]]/6 +  f @@ K[[i, j, 2]]/12 +  
        f @@ K[[i, j, 3]]/12, {i, Length@K}, {j, 6}] // N);

(*load vector version 2*)
loadVector = f @@@ K[[All, 1, 1]]/n^2 // N;

(*center of triangles*)
c = Table[Total@K[[i, j]]/3, {i, Length@K}, {j, 6}] // N;

(*load vector version 3*)
loadVector = A/3 Table[Total@(f @@@ c[[i]]), {i, Length@c}] // N;

(*numeric solution*)
uh = LinearSolve[stiffMatrix, loadVector] // N;

(*exact solution at nodes*)
u = exact @@@ K[[All, 1, 1]] // N;

(*error*)
Sqrt[A Total[#^2 & /@ (u - uh)]]

But when I change n=8 to n=16 error increases which is not correct. Stiffness matrix is right form, I believe problem comes from load vector. That's why I tried couple version of approximation for right hand side but no luck. I follow this website page 219 http://www4.ncsu.edu/~zhilin/TEACHING/MA587/chap9.pdf

Any suggestions? Thanks in advance.

Attachments: