Calculating the Area of a Deformed Domain Using the Jacobian

General
1

Hey there,

I was wondering if anyone had any idea on how one could go about calculating more complex areas in FEnicSx? I’ve been stuck with trying to calculate the area of a deformed object using the Jacobian in FEnicSx.

My problem is that I want to calculate the area of a deformed object (A=?) with the untransformed area being a circle of unitary area (A_0 =1). The value of J will be updated at every time step in my linear elasticity problem.

QUESTION: How do we calculate Area(E), assuming D is the undeformed reference configuration (circle), E is the deformed domain, and J is the Jacobian:

Area(E)=\int \int_E \space dx \space dy =\int \int_D |J| \space du \space dv = ?

In legacy FEnicS, I believe it would have been calculated using:

assemble(J*dx)

However, been unable to understand how to do this with the dolfinx documentation to determine how J=ufl.det(F) and ufl.det objects interact in this new scenario. Potentially need to interpolate into some function space? I’ve also attached some code with the last few lines being roughly how I originally envisioned calculating the area (originally was doing this without J). Moreso here for anyone to see how I’m generating the mesh, measures, function space, and jacobian:

#----------------
# GENERAL IMPORTS
#----------------

import ufl
import numpy as np
import gmsh

from petsc4py import PETSc
from mpi4py import MPI

import dolfinx
from dolfinx import fem, mesh, plot, log, default_scalar_type, geometry
from dolfinx.io import gmshio
from dolfinx.fem.petsc import assemble_matrix, create_vector, assemble_vector, apply_lifting
import basix



#-------------------
# CREATE CIRCLE MESH
#-------------------

# Initialise the GMSH module
gmsh.initialize()

# Create the membrane and start the computations by the GMSH CAD kernel, to generate the relevant underlying data structures
# Create Circular Mesh with area of A0 = 1
membrane = gmsh.model.occ.addDisk(0, 0, 0, np.sqrt(1/np.pi), np.sqrt(1/np.pi)) # x, y, z, x-radius, y-radius
gmsh.model.occ.synchronize()

# Make membrane a physical surface for GMSH to recognise when generating the mesh
gdim = 2 # 2D Geometric Dimension of the Mesh
gmsh.model.addPhysicalGroup(gdim, [membrane], 0) # Dimension, Entity tag, Physical tag

# Generate 2D Mesh with uniform mesh size
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.023)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.023)
gmsh.model.mesh.generate(gdim)

# Create a mesh from the GMSH model
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)
domain.name = "initial_mesh"
gmsh.write("circle.msh")

# Finalize GMSH
gmsh.finalize()



#---------------------------------
# GENERAL MEASURES AND COORDINATES
#---------------------------------

mesh_degree = domain.geometry.cmap.degree

# Global Spatial Coordinates as a Vector-Valued Expression i.e. x[0] = x, x[1] = y
x_spatial = ufl.SpatialCoordinate(domain)

metadata = {"quadrature_degree": 4}

# Define Measure Over Cells of Mesh
dx = ufl.Measure('dx', domain=domain, metadata=metadata)

# Define Measure Over Boundary of Mesh
ds = ufl.Measure('ds', domain=domain, metadata=metadata)

#----------------------------------------
# LINEAR ELASTICITY VECTOR FUNCTION SPACE
#----------------------------------------

# Vector-Valued Function Space for Linear Elasticity
v_cg1 = basix.ufl.element(family="Lagrange", cell=domain.topology.cell_name(), degree=mesh_degree, shape=(domain.geometry.dim, ))
V = fem.functionspace(mesh=domain, element=v_cg1)

# Trial and Test Functions for Linear Elasticity Problem (Mechanics Vectors)
# Initial Conditions: Mechanical System is at Rest
u = ufl.TrialFunction(V)#fem.Function(V)
v = ufl.TestFunction(V)
u_n = fem.Function(V)
u_n.name = "u_n"

u_0 = fem.Constant(domain, 0.0)
u_n.sub(0).interpolate(lambda x: np.full((x.shape[1],), u_0))
u_n.sub(1).interpolate(lambda x: np.full((x.shape[1],), u_0))
u_n.sub(0).collapse()
u_n.sub(1).collapse()
u_n.x.scatter_forward()


#-------------------
# CALCULATE JACOBIAN
#-------------------
# Spatial Dimension
spa_dim = len(x_spatial)

# Identity Tensor
I = ufl.variable(ufl.Identity(spa_dim))

# Deformation Gradient Tensor
F = ufl.variable(I + ufl.grad(u)) 

F_inv = ufl.variable(ufl.inv(F)) # F^-1
F_inv_tra = ufl.variable(F_inv.T) # F^-T

# Jacobian Function
J = ufl.det(F)
J_n = J




#------------------------------
# EXAMPLE INCORRECT CALCULATION
#------------------------------

# Calculate Mesh Area
def mesh_area(mesh_domain, jacobian):
    area = fem.assemble_scalar(fem.form(fem.Constant(mesh_domain, 1.0)*jacobian*dx))
    return area

A=mesh_area(domain, J)

Cheers,
Volkan

2

mesh_area is correct (you can remove the fem.Constant)
You just need to replace u with a fem.Function. This function should be filled with the values that you want to represent the deformation. Being for instance the solution of an elasticity problem.
Note that if the solution is coming from a linearized elasticity, the mesh_area formula will work only at first order, and in this case det(J) = 1+div(u)

1 Like
3

Awesome! This is exactly what I was looking for :smiley: Thanks so much for the quick response and for fixing the problem bleyerj. This will be super useful for others wanting to calculate areas.