Mass Lumping of a coupled system

Question

Relating to http://fenicsproject.org/qa/4284/mass-matrix-lumping , how can I obtaining the mass matrix lumping of a coupled system? If I have two inner form, a1 and a2, how can I obtain the mass lumped matrix A of a = a1 + a2 ?

Thank a lot!
Cristina

Answer 1

Hi,

If you have a system of N equations, you can do:

V = ...
W = MixedFunctionSpace([V]*N)
psi = TrialFunctions(W)  # solution
v   = TestFunctions(W)

a_lumped = 0

for i in range(0, N):
    a_lumped += psi[i]*v[i]*dx(mesh)

# Class for mass matrix lumping
class LumpingClass(Expression):
    def eval(self, values, x):
        for i in range(0, N):
           values[i] = 1.0
    def value_shape(self):
        return (N,)

const = LumpingClass()

mass_action_form = action(a_lumped, const)
M_lumped = assemble(a_lumped) 
M_lumped.zero()
M_lumped.set_diagonal(assemble(mass_action_form)) # M_lumped is now lumped

Comments

Thanks a lot! Very clear!!!