Does Fenics have the ability to project a solution onto a coarse mesh while conserving the integral over the domain

Question

I am currently trying to solve a heating problem, but I need a very refined mesh to compute the heat deposition compared to the coarse mesh in the heat transfer model. I would therefore like to compute a highly detailed local heat deposition profile, using a restriction, and project this onto a coarse mesh for the bioheat calculation.

So far the restriction and projection works, but the integral of the Power is not conserved between the meshes. This is expected due to the piecewise polynomial structure of the fine mesh not being in the space spanned by the coarse mesh. I know that there are ways to do a conservative projection (interpolation) using a supermesh and suspect something like this may be implemented in the automatic meshing capabilities, but I do not know how to access them or use them for my application.

Would anyone be able to point me in the right direction or confirm this feature does not exist?

Thanks,

Sheldon

Answer 1

FEniCS is a programmable environment - you first need to identify the precise mathematical operation that you want to perform. If you know this and present it, then we can suggest how to implement it.

Comments

Hi Garth,

Thank you for replying. I have put together a simple example to explain more clearly my issue:

"""
Example of non-conservative projection:

The function Sin(5*pi*x[0])**2*Sin(5*pi*x[1])**2 is projected onto a basis of linear
Lagrange elements on the unit square with a fine mesh. An integral is
performed of the FEM representation over the whole domain to compute
the analogue of the energy. The FEM approximation is then projected
onto a coarse mesh and the same integral performed. The results of
the integrals do not agree as the algorithm is not conservative.

"""

from dolfin import *

mesh_fine = UnitSquare(20, 20)
mesh_coarse = UnitSquare(3, 3)

V = FunctionSpace(mesh_fine,'Lagrange',1)
W = FunctionSpace(mesh_coarse,'Lagrange',1)
u = Function(V)
f = Expression('pow(sin(5*pi_*x[0]),2)*pow(sin(5*pi_*x[1]),2)',pi_=pi)

u = interpolate(f, V)

energy = u*dx
E = assemble(energy)
print "The fine mesh integral is: %g" % E

v = Function(W)
v = project(u,W)

energy = v*dx
E = assemble(energy)
print "The coarse mesh integral is: %g" % E

If you execute this script the two values output for the integral are different. I would really like to know if there is a function that exists in fenics to compute a projection such that the integrals are equal? Perhaps it is in the adaptive meshing libraries?

Thanks,

Sheldon