lumped mass

On Mon, Aug 20, 2012 at 10:45:38AM -0000, Arnd Flatten wrote:
> New question #206318 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/206318
>
> Hi,
>
> I simply want to achieve the (row-summed) lumped mass from any
> consistent mass matrix defined within FEniCS. That should be some
> simple operation, but i could not find any hint (manual/ faq) on an
> implemented way - can you provide any suggestions or guide me in the
> right direction, please?
>
> the system of eqn to solve for du (solution increment) is: Ml * du =
> dt * rhs(u_n) + (Mc-Ml) * u_n with Ml the lumped mass matrix, Mc the
> consistent mass matrix, dt the timestep and u_n the solutiona t
> previous time increment

If you are using the uBLAS linear algebra backend, it already has
support for lumping: M.lump(m) if M is a uBLASMatrix and m is a
uBLASVector.

In general, create a vector of 1s and do a matrix-vector
multiplication. Here's how it's done in C++ (from the file SingularSolver.cpp):

boost::shared_ptr<GenericVector> ones = A.factory().create_vector();
boost::shared_ptr<GenericVector> m = A.factory().create_vector();
ones->resize(N);
*ones = 1.0;
m->resize(N);
m->zero();
M->mult(*ones, *m);

--
Anders

A better approach is to assemble the action of the matrix on a unit vector. The assembled vector will contain the diagonals.

I've been asked this before, so we should add an example.

On Tue, Aug 21, 2012 at 10:01:09AM -0000, Garth Wells wrote:
> Question #206318 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/206318
>
> Garth Wells proposed the following answer:
> A better approach is to assemble the action of the matrix on a unit
> vector. The assembled vector will contain the diagonals.

Good point.

m = assemble(action(u*v*dx))

--
Anders

On Tue, Aug 21, 2012 at 01:36:50PM +0200, Anders Logg wrote:
> On Tue, Aug 21, 2012 at 10:01:09AM -0000, Garth Wells wrote:
> > Question #206318 on FEniCS Project changed:
> > https://answers.launchpad.net/fenics/+question/206318
> >
> > Garth Wells proposed the following answer:
> > A better approach is to assemble the action of the matrix on a unit
> > vector. The assembled vector will contain the diagonals.
>
> Good point.
>
> m = assemble(action(u*v*dx))

I forgot to set the coefficient. This should work:

m = assemble(action(u*v*dx), coefficients=[Constant(1)])

--
Anders

Thank you v much for your responses.

In either way, I understand to get a vector m representing the diagonal terms of the lumped mass matrix. Can I use this assembled vector to solve for my unknowns? Or do i need to convert the vector back on the diagonal of a matrix which is zero on all off-diagonal terms, just in order to use a solve() statement?! The latter would be somewhat awkward so I hope to avoid this if possible.

On 08/21/2012 02:11 PM, Arnd Flatten wrote:
> Question #206318 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/206318
>
> Status: Answered => Open
>
> Arnd Flatten is still having a problem:
> Thank you v much for your responses.
>
> In either way, I understand to get a vector m representing the diagonal
> terms of the lumped mass matrix. Can I use this assembled vector to
> solve for my unknowns? Or do i need to convert the vector back on the
> diagonal of a matrix which is zero on all off-diagonal terms, just in
> order to use a solve() statement?! The latter would be somewhat awkward
> so I hope to avoid this if possible.

You need a diagonal matrix (not a vector) in order to be added to other
matrices. I fiddled around with adding functionality to the
GenericMatric interface for setting diagonals using a GenericVector. The
code can be checked out here:

  bzr branch lp:~johan-hake/dolfin/dolfin-diagonal

I did not push as in the end I did not need it, and I did not bother
adding the last touch to get it ready (read adding unit tests!).

Johan

Sorry, I still dont get it to work, so this is my problem:

Lets assume some mesh and
V = VectorFunctionSpace(mesh=mesh,family="CG",degree=1,dim=3)
u_tst = TestFunction(V)
u_trl = TrialFunction(V)
u = Function(V)

then, the consistent mass matrix is:
Mc = assemble( dot(n_tst,n_trl)*dx )

Mc turns out to be of type <Matrix wrapper of <uBLASMatrix of size ... x ... >>. Trying to use ml = Mc.lump() or Mc.lump(ml), where ml should be the lumped uBLASVector do not work, because of error: 'Matrix' object has no attribute 'lump'.

also, the variant to use the action, i.e.
n_mL = assemble(action(dot(n_tst,n_trl)*dx), coefficients=[Constant(1)])
results in some error (in \dolfin\cpp.py, line 18049)??

What do I do wrong, and how can I fix this?

Next step would be to get the lumped vector back on a Matrix which has only non-zero entries on its diagonal, thus a very spare matrix, completely diagonal. Say this diagonal matrix is called Ml, then I could use solve(Ml,u.vector(),Rhs), or do some other algebra M = Mc-Ml for example.
I had a look at your code Johan, and i can see how to adopt this for my purposes, but unfortunately I do not understand which files I need to change? Is it enough to just change dolfin/la/uBLASMatrix.h, or do I need to also change other files (eg., dolfin/la/Matrix.h)? Once I change this or these files, do I need to recompile the program, or are these changes considered immediately for the next problem-run?
Alternatively, can I modify the entries of Mc directly within my problem-code written in python? How could I iterate through the entries of Mc within python, or is this a bad idea?

