Dear All,
At one point I needed lumped mass matrix. In GetFEM one can use
so called Lobatto lumping by using GaussLobatto element and
GaussLobato quadratures. Unfortunately this can be done only
for 1D elements or for tensor products of such elements.
Another approach to mass lumping is to sum all row entries
of the consistent mas matrix and put in on diagonal.
I have tried to implement this inside the built-in tensor
manipulation language of getfem::generic_assembly objects.
The result is in attachment.
Trying to use only generic_assembly provided operation
results in the introduction of a fourth order tensor H
which on the second order tensor is to sum element entries
and put them on diagonal.
Well it works, I have tested it on segment and triangular
elements.
There is one thing which bothers me however -- this approach
requires reduction with sparse but potentially very large
tensor H (number of entries is N^4 where N is
number of degrees of freedom. Number of non-zero entries
is only N^2). How this affects performance is yet to be tested.
I see two other approaches :
a) summing the rows to the assembled global matrix.
b) providing additional functions to generic_assembly language
(for instance sum of rows, creation of diagonal matrix,
selection of entries via 0-1 matrix, etc).
Definitely I will look at this issue, but I am writing this
already in case someone of you have dealt with similar problems
or would like to make some comments.
Regards,
Roman
--
Roman Putanowicz, PhD < [email protected] >
Institute for Computational Civil Engng (L-5)
Dept. of Civil Engng, Cracow Univ. of Technology
www.l5.pk.edu.pl, tel. +48 12 628 2569, fax 2034
namespace getfem {
/**
generic lumped mass matrix assembly (on the whole mesh or on the specified
convex set or boundary)
The lumped mass matix is produced by row summation approach.
@ingroup asm
*/
template<typename MAT>
void asm_lumped_mass_matrix(const MAT &M, const mesh_im &mim,
const mesh_fem &mf_u1,
const mesh_region &rg = mesh_region::all_convexes()) {
generic_assembly assem;
int N = mf_u1.nb_dof();
std::stringstream oss;
oss << "h=data$1(" << N << "," << N << "," << N << "," << N << ");";
if (mf_u1.get_qdim() == 1) {
oss << "Mc=sym(comp(Base(#1).Base(#1)));"
<< "M(#1,#1)+=(Mc(i,j).h(:,:,i,j));";
} else {
oss << "Mc=sym(comp(vBase(#1).vBase(#1))(:,i,:,i));"
"M(#1,#1)+=(Mc(i,j).h(:,:,i,j));";
}
assem.set(oss.str());
assem.push_mi(mim);
assem.push_mf(mf_u1);
/* vector v stores 4-th order tensor H which action on consistent
* mass matrix is to produce lumped mass matrix
*/
gmm::rsvector<double> v(N*N*N*N);
size_type N2 = N*N;
size_type N3 = N*N*N;
for (int i=0; i<N; i++) {
size_type is = i + N*i + N2*i;
for (int j=0; j<N; j++) {
v[is + j*N3] = 1.0;
}
}
assem.push_data(v);
assem.push_mat(const_cast<MAT &>(M));
assem.assembly(rg);
}
} /* namespace getfem */
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
