[Getfem-commits] (no subject)

[Liang Jin Lim]

branch: ignore_empty_integration_points
commit 005de207d4436a764ef81f7ef885ac0f3a98ea4e
Author: lj <address@hidden>
Date:   Mon May 14 13:57:51 2018 +0200

    Ignore zero weighted points in the assembly.
---
 src/getfem/getfem_generic_assembly.h            |  3 ++
 src/getfem_generic_assembly_compile_and_exec.cc | 50 +++++++++++++++----------
 src/getfem_generic_assembly_workspace.cc        |  8 ++++
 3 files changed, 42 insertions(+), 19 deletions(-)

diff --git a/src/getfem/getfem_generic_assembly.h 
b/src/getfem/getfem_generic_assembly.h
index 0f7acb7..2bd1478 100644
--- a/src/getfem/getfem_generic_assembly.h
+++ b/src/getfem/getfem_generic_assembly.h
@@ -348,6 +348,7 @@ namespace getfem {
     std::shared_ptr<base_vector> V;
     base_vector unreduced_V;
     base_tensor assemb_t;
+    bool include_empty_int_pts = false;
 
   public:
 
@@ -507,6 +508,8 @@ namespace getfem {
 
     void assembly(size_type order);
 
+    void set_include_empty_int_points(bool include);
+    bool include_empty_int_points() const;
 
     ga_workspace(const getfem::model &md_, bool enable_all_variables = false);
     ga_workspace(bool, const ga_workspace &gaw);
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc 
b/src/getfem_generic_assembly_compile_and_exec.cc
index ae858f8..ae0483a 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -3953,17 +3953,20 @@ namespace getfem {
     bool interpolate;
     virtual int exec() {
       GA_DEBUG_INFO("Instruction: vector term assembly for fem variable");
+      auto empty_weight = abs(coeff) < 1e-15;
       if (ipt == 0 || interpolate) {
+        if (empty_weight) elem.resize(0);
         elem.resize(t.size());
-        auto itt = t.begin(); auto it = elem.begin(), ite = elem.end();
-        size_type nd = ((t.size()) >> 2);
-        for (size_type i = 0; i < nd; ++i) {
-          *it++ = (*itt++) * coeff; *it++ = (*itt++) * coeff;
-          *it++ = (*itt++) * coeff; *it++ = (*itt++) * coeff;
+        if (!empty_weight) {
+          auto itt = t.begin(); auto it = elem.begin(), ite = elem.end();
+          size_type nd = ((t.size()) >> 2);
+          for (size_type i = 0; i < nd; ++i) {
+            *it++ = (*itt++) * coeff; *it++ = (*itt++) * coeff;
+            *it++ = (*itt++) * coeff; *it++ = (*itt++) * coeff;
+          }
+          for (; it != ite;) *it++ = (*itt++) * coeff;
         }
-        for (; it != ite;) *it++ = (*itt++) * coeff;
-        // gmm::copy(gmm::scaled(t.as_vector(), coeff), elem);
-      } else {
+      } else if (!empty_weight) {
         auto itt = t.begin(); auto it = elem.begin(), ite = elem.end();
         size_type nd = ((t.size()) >> 2);
         for (size_type i = 0; i < nd; ++i) {
@@ -4149,18 +4152,21 @@ namespace getfem {
     std::vector<size_type> dofs1, dofs2, dofs1_sort;
     virtual int exec() {
       GA_DEBUG_INFO("Instruction: matrix term assembly");
+      auto empty_weight = abs(coeff < 1e-15);
       if (ipt == 0 || interpolate) {
+        if (empty_weight) elem.resize(0);
         elem.resize(t.size());
-        auto itt = t.begin(); auto it = elem.begin(), ite = elem.end();
-        scalar_type e = coeff*alpha1*alpha2;
-        size_type nd = ((t.size()) >> 2);
-        for (size_type i = 0; i < nd; ++i) {
-          *it++ = (*itt++) * e; *it++ = (*itt++) * e;
-          *it++ = (*itt++) * e; *it++ = (*itt++) * e;
+        if (!empty_weight) {
+          auto itt = t.begin(); auto it = elem.begin(), ite = elem.end();
+          scalar_type e = coeff*alpha1*alpha2;
+          size_type nd = ((t.size()) >> 2);
+          for (size_type i = 0; i < nd; ++i) {
+            *it++ = (*itt++) * e; *it++ = (*itt++) * e;
+            *it++ = (*itt++) * e; *it++ = (*itt++) * e;
+          }
+          for (; it != ite;) *it++ = (*itt++) * e;
         }
-        for (; it != ite;) *it++ = (*itt++) * e;
-        // gmm::copy(gmm::scaled(t.as_vector(), coeff*alpha1*alpha2), elem);
-      } else {
+      } else if (!empty_weight){
         // Faster than a daxpy blas call on my config
         auto itt = t.begin(); auto it = elem.begin(), ite = elem.end();
         scalar_type e = coeff*alpha1*alpha2;
@@ -6801,7 +6807,8 @@ namespace getfem {
                   gmm::clean(gis.Normal, 1e-13);
                 } else gis.Normal.resize(0);
               }
-              gis.coeff = J * pai->coeff(first_ind+gis.ipt);
+              auto ipt_coeff = pai->coeff(first_ind+gis.ipt);
+              gis.coeff = J * ipt_coeff;
               if (first_gp) {
                 for (size_type j = 0; j < gilb.size(); ++j) j+=gilb[j]->exec();
                 first_gp = false;
@@ -6809,7 +6816,12 @@ namespace getfem {
               if (gis.ipt == 0) {
                 for (size_type j = 0; j < gile.size(); ++j) j+=gile[j]->exec();
               }
-              for (size_type j = 0; j < gil.size(); ++j) j+=gil[j]->exec();
+              if (workspace.include_empty_int_points() ||
+                  ipt_coeff > 1e-15 ||
+                  gis.ipt == 0 ||
+                  gis.ipt == gis.nbpt - 1) {
+                for (size_type j = 0; j < gil.size(); ++j) j+=gil[j]->exec();
+              }
               GA_DEBUG_INFO("");
             }
           }
diff --git a/src/getfem_generic_assembly_workspace.cc 
b/src/getfem_generic_assembly_workspace.cc
index ac86f32..a77ebcc 100644
--- a/src/getfem_generic_assembly_workspace.cc
+++ b/src/getfem_generic_assembly_workspace.cc
@@ -829,6 +829,14 @@ namespace getfem {
     }
   }
 
+  void ga_workspace::set_include_empty_int_points(bool include) {
+    include_empty_int_pts = include;
+  }
+
+  bool ga_workspace::include_empty_int_points() const {
+    return include_empty_int_pts;
+  }
+
   void ga_workspace::clear_expressions() { trees.clear(); }
 
   void ga_workspace::print(std::ostream &str) {