[C++] Compute Jacobian matrix of function

Question

Hello,

In my .cpp file I have a function $F:\mathbb{R}^2 \rightarrow \mathbb{R}^2$ to which I would like to compute its Jacobian matrix and determinant in the usual partial derivative sense, i.e.

$$\nabla_x F = \left(\begin{array}{cc}
\frac{\partial F^1}{\partial x_1} & \frac{\partial F^1}{\partial x_2} \ \newline
\frac{\partial F^2}{\partial x_1} & \frac{\partial F^2}{\partial x_2}
\end{array}\right) \in \mathbb{R}^{2\times 2}.$$

Then, I need to access the 4 values separately, as I have to pass them as coefficients in the .ufl file to solve a problem depending on $\nabla_x F$.

So, if my F function is defined as, say,

class F : public Expression
{
  public:
  void eval(Array<double>& values, const Array<double>& x) const
  {     
    R1 = 1.0;
    R2 = 2.0;
    values[0] = R1 * x[0];
    values[1] = R2 * x[1];
  }
};

its Jacobian is extremely easy to be calculated by hand, but for (far) more complicated functions I need a way to make it work for the general case.

Thanks for the help.

Answer 1

Using this question I was able to construct a solution:
Here I used $$F(x,y) = sin(x) + sin(y)$$ and calculated the numerical derivatives.

---

main.cpp

#include <dolfin.h>
#include <mshr.h>
#include "Poisson.h"
#include "jpcGrad.h"

using namespace dolfin;

// For Gradient
class GradClass : public Expression
{
public:
    GradClass(Function &cref)
    {
        c = &cref;
    }
    void eval(Array<double>& values, const Array<double>& x) const
    {
        Array<double> myval(1);
        c->eval(myval, x);
        values[0] = myval[0];
        cout << values[0];
    }
private:
    Function *c;
};

class test : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {

    values[0] = sin(x[0]) + sin(x[1]);
  }
};

class cos_x : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = cos(x[0]);
  }
};

class cos_y : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = cos(x[1]);
  }
};

int main()
{
  // Create mesh and function space
    auto rectangle = std::make_shared<mshr::Rectangle>(dolfin::Point(-5.0, 5.0, 0.0), dolfin::Point(5.0, -5.0, 0.0));
    auto mesh = mshr::generate_mesh(rectangle, 100);

    auto V = std::make_shared<Poisson::FunctionSpace>(mesh);

  //------------------------------------------------------

  test test_;
  cos_x cosx;
  cos_y cosy;

  auto G = std::make_shared<jpcGrad::FunctionSpace>(mesh);
    jpcGrad::LinearForm L_(G);
    jpcGrad::BilinearForm a_(G,G);

    Function w(V);
    Function w1(V);
    Function w2(V);

  w.interpolate(test_);
  w1.interpolate(cosx);
  w2.interpolate(cosy);

    auto gradf = std::make_shared<GradClass>(w);
    L_.solution = gradf;
    Function v(G);
    solve(a_ == L_, v);

  //------------------------------------------------------

  // Plot solution
  plot(w,"u = sin x + sin y");
  plot(v[0],"du/dx");
  plot(v[1],"du/dy");
  plot(w1,"cos(x)");
  plot(w2,"cos(y)");
  interactive();

  return 0;
}

---

jfcGrad.ufl

element     = FiniteElement("Lagrange", triangle, 2)
velement    = VectorElement("Lagrange", triangle, 2)

solution    = Coefficient(element)   
v           = TestFunction(velement)  # Test Function
u           = TrialFunction(velement) # Trial Function

a = inner(u, v)*dx 
L = inner(grad(solution),v)*dx

---

Poisson.ufl

element = FiniteElement("Lagrange", triangle, 2)