Re: [Getfem-users] Lumped mass matrix

[Yves Renard]

Le 08/03/2012 06:14, Roman Putanowicz a écrit :


Dear Roman,

The generic assembly do not allow to define user function to transform tensors. 
You are right that in some situations it could be interesting and may allow 
simpler expressions.

In the present case, yet another alternative is also to compute the result in a 
vector representing the diagonal of the mass matrix with an expression of the 
type 

Mc=sym(comp(Base(#1).Base(#1)));V(#1)+=(Mc(:,i));";

for the scalar case and

Mc=sym(comp(vBase(#1).vBase(#1)));V(#1)+=(Mc(:,i,j,i));";

for the vectorial case. And then copying the resulting vector into the diagonal 
of the matrix.



You are also right that an option to use a lumped matrix in the 
add_basic_d_on_dt brick would be interesting (this mechanism is not completely 
stabilized in Getfem...)


Yves.

> Dear Yves,
>
> Thank you for your remarks.
>> First a remark : as you built your second version, the loop on all the
>> components is very much more expensive than the use of a fourth order
>> tensor in the first version. The best way is the second strategy would
>> be to iterate on the potential nonzero element, i.e. the "neighbhour"
>> degrees of freedom (an exemple exist in the fonction
>> asm_generalized_dirichlet_constraints of getfem_assembling.h, I think).
> Yes, you are right in the above. In fact, mass lumping can be done
> most effectively on the element level. The row summation approach
> or HRZ method (as presented in)
> http://www.colorado.edu/engineering/cas/courses.d/IFEM.d/IFEM.Ch32.d/IFEM.Ch32.pdf
> require short loops over elements of consistent mass matrix and this is
> very local operation. This is equivalent to my first version with the fourth
> order tensor. The only reason I have used such tensor is that I haven't seen
> the way how to manipulate an element matrix __after__ it is calculated
> but __before__ it is assembled into global system. 
>
> The use of a fourth order tensor is the most general but convoluted and 
> inefficient way of expressing a linear transformation on element matrix
> components.
>
> In other words I was seeking a mechanism expressed in the snippet below:
>     generic_assembly assem;
>     assem.set("T = sym(comp(Base(#1).Base(#1)));"
>               "T = UserMap(T);"
>               "M(#1,#1)+= T;");
>
> where UserMap means application of user function passed to assem object
> much in the way a nonlinear calculation routines are passed.
>
> Calculation of lumped matrices would benefit from provision of UserMap
> as above. However, I am not sure where else it could be useful ?
> Maybe in selective/reduced integration techniques? But then, isn't
> NonLinear operation sufficient for this? 
>
>> This is more complex to obtain correct
>> lumped matrices (i.e. non-singular positive matrices and without to much
>> loss of accuracy) for higher order finite element methods.
> Yes, you are again right, and I have found references presenting solution
> tailored to specific cases.
>
> Currently in GetFEM++ the calculation of mass matrix is hidden for instance
> in add_basic_d_on_dt brick. As there are many approaches to mass lumping
> it seems to me, that the way to go would be to enable specification of mass
> matrix on the time derivative brick level.
>
> Regards,
>
> Roman


-- 

  Yves Renard (address@hidden)       tel : (33) 04.72.43.87.08
  Pole de Mathematiques, INSA-Lyon             fax : (33) 04.72.43.85.29
  20, rue Albert Einstein
  69621 Villeurbanne Cedex, FRANCE
  http://math.univ-lyon1.fr/~renard

---------