Question about an issue with assemble() and project() for equations defined as UFL Forms

Legacy DOLFIN
,
2

If you use assemble(strain_energy(u)*dx)) you compute the integral of the strain energy over the whole domain. This will evaluate the strain energy at quadrature points of each cell and accumulate the results across all cells of the domain.

With your second command, you compute the strain-energy of u as piecewise constant function using a projection. When you sum these up, you are not taking the cell area into account, and therefore you get the wrong result.

i.e. if you use:

DG_strain = project(strain_energy(u), FunctionSpace(mesh, "DG", 0))
cell_sizes = project(CellVolume(mesh), FunctionSpace(mesh, "DG", 0))
# same as print(f"strain energy using assemble() quad: {assemble(strain_energy(u)*Measure("dx", mesh, metadata={"quadrature_degree": 0}))}\n")
import numpy as np

print(f"strain energy using assemble(): {assemble(strain_energy(u)*dx)}\n")
print(
    f"strain energy using project() to DG0: {np.dot(cell_sizes.vector().get_local(), DG_strain.vector().get_local())}\n"
)

you get:

strain energy using assemble(): -0.008689827046779201

strain energy using project() to DG0: -0.0086898270467792
1 Like