Dear All,
Continuing my previous post:
> I see two other approaches :
> a) summing the rows to the assembled global matrix.
In the attachment you have the second version of the
assembly routine asm_lumped_mass_matrix_v1() where
the lumping is done after assembly.
The only thing I wonder is, if it makes sense to check for
non-zero entries of the matrix and how to set sensibly the
threshold (which here is set arbitrarily to 1e-10).
Regards,
Roman
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);
}
template<typename MAT>
void asm_lumped_mass_matrix_v1(MAT &M, const mesh_im &mim,
const mesh_fem &mf_u1,
const mesh_region &rg = mesh_region::all_convexes()) {
generic_assembly assem;
if (mf_u1.get_qdim() == 1) {
assem.set("M+=sym(comp(Base(#1).Base(#1)));");
} else {
assem.set("M(#1,#1) +=sym(comp(vBase(#1).vBase(#1))(:,i,:,i));");
}
assem.push_mi(mim);
assem.push_mf(mf_u1);
assem.push_mat(M);
assem.assembly(rg);
for (int i=0; i<gmm::mat_nrows(M); i++) {
typename MAT::value_type s = typename MAT::value_type(0);
for (int j=0; j<gmm::mat_nrows(M); j++) {
if (fabs(M(i,j)) >= 1e-10) {
s = s+M(i,j);
M(i,j) = 0.0;
}
}
M(i,i) = s;
}
}
} /* namespace getfem */
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
