Vortices

Here is another code to solve Euler’s equation (inviscid flow) and plot the pathlines. This method was adapted from Stephen Montgomery-Smith’sEuler2D program. The vortex effects are found using the Biot-Savart law and the differential equations are solved using the Adams Bashforth method. Click here to download some Matlab code for this. Shown below is some Mathematica code:
(* runtime: 19 seconds *)
ToPoint[z_] := {Re[z], Im[z]};dt = 0.01; rcore = 0.1; dz0 = 0; Klist = {1, -1, 1, -1};
zlist = Join[{-1 - 0.5I, -1 + 0.5I, -0.5 - 0.5I, -0.5 + 0.5I}, Flatten[Table[x + I y, {x, -1, 1, 0.1}, {y, -1, 1, 0.1}], 1]];
paths = Map[# Table[1, {10}] &, zlist];
Do[dz = Map[Sum[zj = zlist[[j]]; If[# != zj, r2 = Abs[# - zj]^2; I Klist[[j]](zj - #)/r2 (1 -Exp[-r2/rcore^2]), 0], {j, 1, Length[Klist]}] &, zlist]; zlist += dt (1.5dz - 0.5dz0); dz0 = dz; paths = Table[Prepend[Delete[paths[[i]], -1], zlist[[i]]], {i, 1, Length[zlist]}];Show[Graphics[Map[Line[Map[ToPoint, #]] &, paths], PlotRange -> {{-1, 1}, {-1, 1}}, AspectRatio -> Automatic]], {85}];

Link: math paper by Stephen Montgomery-Smith, seems fairly advanced to me

Unstructured Tetrahedral Meshed Sphere

Here is an unstructured tetrahedral meshed sphere.

Links

• Mesh Generation – interesting paper by Per-Olof Persson
• Free Mesh Generators – click here for a complete list:
  • Gmsh – easy to use, by Christophe Geuzaine, has ability to use pictures to determine mesh size
  • DistMesh – for Matlab, here are some samples
  • EasyMesh – C++ program for 2D meshes
  • Triangle – 2D meshes
    CUBIT (not free) – impressive hexahedral mesh generator

Finite Element Method (FEM) Solution to Poisson’s equation on Triangular Mesh

Finite Element Method (FEM) Solution to Poisson’s equation on Triangular Mesh
solved in Mathematica 4.2, 5/24/06; mesh generated with Gmsh; old version: Matlab and DistMesh

Here is some code to solve Poisson’s equation on an unstructured grid of triangular elements using the Finite Element Method (FEM):
(* runtime: 0.1 second *)
n = 10; nodes = Flatten[Table[{x, y}, {y, -1, 1, 2.0/(n - 1)}, {x, -1, 1, 2.0/(n - 1)}], 1]; elements = Flatten[Table[{{i, i + 1, i + n}, {i + 1, i + n + 1, i + n}} + (j - 1) n, {j, 1, n - 1}, {i, 1, n - 1}], 2]; nnode = Length[nodes]; 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[(xy =Join[Transpose[nodes[[#]]], {{1, 1, 1}}];deta = Inverse[xy][[All, {1, 2}]]; Kglobal[[#, #]] +=Det[xy]deta.Transpose[deta]/2; F[[#]] += Det[xy]/6) &, 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]];

Here is some code for an interpolated plot:
(* runtime: 10 seconds *)
n = 275; x1 = y1 = -1; x2 = y2 = 1; image = Table[0, {n}, {n}];
xIntersect[{{x1_, y1_}, {x2_, y2_}}] := If[y1 == y2, Infinity, x1 + (y - y1)(x2 - x1)/(y2 - y1)];
Scan[(plist = nodes[[#]]; xlist = nodes[[#, 1]]; ylist = nodes[[#, 2]]; p1 = plist[[1]]; {dx2, dy2} = plist[[2]] - p1; {dx3, dy3} = plist[[3]] - p1; Do[y = y1 + (y2 - y1)(i - 1)/(n - 1); {xa, xb,xc} = Map[xIntersect, Partition[Sort[plist, #1[[2]] < #2[[2]] &], 2, 1, 1]]; jlist = Round[n({xc, If[Abs[xb - xc] < Abs[xa - xc], xb, xa]} - x1)/(x2 - x1)]; If[jlist =!= {}, Do[x = x1 + (x2 - x1)(j - 1)/(n - 1); {xi, eta} = ((x - p1[[1]]){dy3, -dy2} + (y - p1[[2]]){-dx3, dx2})/(dx2 dy3 - dx3 dy2); image[[i, j]] = phi[[#]].{1 - xi - eta, xi, eta}, {j, Max[1, Min[jlist]], Min[n, Max[jlist]]}]], {i, Max[1, Round[n(Min[ylist] - y1)/(y2 - y1)]], Min[n, Round[n(Max[ylist] - y1)/(y2 - y1)]]}]) &, elements];
ListDensityPlot[image, Mesh -> False, Frame -> False, ColorFunction -> PlotColor, Epilog -> Map[{Line[Map[({x, y} = nodes[[#, {1, 2}]]; n{(x - x1)/(x2 -x1), (y - y1)/(y2 - y1)}) &, #]]} &, elements]];