[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
From: |
Tetsuo Koyama |
Subject: |
[Getfem-commits] (no subject) |
Date: |
Sun, 13 Jan 2019 01:08:36 -0500 (EST) |
branch: devel-tetsuo-asm_lumped_mass_matrix
commit ff6df87f1ec714f852a337d8483f7e901eb266b2
Author: Tetsuo Koyama <address@hidden>
Date: Thu Jan 10 16:40:45 2019 +0900
Fix by two loops which can be very inefficient for large matrices
---
src/getfem/getfem_assembling.h | 19 ++++++++-----------
1 file changed, 8 insertions(+), 11 deletions(-)
diff --git a/src/getfem/getfem_assembling.h b/src/getfem/getfem_assembling.h
index 9330a97..d958318 100644
--- a/src/getfem/getfem_assembling.h
+++ b/src/getfem/getfem_assembling.h
@@ -831,23 +831,20 @@ namespace getfem {
* @ingroup asm
*/
template<typename MAT>
- inline void asm_lumped_mass_matrix
+ inline void asm_lumped_mass_matrix_for_first_order
(const MAT &M, const mesh_im &mim, const mesh_fem &mf1,
const mesh_region &rg = mesh_region::all_convexes()) {
asm_mass_matrix(M, mim, mf1, rg);
size_type nbd = gmm::mat_ncols(M), nbr = gmm::mat_nrows(M);
- MAT TEMP_M(nbd, nbr);
GMM_ASSERT1(nbd == nbr, "mass matrix is not square");
- gmm::clear(TEMP_M);
- for (size_type i =0; i < nbr; ++i) {
- for (size_type j =0; j < nbd; ++j) {
- gmm::add(
- gmm::sub_matrix(M, gmm::sub_interval(i, 1), gmm::sub_interval(j, 1)),
- gmm::sub_matrix(TEMP_M, gmm::sub_interval(i, 1),
gmm::sub_interval(i, 1))
- );
- }
+ typedef typename gmm::linalg_traits<MAT>::value_type T;
+ std::vector<T> V(nbd), W(nbr);
+ gmm::fill(V, T(1));
+ gmm::mult(M, V, W);
+ gmm::clear(const_cast<MAT &>(M));
+ for (size_type i =0; i < nbd; ++i) {
+ (const_cast<MAT &>(M))(i, i) = W[i];
}
- gmm::copy(TEMP_M, const_cast<MAT &>(M));
}
/**