Calculation of Jacobians for a finite element method simulation

Question

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}}

Answer 1

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}]