# Calculation of Jacobians for a finite element method simulation

I'm attempting to calculate symbolic expressions for the (Newton's method) Jacobian of a Galerkin finite element method using Mathematica. It would seem that this is a perfect application for symbolic differentiation tools but I can't get it working (probably due to my lack of Mathematica knowledge).

# Mathematical statement of the problem

Let $v(\mathbf{x})$ and $\psi_k(\mathbf{x})$ ($k=1,\ldots,N$) be simple(ish) polynomial functions of space. Let $\mathbf{m}(\mathbf{x}) = [m_1(\mathbf{x}),m_2(\mathbf{x}),m_3(\mathbf{x})]$ be a vector field.

Let $$\mathbf{m}(\mathbf{x}) = \sum_{k=1}^N \psi_k(\mathbf{x}) \mathbf{m}_k,$$ and $$r(\mathbf{x}) = \int_V (\mathbf{m} \times \nabla^2 \mathbf{m}) v(\mathbf{x}) dV.$$

(I have used $\nabla^2$ as the vector Laplacian $\nabla^2 \mathbf{m} = [ \nabla^2 m_1, \nabla^2 m_2, \nabla^2 m_3].$)

Calculate $\frac{\partial r}{\partial m_{i,k}}$ for all $i$ (directions of $\mathbf{m}$), and all $k \in {1,\ldots,N}$ (index for the $\psi$ functions).

# My code

I'm sure this can be written much more elegantly but as I said before I'm not too great with Mathematica. Anyway it should be good enough to give you the idea.

psi[x_, y_, z_] := {a[x, y, z], b[x, y, z]};
(* Really we just need to declare that this is a vector somehow...*)

m[x_, y_, z_] := {m1x*psi[x, y, z][[1]] + m2x*psi[x, y, z][[2]],
m1y*psi[x, y, z][[1]] + m2y*psi[x, y, z][[2]],
m1z*psi[x, y, z][[1]] + m2z*psi[x, y, z][[2]]}

VectorLaplacian[f[x,y,z]_] := {
D[f[x,y,z][[1]], {x, 2}] + D[f[x,y,z][[1]], {y, 2}] + D[f[x,y,z][[1]], {z, 2}],
D[f[x,y,z][[2]], {x, 2}] + D[f[x,y,z][[2]], {y, 2}] + D[f[x,y,z][[2]], {z, 2}],
D[f[x,y,z][[3]], {x, 2}] + D[f[x,y,z][[3]], {y, 2}] + D[f[x,y,z][[3]], {z, 2}]
}

r[x_, y_, z_] := (Cross[m[x, y, z], VectorLaplacian[m[x, y, z]]])*v[x,y,z]
(* I've left out the integration for now to keep things simpler *)

J[x_, y_, z_] := {
{D[r[x,y,z][[1]], m1x], D[r[x,y,z][[2]], m1x], D[r[x,y,z][[3]], m1x]},
{D[r[x,y,z][[1]], m1y], D[r[x,y,z][[2]], m1y], D[r[x,y,z][[3]], m1y]},
{D[r[x,y,z][[1]], m1z], D[r[x,y,z][[2]], m1z], D[r[x,y,z][[3]], m1z]},

{D[r[x,y,z][[1]], m2x], D[r[x,y,z][[2]], m2x], D[r[x,y,z][[3]], m2x]},
{D[r[x,y,z][[1]], m2y], D[r[x,y,z][[2]], m2y], D[r[x,y,z][[3]], m2y]},
{D[r[x,y,z][[1]], m2z], D[r[x,y,z][[2]], m2z], D[r[x,y,z][[3]], m2z]}
}

J[1, 1, 1]


But all I get is:

{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}}

-
There are some syntax issues here. If you change all the terms like psi[[1]][x, y, z] to psi[x, y, z][[1]] then at least it makes sense (this would be a[x,y,z]) in your example. Just making this change in the m[] definition gives values for r, i.e., you can ask for r[x,y,z] and you get something fairly large. The J function has a similar issue, perhaps m[[1]] should be m[x,y,z][[1]] (presumably you mean the first element of the polynomial m[x,y,z]). – bill s Aug 15 '12 at 3:46
Thanks, I've fixed those mistakes now but I'm still not getting a real answer. – dshepherd Aug 15 '12 at 10:45
@dshepherd, note that the definiton of your VectorLaplacian is not quite right. You may want to define a function and then test it to see if it works. e.g. evaluate VectorLaplacian[m[x, y, z]] to see that it does what you want. – user21 Aug 15 '12 at 14:21

If you want to do this for a specific N it is straightforward.

Let us load the calculus package.

<<VectorAnalysis


For readability let us use the formatting rules

 Format[Subscript[\[Psi], a__][x, y, z]] = Subscript[\[Psi], a];
Format[\[Nu][x, y, z]] = \[Nu];


Use cartesian coordinate for Laplacian and Cross product

  SetCoordinates[Cartesian[x, y, z]]


Let us define

 psi[k_][x_, y_, z_] = {Subscript[\[Psi], k, 1][x, y, z], Subscript[\[Psi], k, 2][x, y, z] ,
Subscript[\[Psi], k, 3][x, y, z] };


Now let's say we have N=2 for the sake of this example; Define

 m[x_, y_, z_] = Sum[psi[k][x, y, z]*{Subscript[m, 1, k], Subscript[m, 2, k],
Subscript[m, 3, k]}, {k, 1, 2}];


Defining the integrant:

 rint = \[Nu][x, y, z] Cross[m[x, y, z], Map[Laplacian, m[x, y, z]]];


Now mm is the set of variables:

 mm =
Table[{Subscript[m, 1, k], Subscript[m, 2, k], Subscript[m, 3,
k]}, {k, 2}]


This table will be your 'Jacobian'

  jac=Map[D[rint, #] &, mm, {2}];
jac // TableForm


If you want an expression for arbitrary N its also possible.

EDIT

You might want to consider the case where the fields Subscript[m,i,k] depends on {x,y,z} as well.

One should then add the formatting rule

Format[Subscript[m, a__][x, y, z]] = Subscript[m, a];


define

m[x_, y_, z_] = Sum[psi[k][x, y, z]*{Subscript[m, 1, k][x, y, z],
Subscript[m, 2, k][x, y, z], Subscript[m, 3, k][x, y, z]}, {k, 1, 2}]


And deal with derivatives by integration by part (subject to proper boundary conditions...); e.g.

Subscript[\[Psi], 1, 1][x, y, z] D[D[Subscript[m, 1, 1][x, y, z], {x,2}],
Subscript[m, 1, 1][x, y, z]] -> D[Subscript[\[Psi], 1, 1][x, y, z], {x,2}]

-
Of course, these days, one loads <<VectorAnalysis  instead of the <<Calculus ` context... – J. M. Aug 21 '12 at 20:21
No point reminding me how old I am! :-) – chris Aug 21 '12 at 20:47
Don't worry about it, it took me a while to get used to the new way they're doing add-on packages myself... :) – J. M. Aug 21 '12 at 20:51
Thanks, I actually ended up rechecking my solutions by hand and then forgot about this question. But it will definitely be useful later on so thank you very much :) – dshepherd Sep 12 '12 at 17:43