Thank you very much for your help!
Arnd

On 08/23/2012 08:20 PM, Arnd Flatten wrote:
> Question #206318 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/206318
>
> Status: Answered => Open
>
> Arnd Flatten is still having a problem:
> Sorry, I still dont get it to work, so this is my problem:
>
> Lets assume some mesh and
> V = VectorFunctionSpace(mesh=mesh,family="CG",degree=1,dim=3)
> u_tst = TestFunction(V)
> u_trl = TrialFunction(V)
> u = Function(V)
>
> then, the consistent mass matrix is:
> Mc = assemble( dot(n_tst,n_trl)*dx )
>
> Mc turns out to be of type <Matrix wrapper of <uBLASMatrix of size ... x
> ... >>. Trying to use ml = Mc.lump() or Mc.lump(ml), where ml should be
> the lumped uBLASVector do not work, because of error: 'Matrix' object
> has no attribute 'lump'.
>
> also, the variant to use the action, i.e.
> n_mL = assemble(action(dot(n_tst,n_trl)*dx), coefficients=[Constant(1)])
> results in some error (in \dolfin\cpp.py, line 18049)??
>
> What do I do wrong, and how can I fix this?
>
>
> Next step would be to get the lumped vector back on a Matrix which has only non-zero entries on its diagonal, thus a very spare matrix, completely diagonal. Say this diagonal matrix is called Ml, then I could use solve(Ml,u.vector(),Rhs), or do some other algebra M = Mc-Ml for example.
> I had a look at your code Johan, and i can see how to adopt this for my purposes, but unfortunately I do not understand which files I need to change? Is it enough to just change dolfin/la/uBLASMatrix.h, or do I need to also change other files (eg., dolfin/la/Matrix.h)? Once I change this or these files, do I need to recompile the program, or are these changes considered immediately for the next problem-run?
> Alternatively, can I modify the entries of Mc directly within my problem-code written in python? How could I iterate through the entries of Mc within python, or is this a bad idea?

I can have a look at the code onces more and see if I can port it to
trunk. Most of the functionality should be there. I just need to check
that it works in parallel too. Might take some days...

And, yes you need to recompile the whole DOLFIN library to get any
changes inside the library to apply.

Johan

> Thank you very much for your help!
> Arnd
>

Thank you for looking further into this!

In the meantime, I came up with the following attempt - it seems to be working OK, but presumably there is a nicer way than looping over each and every row:

def lumpMatrix(sparseMatrix):
 import numpy as np

 assert(sparseMatrix.rank()==2)
 assert(sparseMatrix.size(0)==sparseMatrix.size(1))
 assert(parameters["linear_algebra_backend"]=="uBLAS")

 ### create new uBLASSparseMatrix to be returned
 matSize = sparseMatrix.size(0)
 lumpedMatrix = uBLASSparseMatrix(matSize,matSize)
 ### assign diagonal values of this matrix
 for rowno in range(matSize):
  #colnos,vals = sparseMatrix.getrow(rowno)
  #rowsum = sum(vals)
  #newcolnos = np.array([rowno],dtype='I') # np.uint
  #newvals = np.array([rowsum],dtype='d') # np.float
  #lumpedMatrix.setrow(rowno,newcolnos,newvals)
  rowsum = sum(sparseMatrix.getrow(rowno)[1])
  lumpedMatrix.setrow(rowno, \
                  np.array([rowno],dtype='I'),\
    np.array([rowsum],dtype='d'))
 ### return the lumped matrix
 return lumpedMatrix

Mc = assemble(dot(u_tst,u_trl)*dx)
Ml = lumpMatrix(Mc)

Since lumped mass matrix is just the result of applying a different quadrature rule when assembling the mass matrix, I think a better way to get diagonal lumped mass matrix is to change quadrature rule. For example, let's use "a= inner(u,v)*dx" in "Poisson.ufl". After "ffc -l dolfin Poisson.ufl" we can use the following code in function poisson_cell_integral_0_otherwise::tabulate_tensor to get lumped mass matrix.

    // Compute element tensor
    // A[0] = 0.0833333333333334*G0_;
    // A[1] = 0.0416666666666667*G0_;
    // A[2] = 0.0416666666666667*G0_;
    // A[3] = 0.0416666666666667*G0_;
    // A[4] = 0.0833333333333333*G0_;
    // A[5] = 0.0416666666666666*G0_;
    // A[6] = 0.0416666666666667*G0_;
    // A[7] = 0.0416666666666666*G0_;
    // A[8] = 0.0833333333333333*G0_;

    // new way
    A[0] = 0.5*det;
    A[4] = 0.5*det;
    A[8] = 0.5*det;