Problem with the determinant of the Jacobian

General
1

I am using dolfin 2019.

I tried to calculated the det of jacobian following these post,

And combine it with the computation of quad points

My code to generate the mesh, it is a 1d interval with two elements [0,0.75] and [0.75,1.0]:

##################
# Mesh, vertices and cells
A = 0
B = 0.75
C = 1
N_AB = 1  # Number of cells in AB
N_BC = 1  # Number of cells in BC

topological_dim = 1
geometrical_dim = 1

mesh = Mesh()
editor = MeshEditor()

vertices_a = np.linspace(A, B, N_AB + 1)
vertices_b = np.linspace(B, C, N_BC + 1)[1:]
vertices = np.hstack([vertices_a, vertices_b])


num_local_vertices = N_AB + N_BC + 1
num_local_cells = N_AB + N_BC

num_global_vertices = num_local_vertices
num_global_cells = num_local_cells

editor.open(mesh, "interval", topological_dim, geometrical_dim)
editor.init_vertices_global(num_local_vertices, num_global_vertices)
editor.init_cells_global(num_local_cells, num_global_cells)

# Add vertices
for i, vertex in enumerate(vertices):
    editor.add_vertex(i, np.array([vertex], dtype="double"))
# Add cells
for i in range(num_local_cells):
    editor.add_cell(i, np.array([i, i + 1], dtype="uint"))

# Close editor
editor.close()

mesh1d=mesh

cells_ = mesh.cells()
vertices = mesh.coordinates()

print("vertices")
for vertex in vertices:
    print(vertex)

print("cells")
for cell in cells_:
    print(cell)

output two cells, three vertices

vertices
[0.]
[0.75]
[1.]
cells
[0 1]
[1 2]

The code to see the quad points, weights and det of Jacobian:

# fristly get the quad points
import ffc
import dolfin
from ufl import Jacobian

degree=1

# mesh = UnitSquareMesh(1, 1)

el = FiniteElement("Quadrature", mesh.ufl_cell(),  degree=degree, quad_scheme="default")
V = FunctionSpace(mesh, el)
x = V.tabulate_dof_coordinates()

print("\n physical cel quad points")

for cell in cells(mesh):
    print("cell",cell.index(),"cell nodes",V.dofmap().cell_dofs(cell.index()),'\n', x[V.dofmap().cell_dofs(cell.index())],'\n')

print("\nreference cel, cell name is  " ,mesh.cell_name())
points, weights = ffc.fiatinterface.create_quadrature(mesh.cell_name(), degree=degree, scheme="default")
print("weights of ref cel: \n", weights)

print("\nJacobian")

Q = FunctionSpace(mesh, "Lagrange", degree)

J = abs(det(Jacobian(mesh)))
cell_jac =dolfin.project(J, Q)
print(cell_jac.vector().get_local())

output

 physical cel quad points
cell 0 cell nodes [0] 
 [[0.375]] 

cell 1 cell nodes [1] 
 [[0.875]] 


reference cel, cell name is   interval
weights of ref cel: 
 [1.]

Jacobian
[0.0625 0.625  0.8125]

My first question is , there are two cells in the mesh, why are there three Jacobian dets?

Second question is, since the standard element has length 1, and length of two cells are 0.75, 0.25, the Jacobian of these two cells should be 0.75 and 0.25, is it right?

2

You are choosing to project the jacobian, a cellwise constant into a first order continuous Lagrange space. This space has degrees of freedom at each vertex, and you introduce a numerical error when trying to approximate a discontinuous function with a continuous basis.

If you use Q = FunctionSpace(mesh, "DG", 0) you should Get the two values you are expecting.

This was the function space chosen in the original post.

1 Like
3

Thanks very much! Got it.