[Getfem-commits] [getfem-commits] branch devel-logari81-internal-variables updated: Explicit recording of terms that require post-assembly reduction

[Konstantinos Poulios]

This is an automated email from the git hooks/post-receive script.

logari81 pushed a commit to branch devel-logari81-internal-variables
in repository getfem.

The following commit(s) were added to 
refs/heads/devel-logari81-internal-variables by this push:
     new 4c9a7e9  Explicit recording of terms that require post-assembly 
reduction
4c9a7e9 is described below

commit 4c9a7e954d691be0df8edd731042ffe558ec9c5e
Author: Konstantinos Poulios <address@hidden>
AuthorDate: Sat Feb 1 09:53:23 2020 +0100

    Explicit recording of terms that require post-assembly reduction
---
 .../getfem_generic_assembly_compile_and_exec.h     |   1 +
 src/getfem_generic_assembly_compile_and_exec.cc    |  21 +++-
 src/getfem_generic_assembly_workspace.cc           | 125 ++++++++++-----------
 3 files changed, 79 insertions(+), 68 deletions(-)

diff --git a/src/getfem/getfem_generic_assembly_compile_and_exec.h 
b/src/getfem/getfem_generic_assembly_compile_and_exec.h
index 84d28ea..19c07e9 100644
--- a/src/getfem/getfem_generic_assembly_compile_and_exec.h
+++ b/src/getfem/getfem_generic_assembly_compile_and_exec.h
@@ -93,6 +93,7 @@ namespace getfem {
     bgeot::geotrans_precomp_pool gp_pool;
     fem_precomp_pool fp_pool;
     std::map<gauss_pt_corresp, bgeot::pstored_point_tab> neighbour_corresp;
+    std::set<std::pair<std::string,std::string>> unreduced_terms;
 
     using region_mim_tuple = std::tuple<const mesh_im *, const mesh_region *, 
psecondary_domain>;
     struct region_mim : public region_mim_tuple {
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc 
b/src/getfem_generic_assembly_compile_and_exec.cc
index f033698..0a8f35a 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -7384,6 +7384,7 @@ namespace getfem {
                   ga_instruction_set &gis, size_type order, bool condensation) 
{
     gis.transformations.clear();
     gis.all_instructions.clear();
+    gis.unreduced_terms.clear();
     workspace.clear_temporary_variable_intervals();
 
     std::map<const ga_instruction_set::region_mim, condensation_description>
@@ -7579,6 +7580,9 @@ namespace getfem {
                            (root->tensor(), Vr, Vu, ctx,
                             vgi.I, vgi.mf, vgi.reduced_mf,
                             gis.coeff, gis.nbpt, gis.ipt, interpolate);
+                    for (const std::string &name
+                         : workspace.variable_group(root->name_test1))
+                      gis.unreduced_terms.emplace(name, "");
                   } else {
                     base_vector &V = mf->is_reduced() ? Vu : Vr;
                     const gmm::sub_interval
@@ -7589,6 +7593,8 @@ namespace getfem {
                     pgai = std::make_shared<ga_instruction_vector_assembly_mf>
                            (root->tensor(), V, ctx, I, *mf,
                             gis.coeff, gis.nbpt, gis.ipt, interpolate);
+                    if (mf->is_reduced())
+                      gis.unreduced_terms.emplace(root->name_test1, "");
                   }
                 } else if (imd) {
                   GMM_ASSERT1(root->interpolate_name_test1.size() == 0,
@@ -7843,6 +7849,12 @@ namespace getfem {
                   }
                 } else if (!workspace.is_internal_variable(root->name_test1) &&
                            !workspace.is_internal_variable(root->name_test2)) {
+
+                  if ((mf1 && mf1->is_reduced()) || (mf2 && mf2->is_reduced())
+                      || has_var_group1 || has_var_group2)
+                    gis.unreduced_terms.emplace(root->name_test1,
+                                                root->name_test2);
+
                   auto &Kxu = (mf1 && mf1->is_reduced()) ? Kuu : Kru;
                   auto &Kxr = (mf1 && mf1->is_reduced()) ? Kur : Krr;
                   auto &Kux = (mf2 && mf2->is_reduced()) ? Kuu : Kur;
@@ -8021,13 +8033,15 @@ namespace getfem {
                   "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)
+                if (mf2) {
                   pgai =
                     std::make_shared<ga_instruction_matrix_assembly_imd_mf>
                     (Kq1j2pr, KQJpr, gis.ctx, gis.ctx,
                      I1, imd1, alpha1, I2, *mf2, alpha2, gis.ipt); // 
constructor without gis.coeff
                     // TODO: name_test2 variable group
-                else // for global variable imd2 == 0
+                    if (mf2->is_reduced())
+                      gis.unreduced_terms.emplace(name_test1, name_test2);
+                } else // for global variable imd2 == 0
                   pgai =
                     std::make_shared<ga_instruction_matrix_assembly_imd_imd>
                     (Kq1j2pr, KQJpr, gis.ctx, gis.ctx,
@@ -8098,6 +8112,9 @@ namespace getfem {
                 auto &Krx = (mf2 && mf2->is_reduced()) ? Kru : Krr;
                 auto &Kxx = (mf2 && mf2->is_reduced()) ? Kxu : Kxr;
 
+                if ((mf1 && mf1->is_reduced()) || (mf2 && mf2->is_reduced()))
+                  gis.unreduced_terms.emplace(name_test1, name_test2);
+
                 if (mf1 && mf2) // TODO: name_test1 or name_test2 variable 
group
                   pgai = std::make_shared
                          <ga_instruction_matrix_assembly_mf_mf>
diff --git a/src/getfem_generic_assembly_workspace.cc 
b/src/getfem_generic_assembly_workspace.cc
index 41821b7..221088e 100644
--- a/src/getfem_generic_assembly_workspace.cc
+++ b/src/getfem_generic_assembly_workspace.cc
@@ -830,75 +830,68 @@ namespace getfem {
     if (order > 0) {
       std::set<std::string> vars_vec_done;
       std::set<std::pair<std::string, std::string> > vars_mat_done;
-      for (ga_tree &tree : gis.trees) {
-        if (tree.root) {
-          std::string &name1 = tree.root->name_test1;
-          std::string &name2 = tree.root->name_test2;
-          const std::vector<std::string> vnames1_(1,name1),
-                                         vnames2_(1,name2);
-          const std::vector<std::string>
-            &vnames1 = variable_group_exists(name1) ? variable_group(name1)
-                                                    : vnames1_,
-            &vnames2 = variable_group_exists(name2) ? variable_group(name2)
-                                                    : vnames2_;
-          if (order == 1) {
-            for (const std::string &vname1 : vnames1) {
-              const mesh_fem *mf1 = associated_mf(vname1);
-              if (mf1 && mf1->is_reduced() && vars_vec_done.count(vname1) == 0)
-              {
-                gmm::sub_interval uI1 = temporary_interval_of_variable(vname1),
-                                  I1 = interval_of_variable(vname1);
-                gmm::mult_add(gmm::transposed(mf1->extension_matrix()),
-                              gmm::sub_vector(unreduced_V, uI1),
-                              gmm::sub_vector(*V, I1));
-                vars_vec_done.insert(vname1);
-              }
+      for (const auto &term : gis.unreduced_terms) {
+        const std::string &name1 = term.first;
+        const std::string &name2 = term.second;
+        const std::vector<std::string>
+          vg1_(1,name1), vg2_(1,name2),
+          &vg1 = variable_group_exists(name1) ? variable_group(name1) : vg1_,
+          &vg2 = variable_group_exists(name2) ? variable_group(name2) : vg2_;
+        if (order == 1) {
+          for (const std::string &vname1 : vg1) {
+            const mesh_fem *mf1 = associated_mf(vname1);
+            if (mf1 && mf1->is_reduced() && vars_vec_done.count(vname1) == 0) {
+              gmm::sub_interval uI1 = temporary_interval_of_variable(vname1),
+                                I1 = interval_of_variable(vname1);
+              gmm::mult_add(gmm::transposed(mf1->extension_matrix()),
+                            gmm::sub_vector(unreduced_V, uI1),
+                            gmm::sub_vector(*V, I1));
+              vars_vec_done.insert(vname1);
             }
-          } else {
-            for (const std::string &vname1 : vnames1) {
-              for (const std::string &vname2 : vnames2) {
-                const mesh_fem *mf1 = associated_mf(vname1),
-                               *mf2 = associated_mf(vname2);
-                if (((mf1 && mf1->is_reduced()) || (mf2 && mf2->is_reduced()))
-                    && vars_mat_done.count(std::make_pair(vname1,vname2)) == 0)
-                {
-                  gmm::sub_interval
-                    uI1 = temporary_interval_of_variable(vname1),
-                    uI2 = temporary_interval_of_variable(vname2),
-                    I1 = interval_of_variable(vname1),
-                    I2 = interval_of_variable(vname2);
-                  if (mf1 && mf1->is_reduced() && mf2 && mf2->is_reduced()) {
-                    model_real_sparse_matrix aux(I1.size(), uI2.size());
-                    model_real_row_sparse_matrix M(I1.size(), I2.size());
-                    gmm::mult(gmm::transposed(mf1->extension_matrix()),
-                              gmm::sub_matrix(row_col_unreduced_K, uI1, uI2),
-                              aux);
-                    gmm::mult(aux, mf2->extension_matrix(), M);
-                    gmm::add(M, gmm::sub_matrix(*K, I1, I2));
-                  } else if (mf1 && mf1->is_reduced()) {
-                    if (condensation && vars_vec_done.count(vname1) == 0) {
-                      gmm::mult_add(gmm::transposed(mf1->extension_matrix()),
-                                    gmm::sub_vector(unreduced_V, uI1),
-                                    gmm::sub_vector(*V, I1));
-                      vars_vec_done.insert(vname1);
-                    }
-                    model_real_sparse_matrix M(I1.size(), I2.size());
-                    gmm::mult(gmm::transposed(mf1->extension_matrix()),
-                              gmm::sub_matrix(row_unreduced_K, uI1, I2), M);
-                    gmm::add(M, gmm::sub_matrix(*K, I1, I2));
-                  } else {
-                    model_real_row_sparse_matrix M(I1.size(), I2.size());
-                    gmm::mult(gmm::sub_matrix(col_unreduced_K, I1, uI2),
-                              mf2->extension_matrix(), M);
-                    if (I1.first() < nb_prim_dof) {
-                      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
-                      gmm::add(M, gmm::sub_matrix(*KQJpr, I1, I2)); // -> 
*KQJpr
-                    }
+          }
+        } else {
+          for (const std::string &vname1 : vg1) {
+            for (const std::string &vname2 : vg2) {
+              const mesh_fem *mf1 = associated_mf(vname1),
+                             *mf2 = associated_mf(vname2);
+              if (((mf1 && mf1->is_reduced()) || (mf2 && mf2->is_reduced()))
+                  && vars_mat_done.count(std::make_pair(vname1,vname2)) == 0) {
+                gmm::sub_interval
+                  uI1 = temporary_interval_of_variable(vname1),
+                  uI2 = temporary_interval_of_variable(vname2),
+                  I1 = interval_of_variable(vname1),
+                  I2 = interval_of_variable(vname2);
+                if (mf1 && mf1->is_reduced() && mf2 && mf2->is_reduced()) {
+                  model_real_sparse_matrix aux(I1.size(), uI2.size());
+                  model_real_row_sparse_matrix M(I1.size(), I2.size());
+                  gmm::mult(gmm::transposed(mf1->extension_matrix()),
+                            gmm::sub_matrix(row_col_unreduced_K, uI1, uI2),
+                            aux);
+                  gmm::mult(aux, mf2->extension_matrix(), M);
+                  gmm::add(M, gmm::sub_matrix(*K, I1, I2));
+                } else if (mf1 && mf1->is_reduced()) {
+                  if (condensation && vars_vec_done.count(vname1) == 0) {
+                    gmm::mult_add(gmm::transposed(mf1->extension_matrix()),
+                                  gmm::sub_vector(unreduced_V, uI1),
+                                  gmm::sub_vector(*V, I1));
+                    vars_vec_done.insert(vname1);
+                  }
+                  model_real_sparse_matrix M(I1.size(), I2.size());
+                  gmm::mult(gmm::transposed(mf1->extension_matrix()),
+                            gmm::sub_matrix(row_unreduced_K, uI1, I2), M);
+                  gmm::add(M, gmm::sub_matrix(*K, I1, I2));
+                } else {
+                  model_real_row_sparse_matrix M(I1.size(), I2.size());
+                  gmm::mult(gmm::sub_matrix(col_unreduced_K, I1, uI2),
+                            mf2->extension_matrix(), M);
+                  if (I1.first() < nb_prim_dof) {
+                    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
+                    gmm::add(M, gmm::sub_matrix(*KQJpr, I1, I2)); // -> *KQJpr
                   }
-                  vars_mat_done.insert(std::make_pair(vname1,vname2));
                 }
+                vars_mat_done.insert(std::make_pair(vname1,vname2));
               }
             }
           }