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 */