getfem-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Getfem-commits] r5379 - in /trunk/getfem: contrib/opt_assembly/ src/ sr


From: Yves . Renard
Subject: [Getfem-commits] r5379 - in /trunk/getfem: contrib/opt_assembly/ src/ src/getfem/ src/gmm/ tests/
Date: Mon, 03 Oct 2016 16:09:23 -0000

Author: renard
Date: Mon Oct  3 18:09:22 2016
New Revision: 5379

URL: http://svn.gna.org/viewcvs/getfem?rev=5379&view=rev
Log:
a small optimization

Modified:
    trunk/getfem/contrib/opt_assembly/opt_assembly.cc
    trunk/getfem/src/getfem/getfem_fem.h
    trunk/getfem/src/getfem/getfem_mesh_fem.h
    trunk/getfem/src/getfem_generic_assembly.cc
    trunk/getfem/src/getfem_mesh_fem.cc
    trunk/getfem/src/gmm/gmm_vector.h
    trunk/getfem/tests/test_assembly.cc

Modified: trunk/getfem/contrib/opt_assembly/opt_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/contrib/opt_assembly/opt_assembly.cc?rev=5379&r1=5378&r2=5379&view=diff
==============================================================================
--- trunk/getfem/contrib/opt_assembly/opt_assembly.cc   (original)
+++ trunk/getfem/contrib/opt_assembly/opt_assembly.cc   Mon Oct  3 18:09:22 2016
@@ -138,8 +138,8 @@
   ch.init(); ch.tic(); old_asm; ch.toc();                              \
   gmm::copy(K, K2);                                                     \
   cout << "Elapsed time for old assembly " << ch.elapsed() << endl;     \
-  gmm::add(gmm::scaled(gmm::sub_matrix(K3,                             \
-                                       I1_, I2_), scalar_type(-1)), K); \
+  gmm::add(gmm::scaled(gmm::sub_matrix(K3,I1_, I2_),                   \
+                       scalar_type(-1)), K);                           \
   scalar_type norm_error = gmm::mat_norminf(K);                         \
   cout << "Error : " << norm_error << endl;
 
@@ -147,12 +147,13 @@
 #define MAT_TEST_2(nbdof1, nbdof2, expr, mim_, I1_, I2_)                \
   workspace.clear_expressions();                                        \
   workspace.add_expression(expr, mim_);                                 \
+  gmm::clear(K3);                                                      \
   ch.init(); ch.tic(); workspace.assembly(2);   ch.toc();               \
   cout << "Elapsed time for new assembly, alternative expression "      \
           << ch.elapsed() << endl;                                      \
   gmm::copy(K2, K);                                                     \
-  gmm::add(gmm::scaled(gmm::sub_matrix(workspace.assembled_matrix(),    \
-                                       I1_, I1_), scalar_type(-1)), K); \
+  gmm::add(gmm::scaled(gmm::sub_matrix(K3, I1_, I1_),                  \
+                      scalar_type(-1)), K);                            \
   norm_error = gmm::mat_norminf(K);                                     \
   cout << "Error : " << norm_error << endl;
 
@@ -422,32 +423,32 @@
   // storage estimate part for the new assembly, global assembly part,
   // ga_exec cost (instructions not executed), J computation, resizing
   // instruction cost.
-  //                        new  | old  | sto  | asse | exec |  J   |resize|
+  //                        new  | old  | sto  | asse | exec |  J   |
   test_new_assembly(2, 400, 1); // ndofu = 321602 ndofp = 160801 ndofchi = 1201
-  // Mass                 : 0.87 | 0.86 | 0.19 | 0.32 | 0.26 | 0.09 | 0.16 |
-  // Laplacian            : 0.47 | 0.85 | 0.10 | 0.16 | 0.19 | 0.08 | 0.03 |
-  // Homogeneous elas     : 0.77 | 1.91 | 0.23 | 0.30 | 0.18 | 0.07 | 0.08 |
-  // Non-homogeneous elast: 0.94 | 2.32 | 0.26 | 0.32 | 0.18 | 0.08 | 0.08 |
+  // Mass                 : 0.79 | 0.86 | 0.19 | 0.32 | 0.26 | 0.09 |
+  // Laplacian            : 0.43 | 0.85 | 0.10 | 0.16 | 0.19 | 0.08 |
+  // Homogeneous elas     : 0.67 | 1.91 | 0.23 | 0.30 | 0.18 | 0.07 |
+  // Non-homogeneous elast: 0.83 | 2.32 | 0.26 | 0.32 | 0.18 | 0.08 |
   test_new_assembly(3, 36, 1);  // ndofu = 151959 ndofp = 50653 ndofchi = 6553
