06
Jun
07

### FEM Solution to Poisson’s equation on Quad Mesh

FEM Solution to Poisson’s equation on Quad Mesh – Mathematica 4.2, AutoCAD 2000, 6/10/07
Here is some code to solve Poisson’s equation on an unstructured grid of quadrilateral elements:
```(* runtime: 0.05 second *) n = 10; nodes = Flatten[Table[{x, y}, {x, -1, 1,2.0/(n - 1)}, {y, -1, 1, 2.0/(n - 1)}], 1]; nnode = Length[nodes]; elements =Flatten[Table[{i, i + 1, i + n + 1, i + n} + (j - 1) n, {j, 1, n - 1}, {i, 1, n - 1}], 1]; nelem = Length[elements]; fixed = Map[Position[nodes, #][[1, 1]] &, Select[nodes, Max[Abs[#]] == 1 &]]; Kglobal = Table[0, {nnode}, {nnode}]; F = phi = Table[0, {nnode}]; Scan[(plist = nodes[[#]]; dphi = {plist[[2]] - plist[[1]],plist[[4]] - plist[[1]]}; B = Inverse[Transpose[dphi].dphi]; C1 = {{2, -2}, {-2, 2}}B[[1, 1]] + {{3, 0}, {0, -3}}B[[1, 2]] + {{2, 1}, {1,2}}B[[2, 2]];C2 = {{-1, 1}, {1, -1}}B[[1, 1]] + {{-3, 0}, {0, 3}} B[[1, 2]] + {{-1, -2}, {-2, -1}}B[[2, 2]]; Kglobal[[#, #]] += Det[dphi]Join[Thread[Join[C1, C2]], Thread[Join[C2, C1]]]/6; F[[#]] += Det[Join[{{1, 1, 1}}, Transpose[nodes[[#[[{1, 2, 3}]], {1, 2}]]]]]/4) &, elements]; free = Complement[Range[nnode], fixed]; phi[[free]] = LinearSolve[Kglobal[[free, free]], F[[free]]]; Mean[x_] := Plus @@ x/Length[x]; PlotColor[x_] := Hue[2(1 - x)/3]; Show[Graphics[Table[{PlotColor[Mean[phi[[elements[[i]]]]]/Max[phi]], Polygon[nodes[[elements[[i]]]]]}, {i, 1, nelem}], AspectRatio -> 1]];```

## Welcome !

You will find here some of my favorite hobbies and interests, especially science and art.

I hope you enjoy it!

Subscribe to the RSS feed to stay informed when I publish something new here.

I would love to hear from you! Please feel free to send me an email : bugman123-at-gmail-dot-com