[Getfem-commits] (no subject)

[Konstantinos Poulios]

branch: devel-logari81-internal-variables
commit 232385db09a7ff588cb8144fd1e8aa79ce6a32b6
Author: Konstantinos Poulios <address@hidden>
AuthorDate: Wed Jan 29 09:47:03 2020 +0100

    Fix size of internal variable condensation matrix
    
      - Including empty rows for keeping the implementation simpler
---
 src/getfem/getfem_generic_assembly.h            | 4 +++-
 src/getfem_generic_assembly_compile_and_exec.cc | 7 +++++++
 src/getfem_generic_assembly_workspace.cc        | 6 ++----
 3 files changed, 12 insertions(+), 5 deletions(-)

diff --git a/src/getfem/getfem_generic_assembly.h 
b/src/getfem/getfem_generic_assembly.h
index 16b77b6..367d227 100644
--- a/src/getfem/getfem_generic_assembly.h
+++ b/src/getfem/getfem_generic_assembly.h
@@ -413,7 +413,9 @@ namespace getfem {
       KQJpr = std::shared_ptr<model_real_sparse_matrix>
               (std::shared_ptr<model_real_sparse_matrix>(), &KQJpr_); // alias
     }
-    // getter functions for condensation matrix/vectors
+    /** getter function for condensation matrix
+        Its size is (nb_primary_dof()+nb_internal_dof()) x nb_primary_dof()
+        but its first nb_primary_dof() rows are empty */
     const model_real_sparse_matrix &internal_coupling_matrix() const
     { return *KQJpr; }
     model_real_sparse_matrix &internal_coupling_matrix() { return *KQJpr; }
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc 
b/src/getfem_generic_assembly_compile_and_exec.cc
index 5c301ba..11c2a72 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -8000,6 +8000,13 @@ namespace getfem {
                 const base_tensor &Kq1j2pr = *(CC.KQJpr(q1,j2)); // <- input
                 model_real_sparse_matrix
                   &KQJpr = workspace.internal_coupling_matrix(); // <- output
+                GMM_ASSERT1(
+                  gmm::mat_nrows(KQJpr) == workspace.nb_primary_dof()
+                                          +workspace.nb_internal_dof() &&
+                  gmm::mat_ncols(KQJpr) == workspace.nb_primary_dof(),
+                  "The internal coupling matrix needs to be allocated with "
+                  "nb_primary_dof()+nb_internal_dof() number of rows, even if "
+                  "only the last nb_internal_dof() rows are filled in.");
                 if (mf2)
                   pgai =
                     std::make_shared<ga_instruction_matrix_assembly_imd_mf>
diff --git a/src/getfem_generic_assembly_workspace.cc 
b/src/getfem_generic_assembly_workspace.cc
index 50c0db1..4b7517f 100644
--- a/src/getfem_generic_assembly_workspace.cc
+++ b/src/getfem_generic_assembly_workspace.cc
@@ -789,9 +789,9 @@ namespace getfem {
       //              gmm::mat_ncols(*K) == nb_prim_dof, "Wrong sizes");
       if (KQJpr.use_count()) {
         gmm::clear(*KQJpr);
-        gmm::resize(*KQJpr, nb_intern_dof, nb_prim_dof); // redundant if 
condensation == false
+        gmm::resize(*KQJpr, nb_prim_dof+nb_intern_dof, nb_prim_dof); // 
redundant if condensation == false
       } else if (condensation)
-        GMM_ASSERT1(gmm::mat_nrows(*KQJpr) == nb_intern_dof &&
+        GMM_ASSERT1(gmm::mat_nrows(*KQJpr) == nb_prim_dof+nb_intern_dof &&
                     gmm::mat_ncols(*KQJpr) == nb_prim_dof, "Wrong sizes");
       gmm::clear(col_unreduced_K);
       gmm::clear(row_unreduced_K);
@@ -894,8 +894,6 @@ namespace getfem {
                       GMM_ASSERT1(I1.last() <= nb_prim_dof, "Internal error");
                       gmm::add(M, gmm::sub_matrix(*K, I1, I2)); // -> *K
                     } else { // vname1 is an internal variable
-                      I1.min -= first_internal_dof();
-                      I1.max -= first_internal_dof();
                       gmm::add(M, gmm::sub_matrix(*KQJpr, I1, I2)); // -> 
*KQJpr
                     }
                   }