-  // Mass                 : 1.46 | 1.68 | 0.34 | 0.54 | 0.31 | 0.15 | 0.12 |
-  // Laplacian            : 0.90 | 1.51 | 0.10 | 0.17 | 0.24 | 0.14 | 0.06 |
-  // Homogeneous elas     : 1.94 | 4.77 | 0.88 | 0.95 | 0.24 | 0.14 | 0.11 |
-  // Non-homogeneous elast: 2.06 | 6.81 | 0.74 | 0.86 | 0.24 | 0.14 | 0.11 |
+  // Mass                 : 1.36 | 1.68 | 0.34 | 0.54 | 0.31 | 0.15 |
+  // Laplacian            : 0.89 | 1.51 | 0.10 | 0.17 | 0.24 | 0.14 |
+  // Homogeneous elas     : 1.92 | 4.77 | 0.88 | 0.95 | 0.24 | 0.14 |
+  // Non-homogeneous elast: 2.05 | 6.81 | 0.74 | 0.86 | 0.24 | 0.14 |
   test_new_assembly(2, 200, 2); // ndofu = 321602 ndofp = 160801 ndofchi = 1201
-  // Mass                 : 0.49 | 0.45 | 0.14 | 0.22 | 0.07 | 0.03 | 0.01 |
-  // Laplacian            : 0.24 | 0.38 | 0.06 | 0.10 | 0.06 | 0.03 | 0.03 |
-  // Homogeneous elas     : 0.62 | 1.28 | 0.22 | 0.25 | 0.05 | 0.02 | 0.11 |
-  // Non-homogeneous elast: 0.72 | 2.43 | 0.23 | 0.28 | 0.05 | 0.02 | 0.11 |
+  // Mass                 : 0.49 | 0.45 | 0.14 | 0.22 | 0.07 | 0.03 |
+  // Laplacian            : 0.21 | 0.38 | 0.06 | 0.10 | 0.06 | 0.03 |
+  // Homogeneous elas     : 0.53 | 1.28 | 0.22 | 0.25 | 0.05 | 0.02 |
+  // Non-homogeneous elast: 0.64 | 2.43 | 0.23 | 0.28 | 0.05 | 0.02 |
   test_new_assembly(3, 18, 2);  // ndofu = 151959 ndofp = 50653 ndofchi = 6553
-  // Mass                 : 2.00 | 0.90 | 0.27 | 0.63 | 0.05 | 0.02 | 0.27 |
-  // Laplacian            : 0.30 | 0.58 | 0.16 | 0.17 | 0.04 | 0.02 | 0.02 |
-  // Homogeneous elas     : 2.50 | 3.48 | 1.53 | 1.72 | 0.04 | 0.02 | 0.35 |
-  // Non-homogeneous elast: 2.55 | 9.32 | 1.48 | 1.68 | 0.04 | 0.02 | 0.37 |
+  // Mass                 : 1.95 | 0.90 | 0.27 | 0.63 | 0.05 | 0.02 |
+  // Laplacian            : 0.29 | 0.58 | 0.16 | 0.17 | 0.04 | 0.02 |
+  // Homogeneous elas     : 2.48 | 3.48 | 1.53 | 1.72 | 0.04 | 0.02 |
+  // Non-homogeneous elast: 2.55 | 9.32 | 1.48 | 1.68 | 0.04 | 0.02 |
   test_new_assembly(3, 9, 4);   // ndofu = 151959 ndofp = 50653 ndofchi = 6553
-  // Mass                 : 6.64 | 1.33 | 0.65 | 1.77 | 0.01 | .005 | 0.21 |
-  // Laplacian            : 0.78 | 0.79 | 0.32 | 0.43 | 0.01 | .005 | 0.07 |
-  // Homogeneous elas     : 10.2 | 5.52 | 0.90 | 1.69 | 0.01 | .005 | 2.50 |
-  // Non-homogeneous elast: 10.1 | 48.0 | 0.95 | 1.48 | 0.01 | .005 | 2.45 |
+  // Mass                 : 6.64 | 1.33 | 0.65 | 1.77 | 0.01 | .005 |
+  // Laplacian            : 0.78 | 0.79 | 0.32 | 0.43 | 0.01 | .005 |
+  // Homogeneous elas     : 10.2 | 5.52 | 0.90 | 1.69 | 0.01 | .005 |
+  // Non-homogeneous elast: 10.1 | 48.0 | 0.95 | 1.48 | 0.01 | .005 |
 
   // Conclusions :
   // Desactivation of debug test has no sensible effect.

