# How to obtain Jacobian matrix for domain mapping purpose?

**URL:** <https://fenicsproject.discourse.group/t/how-to-obtain-jacobian-matrix-for-domain-mapping-purpose/7961>\
**Category:** variational formulation\
**Created:** [March 31, 2022, 9:22am UTC](https://fenicsproject.discourse.group/t/how-to-obtain-jacobian-matrix-for-domain-mapping-purpose/7961 "2022-03-31T09:22:54Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![CuriousGuy](https://avatars.discourse-cdn.com/v4/letter/c/b5a626/32.png) [@CuriousGuy](https://fenicsproject.discourse.group/u/CuriousGuy)\
**Post date:** [March 31, 2022, 9:22am UTC](https://fenicsproject.discourse.group/t/how-to-obtain-jacobian-matrix-for-domain-mapping-purpose/7961/1 "2022-03-31T09:22:54Z")

</div>

Dear FEniCS community,

I am working on a problem where my domain, i.e., the geometric model has to be parameterized. To tackle such a challenge, I am trying to map different versions/variants of this geometric model onto a reference geometry. Mapping function is a Laplace’s equation, so by solving a Laplace’s equation on reference domain, I get the new version of my domain. Now, the goal is to solve some fluid dynamics problem on this new version of domain, however, the idea is to formulate everything on reference domain to save the cost of meshing again and again for every new domain.

So, in order to formulate everything in reference domain, I need Jacobian matrix, which in this case is simply the gradient of solution of Laplace’s equation. Please have a look at the MWE to have an idea

```auto
from dolfin import *

geometry = Rectangle(Point(0., 0.), Point(4., 3.))  
mesh = generate_mesh(geometry, 10)

class Top(SubDomain):
    def inside(xself,x,on_boundary):
        return x[1] > 3.0 - DOLFIN_EPS
class Bottom(SubDomain):
    def inside(self,x,on_boundary):
        return x[1] < DOLFIN_EPS
class Left(SubDomain):
    def inside(xself,x,on_boundary):
        return x[0] < DOLFIN_EPS
class Right(SubDomain):
    def inside(xself,x,on_boundary):
        return x[0] > 4.0 - DOLFIN_EPS

boundaries = MeshFunction("size_t", mesh, 1)
boundaries.set_all(0)

Left().mark(boundaries, 0)
Right().mark(boundaries, 1)
Top().mark(boundaries, 3)
Bottom().mark(boundaries, 2)

n = FacetNormal(mesh)

ds = Measure("ds", domain=mesh, subdomain_data=boundaries)

V = VectorFunctionSpace(mesh, "Lagrange", 1)

u_trial = TrialFunction(V)
v_test = TestFunction(V)

exp = Expression(('0','0'),element=V.ufl_element())
g = exp

# Define variational problem 
a = inner(grad(u_trial), grad(v_test))*dx
L = inner(g,v_test)*dx
# Boundary Conditions
bc = [(DirichletBC(V,Constant((0,0)),boundaries,0)),(DirichletBC(V,Constant((1,0)),boundaries,1))]

# Compute solution
sol = Function(V)
solve(a == L, sol, bc, solver_parameters={'linear_solver':'mumps'})

xderivative = sol.dx(0)
yderivative = sol.dx(1)

duoverdx = project(xderivative,V)
duoverdy = project(yderivative,V)

J = PETSc.Mat().createDense((len(u_trial),len(sol.vector())),array=([duoverdx.vector().vec().getArray()],[duoverdy.vector().vec().getArray()])).transpose()

```

This was my approach for obtaining the Jacobian matrix, however, the problem with this approach is that it is not a square matrix. I believe that the values in the matrix are correct, but this matrix is not sparse and square. What I could do is, I could manually make a square matrix of zeros and assign these values at the right positions according to the node numbers. But, this approach would very cumbersome and inefficient.

So, my question is, is there an efficient way to obtain the Jacobian matrix of such a system?

Thanks!

---

<div class="post-metadata">

**Author:** ![CuriousGuy](https://avatars.discourse-cdn.com/v4/letter/c/b5a626/32.png) [@CuriousGuy](https://fenicsproject.discourse.group/u/CuriousGuy)\
**Post date:** [April 1, 2022, 8:08am UTC](https://fenicsproject.discourse.group/t/how-to-obtain-jacobian-matrix-for-domain-mapping-purpose/7961/2 "2022-04-01T08:08:12Z")

</div>

Hi @dokken ,

Do you think maybe you could have a look at it? I will really appreciate it

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [April 1, 2022, 10:01pm UTC](https://fenicsproject.discourse.group/t/how-to-obtain-jacobian-matrix-for-domain-mapping-purpose/7961/3 "2022-04-01T22:01:05Z")

</div>

It sounds like what you’re trying to do is similar to what I covered in a short tutorial that I prepared a while back. You might take a look at the notes, code, and linked video from [this repository](https://github.com/david-kamensky/uli-workshop-examples-2022). The only non-obvious point that is missing is transformation of area elements, e.g., if you have a `ds` term in the formulation that you want to solve on the mapped domain. Area elements would need to transform using Nanson’s formula (as given, e.g., [here](https://en.wikiversity.org/wiki/Continuum_mechanics/Volume_change_and_area_change), in the context of referring solid mechanics equations back to a reference domain). My graduate course notes, available [here](https://david-kamensky.eng.ucsd.edu/teaching/mae-207-fea-for-coupled-problems), cover fluid–structure interaction using a similar formulation, but this is significantly more complicated.

---

<div class="post-metadata">

**Author:** ![CuriousGuy](https://avatars.discourse-cdn.com/v4/letter/c/b5a626/32.png) [@CuriousGuy](https://fenicsproject.discourse.group/u/CuriousGuy)\
**Post date:** [April 5, 2022, 7:52am UTC](https://fenicsproject.discourse.group/t/how-to-obtain-jacobian-matrix-for-domain-mapping-purpose/7961/4 "2022-04-05T07:52:42Z")

</div>

Dear @kamensky ,

Thank you for your guidance. The material you have referred me to was helpful! 🙂

---

<div class="post-metadata">

**Author:** ![CuriousGuy](https://avatars.discourse-cdn.com/v4/letter/c/b5a626/32.png) [@CuriousGuy](https://fenicsproject.discourse.group/u/CuriousGuy)\
**Post date:** [April 8, 2022, 8:28am UTC](https://fenicsproject.discourse.group/t/how-to-obtain-jacobian-matrix-for-domain-mapping-purpose/7961/5 "2022-04-08T08:28:04Z")

</div>

Dear @kamensky ,

One short query though. In your [notes](https://github.com/david-kamensky/uli-workshop-examples-2022), where you mapped domain from polar (reference config.) to rectangular coordinates (current config.), you formulate everything in the reference config. but in terms of Jacobian, in order to solve for the current config. So in principle, the solution obtained on this reference config. should be the same solution if I were to formulate everything on current config and solve it .

For the purpose of verifying the results, I recently formulated everything in current config., where the domain (after mapping) was bounded by the points `(1,0),(0,1)(2,0)(0,2)`. However, the solution on this domain does not match the solution on reference domain. Here are the results from both the domains:

`Note : For these results, kappa = Constant(2), f = Constant(2). The rest is same`

![polar](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/d/df903bcf147e1ceb0c961cbcb68c15538a398c3e.jpeg)  
 ![rectangular](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/5/5e805579ea281f3f0bb2d177cf59198274b752f0.jpeg)

In my opinion, the problem could be and as you also mentioned, the transformation of area. For your formulation in reference config.(or polar coordinates), you used only the determinant of deformation gradient for replacing `dx`(which for this 2D problem is the area integral and not volume integral). However, according to Nanson’s formula, we need to have other terms as well. I did try to account for those terms, but was not able to solve and was getting error in formulating them. My questions are:

1- Am I doing\thinking anything fundamentally wrong here?  
2- If no, then could you please direct me to formulation of area transformations in FEniCS?

Thank you!

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [April 8, 2022, 12:45pm UTC](https://fenicsproject.discourse.group/t/how-to-obtain-jacobian-matrix-for-domain-mapping-purpose/7961/6 "2022-04-08T12:45:08Z")

</div>

It’s not clear what you’ve done for the plots in your previous post, but a correct modification of the example code in `poisson-polar.py` to solve on the physical annular domain would look something like the following:

```python
from fenics import *
import math
M = 16; N = 32
mesh = RectangleMesh(Point(1,0),
                     Point(2,math.pi/2),
                     M,N)
xi = SpatialCoordinate(mesh)
r = xi[0]
theta = xi[1]
x = as_vector([r*cos(theta),
               r*sin(theta)])

# Deform the mesh to the physical domain
# using a projection of the analytical
# mapping's displacement onto a piecewise
# linear vector-valued function space.
Vvec = VectorFunctionSpace(mesh,"CG",1)
u_hat = project(x-xi,Vvec)
ALE.move(mesh,u_hat)

# Solve with no change of variables.
kappa = Constant(2)
f = Constant(2)
V = FunctionSpace(mesh,"CG",1)
u = TrialFunction(V)
v = TestFunction(V)
a = kappa*dot(grad(u),grad(v))*dx
L = f*v*dx

bc = DirichletBC(V,Constant(0),"on_boundary")
u = Function(V)
solve(a==L,u,bc)
u.rename("u","u")
File("u.pvd") << u

```

The resulting solution is nearly identical to the one on the reference domain (up to only slight differences from the piecewise-linear approximation of the mapping needed to mesh the curved domain geometry in physical space).

> For your formulation in reference config.(or polar coordinates), you used only the determinant of deformation gradient for replacing `dx` (which for this 2D problem is the area integral and not volume integral).

Note that I also transformed spatial derivatives, using the chain rule.

> However, according to Nanson’s formula, we need to have other terms as well. I did try to account for those terms, but was not able to solve and was getting error in formulating them.

Nanson’s formula would only be necessary in formulations with boundary integrals (i.e., using the `ds` measure). It looks like you’re still solving with strongly-enforced Dirichlet boundary conditions, so there is no need to use Nanson’s formula.

---

<div class="post-metadata">

**Author:** ![CuriousGuy](https://avatars.discourse-cdn.com/v4/letter/c/b5a626/32.png) [@CuriousGuy](https://fenicsproject.discourse.group/u/CuriousGuy)\
**Post date:** [April 8, 2022, 1:02pm UTC](https://fenicsproject.discourse.group/t/how-to-obtain-jacobian-matrix-for-domain-mapping-purpose/7961/7 "2022-04-08T13:02:40Z")

</div>

Yup, the solution now makes sense. I made a mistake in mapping to the rectangular coordinate system. So my domain was incorrect and hence the result was wrong.

Regarding the application of Nanson’s formula, what I have understood from your explanation is that we need it only for boundary integrals, e.g., for implementing Neumann BCs in weak form. In other words, it does not matter if we are working with a 2D or 3D domain, what matters is that it has to be used only for boundary intergrals, i.e, for a line in case of 2D or for a surface/plane in case of 3D.

Am I right?

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [April 8, 2022, 2:25pm UTC](https://fenicsproject.discourse.group/t/how-to-obtain-jacobian-matrix-for-domain-mapping-purpose/7961/8 "2022-04-08T14:25:01Z")

</div>

Yes, Nanson’s formula would be used for boundary (or interior facet, i.e., `dS`) integrals.

---

<div class="post-metadata">

**Author:** ![CuriousGuy](https://avatars.discourse-cdn.com/v4/letter/c/b5a626/32.png) [@CuriousGuy](https://fenicsproject.discourse.group/u/CuriousGuy)\
**Post date:** [April 8, 2022, 2:36pm UTC](https://fenicsproject.discourse.group/t/how-to-obtain-jacobian-matrix-for-domain-mapping-purpose/7961/9 "2022-04-08T14:36:27Z")

</div>

Thank you very much! 🙂

Just one last thing, could you also direct me to implementation of Nanson’s formula for boundary integrals in FEniCS?

Thanks a lot again!

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [April 8, 2022, 3:12pm UTC](https://fenicsproject.discourse.group/t/how-to-obtain-jacobian-matrix-for-domain-mapping-purpose/7961/10 "2022-04-08T15:12:30Z")

</div>

For a simple demo with pure Neumann BCs, we can tack on a reaction term to keep things well-posed, i.e., solve

-\nabla\cdot(\kappa\nabla u) + u = f\text{ ,}

in \Omega, subject to the boundary condition \nabla u\cdot\mathbf{n} = g on \partial \Omega. Posed directly on the deformed mesh, this would involve the following modification of the previous example:

```python
kappa = Constant(2)
f = Constant(2)
g = Constant(2)
V = FunctionSpace(mesh,"CG",1)
u = TrialFunction(V)
v = TestFunction(V)
a = (kappa*dot(grad(u),grad(v)) + u*v)*dx
L = f*v*dx + g*v*ds

u = Function(V)
solve(a==L,u)

```

To pull this problem back to the rectangular reference domain, one use:

```python
kappa = Constant(2)
f = Constant(2)
g = Constant(2)
V = FunctionSpace(mesh,"CG",1)
u = TrialFunction(V)
v = TestFunction(V)
a = (kappa*dot(grad_x(u),grad_x(v)) + u*v)*det(dx_dxi)*dxi
L = f*v*det(dx_dxi)*dxi

# Nanson's formula with "deformation gradient"
# of dx_dxi:
N = FacetNormal(mesh)
F = dx_dxi
J = det(dx_dxi)
da_dA_n = J*(inv(F).T)*N
da_dA = sqrt(dot(da_dA_n,da_dA_n))

# Use to trasnform surface integral over
# reference boundary (dA, corresponding to
# UFL ds in this setup) to a surface
# integral over the deformed boundary:
L += g*v*da_dA*ds

u = Function(V)
solve(a==L,u)

```

---

<div class="post-metadata">

**Author:** ![CuriousGuy](https://avatars.discourse-cdn.com/v4/letter/c/b5a626/32.png) [@CuriousGuy](https://fenicsproject.discourse.group/u/CuriousGuy)\
**Post date:** [April 11, 2022, 7:53am UTC](https://fenicsproject.discourse.group/t/how-to-obtain-jacobian-matrix-for-domain-mapping-purpose/7961/11 "2022-04-11T07:53:14Z")

</div>

Ah, I see! So the `n` of current config. doesn’t show up at all in formulation. We just eliminate it by taking the dot product of `da_dA_n` by itself to get rid of `n`. And that’s how we get the ratio two different areas, which is the scaling factor for change of variable/domain for boundary term integration in the formulation.

Many thanks! Really appreciate your help 🙂