Modified: trunk/getfem/src/getfem/getfem_fem.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_fem.h?rev=5379&r1=5378&r2=5379&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_fem.h        (original)
+++ trunk/getfem/src/getfem/getfem_fem.h        Mon Oct  3 18:09:22 2016
@@ -699,6 +699,7 @@
     /** get the current convex number */
     size_type convex_num() const;
     bool is_convex_num_valid() const;
+    void invalid_convex_num() { convex_num_ = size_type(-1); }
     /** set the current face number */
     void set_face_num(short_type f) { face_num_ = f; }
     /** get the current face number */

Modified: trunk/getfem/src/getfem/getfem_mesh_fem.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_mesh_fem.h?rev=5379&r1=5378&r2=5379&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_mesh_fem.h   (original)
+++ trunk/getfem/src/getfem/getfem_mesh_fem.h   Mon Oct  3 18:09:22 2016
@@ -160,6 +160,7 @@
     EXTENSION_MATRIX E_;
     mutable bgeot::mesh_structure dof_structure;
     mutable bool dof_enumeration_made;
+    mutable bool is_uniform_;
     mutable size_type nb_total_dof;
     pfem auto_add_elt_pf; /* fem for automatic addition                   */
                        /* of element option. (0 = no automatic addition)  */
@@ -275,6 +276,8 @@
 
     /// Return a reference to the underlying mesh.
     const mesh &linked_mesh(void) const { return *linked_mesh_; }
+
+    bool is_uniform(void) const;
 
     /** Set the degree of the fem for automatic addition
      *  of element option. K=-1 disables the automatic addition.

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=5379&r1=5378&r2=5379&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Mon Oct  3 18:09:22 2016
@@ -9515,29 +9515,45 @@
       }
     }
 
+    bool is_uniform = false;
     if (pnode->test_function_type == 1) {
       if (mf1 || mfg1)
         pgai = std::make_shared<ga_instruction_first_ind_tensor>
           (pnode->t, *pctx1, pnode->qdim1, mf1, mfg1);
+      if (mf1 && mf1->is_uniform())
+       { is_uniform = true; pctx1->invalid_convex_num(); }
     } else if (pnode->test_function_type == 2) {
       if (mf2 || mfg2)
         pgai = std::make_shared<ga_instruction_first_ind_tensor>
           (pnode->t, *pctx2, pnode->qdim2, mf2, mfg2);
+      if (mf2 && mf2->is_uniform())
+       { is_uniform = true; pctx2->invalid_convex_num(); }
     } else if (pnode->test_function_type == 3) {
-      if ((mf1 || mfg1) && (mf2 || mfg2))
+      if ((mf1 || mfg1) && (mf2 || mfg2)) {
         pgai = std::make_shared<ga_instruction_two_first_ind_tensor>
           (pnode->t, *pctx1, *pctx2, pnode->qdim1, mf1, mfg1,
            pnode->qdim2, mf2, mfg2);
-      else if (mf1 || mfg1)
+       if (mf1 && mf1->is_uniform() && mf2 && mf2->is_uniform()) {
+         is_uniform=true;
+         pctx1->invalid_convex_num();
+         pctx2->invalid_convex_num();
+       } 
+      } else if (mf1 || mfg1) {
         pgai = std::make_shared<ga_instruction_first_ind_tensor>
           (pnode->t, *pctx1, pnode->qdim1, mf1, mfg1);
-      else if (mf2 || mfg2)
+       if (mf1 && mf1->is_uniform())
+         { is_uniform=true; pctx1->invalid_convex_num(); }
+      } else if (mf2 || mfg2) {
         pgai = std::make_shared<ga_instruction_second_ind_tensor>
           (pnode->t, *pctx2, pnode->qdim2, mf2, mfg2);
-    }
-    // if (pgai) { pgai->exec(); }
-    if (pgai) rmi.instructions.push_back(std::move(pgai));
-
+       if (mf2 && mf2->is_uniform())
+         { is_uniform=true; pctx2->invalid_convex_num(); }
+      }
+    }
+    if (pgai) {
+      if (is_uniform) pgai->exec();
+      else rmi.instructions.push_back(std::move(pgai));
+    }
 
     // Optimization: detect if an equivalent node has already been compiled
     if (rmi.node_list.find(pnode->hash_value) != rmi.node_list.end()) {

Modified: trunk/getfem/src/getfem_mesh_fem.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_mesh_fem.cc?rev=5379&r1=5378&r2=5379&view=diff
==============================================================================
--- trunk/getfem/src/getfem_mesh_fem.cc (original)
+++ trunk/getfem/src/getfem_mesh_fem.cc Mon Oct  3 18:09:22 2016
@@ -290,13 +290,21 @@
     }
   }
 
+  bool mesh_fem::is_uniform() const {
+    context_check(); if (!dof_enumeration_made) enumerate_dof();
+    return is_uniform_;
+  }
+
   /// Enumeration of dofs
   void mesh_fem::enumerate_dof(void) const {
     bgeot::index_node_pair ipt;
+    is_uniform_ = true;
     GMM_ASSERT1(linked_mesh_ != 0, "Uninitialized mesh_fem");
     context_check();
     if (fe_convex.card() == 0)
       { dof_enumeration_made = true; nb_total_dof = 0; return; }
+    pfem first_pf = f_elems[fe_convex.first_true()];
+    if (first_pf->is_on_real_element()) is_uniform_ = false;
 
     // Gives the Cuthill McKee ordering to iterate on elements
     const std::vector<size_type> &cmk = linked_mesh().cuthill_mckee_ordering();
@@ -346,6 +354,7 @@
     for (size_type cv : cmk) { // Loop on elements
       if (!fe_convex.is_in(cv)) continue;
       pfem pf = fem_of_element(cv);
+      if (pf != first_pf) is_uniform_ = false;
       bgeot::pgeometric_trans pgt = linked_mesh().trans_of_convex(cv);
       bgeot::pstored_point_tab pspt = pf->node_tab(cv);
       if (pgt != pgt_old || pspt != pspt_old)
@@ -436,6 +445,7 @@
   void mesh_fem::clear(void) {
     fe_convex.clear();
     dof_enumeration_made = false;
+    is_uniform_ = true;
     touch(); v_num = act_counter();
     dof_structure.clear();
     use_reduction = false;
@@ -446,6 +456,7 @@
   void mesh_fem::init_with_mesh(const mesh &me, dim_type Q) {
     GMM_ASSERT1(linked_mesh_ == 0, "Mesh level set already initialized");
     dof_enumeration_made = false;
+    is_uniform_ = false;
     auto_add_elt_pf = 0;
     auto_add_elt_K = dim_type(-1);
     Qdim = Q;
@@ -467,6 +478,7 @@
     E_ = mf.E_;
     dof_structure = mf.dof_structure;
     dof_enumeration_made = mf.dof_enumeration_made;
+    is_uniform_ = mf.is_uniform_;
     nb_total_dof = mf.nb_total_dof;
     auto_add_elt_pf = mf.auto_add_elt_pf;
     auto_add_elt_K = mf.auto_add_elt_K;
@@ -491,8 +503,12 @@
   mesh_fem::mesh_fem(const mesh &me, dim_type Q)
     { linked_mesh_ = 0; init_with_mesh(me, Q); }
 
-  mesh_fem::mesh_fem(void)
-    { linked_mesh_ = 0; dof_enumeration_made = false; set_qdim(1); }
+  mesh_fem::mesh_fem(void) {
+    linked_mesh_ = 0;
+    dof_enumeration_made = false;
+    is_uniform_ = true;
+    set_qdim(1);
+  }
 
   mesh_fem::~mesh_fem() { }
 
@@ -538,6 +554,8 @@
         } else if (bgeot::casecmp(tmp, "DOF_ENUMERATION") == 0) {
           dal::bit_vector doflst;
           dof_structure.clear(); dof_enumeration_made = false;
+         is_uniform_ = true;
+         size_type nbdof_unif = size_type(-1);
           touch(); v_num = act_counter();
           while (true) {
             bgeot::get_token(ist, tmp);
@@ -550,6 +568,11 @@
             std::vector<size_type> tab;
             if (convex_index().is_in(ic) && tmp.size() &&
                 isdigit(tmp[0]) && tmp2 == ":") {
+             size_type nbd = nb_basic_dof_of_element(ic);
+             if (nbdof_unif == size_type(-1))
+               nbdof_unif = nbd;
+             else if (nbdof_unif != nbd)
+               is_uniform_ = false;
               tab.resize(nb_basic_dof_of_element(ic));
               for (size_type i=0; i < fem_of_element(ic)->nb_dof(ic);
                    i++) {

Modified: trunk/getfem/src/gmm/gmm_vector.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/gmm/gmm_vector.h?rev=5379&r1=5378&r2=5379&view=diff
==============================================================================
--- trunk/getfem/src/gmm/gmm_vector.h   (original)
+++ trunk/getfem/src/gmm/gmm_vector.h   Mon Oct  3 18:09:22 2016
@@ -573,7 +573,7 @@
     }
 
     inline void wa(size_type c, const T &e)
-    { if (e != T(0)) { if (read_access(c)) *(write_access(c)) += e; } }
+    { if (e != T(0)) { *(write_access(c)) += e; } }
 
     inline T r(size_type c) const
     { const T *p = read_access(c); if (p) return *p; else return T(0); }
@@ -750,7 +750,7 @@
       GMM_ASSERT2(c < nbl, "out of range");
       if (e != T(0)) {
        iterator it = this->lower_bound(c);
-       if (it != this->end()) it->second += e;
+       if (it != this->end() && it->first == c) it->second += e;
        else base_type::operator [](c) = e;
       }
     }
@@ -1440,6 +1440,7 @@
       { return data.end(); }
 
     void w(size_type c, const T &e);
+    void wa(size_type c, const T &e);
     T r(size_type c) const {
       GMM_ASSERT2(c < size_, "out of range");
       if (c < shift || c >= shift + data.size()) return T(0);
@@ -1483,11 +1484,35 @@
       shift = c;
     }
     else if (c >= shift + s) {
-      data.resize(c - shift + 1);
-      std::fill(data.begin() + s, data.end(), T(0));
+      data.resize(c - shift + 1, T(0));
+      // std::fill(data.begin() + s, data.end(), T(0));
     }
     data[c - shift] = e;
   }
+
+  template<typename T>  void slvector<T>::wa(size_type c, const T &e) {
+    GMM_ASSERT2(c < size_, "out of range");
+    size_type s = data.size();
+    if (!s) { data.resize(1, e); shift = c; return; }
+    else if (c < shift) {
+      data.resize(s + shift - c); 
+      typename std::vector<T>::iterator it = data.begin(),it2=data.end()-1;
+      typename std::vector<T>::iterator it3 = it2 - shift + c;
+      for (; it3 >= it; --it3, --it2) *it2 = *it3;
+      std::fill(it, it + shift - c, T(0));
+      shift = c;
+      data[c - shift] = e;
+      return;
+    }
+    else if (c >= shift + s) {
+      data.resize(c - shift + 1, T(0));
+      data[c - shift] = e;
+      return;
+      // std::fill(data.begin() + s, data.end(), T(0));
+    }
+    data[c - shift] += e;
+  }
+  
   
   template <typename T> struct linalg_traits<slvector<T> > {
     typedef slvector<T> this_type;

Modified: trunk/getfem/tests/test_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/tests/test_assembly.cc?rev=5379&r1=5378&r2=5379&view=diff
==============================================================================
--- trunk/getfem/tests/test_assembly.cc (original)
+++ trunk/getfem/tests/test_assembly.cc Mon Oct  3 18:09:22 2016
@@ -471,6 +471,7 @@
       gmm::copy(gmm::mat_const_row(M1,i),r1);
       gmm::copy(gmm::mat_const_row(M2,i),r2);    
       cout << "\nrow(" << i+1 << "),\nM1=" << r1 << "\nM2=" << r2 << endl;
+      cout << "mx = " << mx << " d = " << d << endl;
       fail_cnt++;
       GMM_ASSERT1(false, "Failed ! ");
       break;
@@ -658,6 +659,7 @@
     cout << "done " << c << endl;
 
   }
+  cout << "M1 = " << M1 <<  endl;
   if (do_old && do_new) comp_mat(M1,M2);
   }
 }




reply via email to

[Prev in Thread] Current Thread [Next in Thread]