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

[Yves . Renard]

Author: renard
Date: Mon Oct 10 17:54:01 2016
New Revision: 5399

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

Modified:
    trunk/getfem/contrib/opt_assembly/opt_assembly.cc
    trunk/getfem/src/getfem/getfem_mesh.h
    trunk/getfem/src/getfem/getfem_mesh_fem.h
    trunk/getfem/src/getfem_generic_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=5399&r1=5398&r2=5399&view=diff
==============================================================================
--- trunk/getfem/contrib/opt_assembly/opt_assembly.cc   (original)
+++ trunk/getfem/contrib/opt_assembly/opt_assembly.cc   Mon Oct 10 17:54:01 2016
@@ -255,7 +255,7 @@
               Iu, getfem::asm_source_term(V, mim, mf_u, mf_u, U));
     
   }
-  
+
   if (false) {
     
     {VEC_TEST_1("Test for Neumann term", ndofu, "u.Test_u",
@@ -428,40 +428,40 @@
   // - old one,
   // - estimate of the storage in sparse matrices part for the new assembly,
   // - global assembly part (assembly instruction),
-  // - ga_exec cost (instructions not executed),
+  // - ga_exec cost (instructions not executed, includes the compilation),
   // - J computation.
   // - Instructions execution except for assembly ones
   //                        new  | old  | sto  | asse | exec |  J   | Ins  |
   test_new_assembly(2, 400, 1); // ndofu = 321602 ndofp = 160801 ndofchi = 1201
-  // Mass (scalar)        : 0.45 | 0.62 | 0.09 | 0.14 | 0.26 | 0.08 | 0.06 |
-  // Mass (vector)        : 0.51 | 0.83 | 0.14 | 0.23 | 0.20 | 0.09 | 0.08 |
-  // Laplacian            : 0.35 | 0.83 | 0.05 | 0.11 | 0.18 | 0.08 | 0.06 |
-  // Homogeneous elas     : 0.53 | 1.88 | 0.14 | 0.21 | 0.18 | 0.07 | 0.14 |
-  // Non-homogeneous elast: 0.68 | 2.26 | 0.15 | 0.23 | 0.18 | 0.08 | 0.27 |
+  // Mass (scalar)        : 0.37 | 0.61 | 0.07 | 0.10 | 0.21 | 0.06 | 0.06 |
+  // Mass (vector)        : 0.46 | 0.82 | 0.12 | 0.22 | 0.15 | 0.05 | 0.09 |
+  // Laplacian            : 0.26 | 0.83 | 0.05 | 0.07 | 0.14 | 0.06 | 0.05 |
+  // Homogeneous elas     : 0.47 | 1.88 | 0.14 | 0.22 | 0.14 | 0.05 | 0.11 |
+  // Non-homogeneous elast: 0.58 | 2.26 | 0.14 | 0.22 | 0.14 | 0.05 | 0.21 |
   test_new_assembly(3, 36, 1);  // ndofu = 151959 ndofp =  50653 ndofchi = 6553
-  // Mass (scalar)        : 0.53 | 0.97 | 0.09 | 0.16 | 0.30 | 0.14 | 0.07 |
-  // Mass (vector)        : 1.00 | 1.68 | 0.20 | 0.40 | 0.26 | 0.15 | 0.34 |
-  // Laplacian            : 0.79 | 1.49 | 0.05 | 0.11 | 0.24 | 0.14 | 0.44 |
-  // Homogeneous elas     : 1.72 | 4.73 | 0.56 | 0.65 | 0.24 | 0.14 | 0.83 |
-  // Non-homogeneous elast: 1.90 | 6.81 | 0.56 | 0.66 | 0.24 | 0.14 | 1.05 |
+  // Mass (scalar)        : 0.44 | 0.77 | 0.08 | 0.12 | 0.24 | 0.11 | 0.08 |
+  // Mass (vector)        : 0.94 | 1.57 | 0.18 | 0.39 | 0.20 | 0.11 | 0.35 |
+  // Laplacian            : 0.42 | 1.38 | 0.06 | 0.09 | 0.19 | 0.11 | 0.14 |
+  // Homogeneous elas     : 1.28 | 4.58 | 0.48 | 0.59 | 0.19 | 0.11 | 0.51 |
+  // Non-homogeneous elast: 1.44 | 6.81 | 0.49 | 0.63 | 0.19 | 0.11 | 0.62 |
   test_new_assembly(2, 200, 2); // ndofu = 321602 ndofp = 160801 ndofchi = 1201
-  // Mass (scalar)        : 0.17 | 0.25 | 0.04 | 0.07 | 0.07 | 0.03 | 0.04 |
-  // Mass (vector)        : 0.35 | 0.45 | 0.08 | 0.15 | 0.06 | 0.03 | 0.14 |
-  // Laplacian            : 0.16 | 0.38 | 0.03 | 0.05 | 0.06 | 0.03 | 0.05 |
-  // Homogeneous elas     : 0.39 | 1.28 | 0.15 | 0.19 | 0.05 | 0.02 | 0.16 |
-  // Non-homogeneous elast: 0.49 | 2.40 | 0.14 | 0.19 | 0.05 | 0.02 | 0.25 |
+  // Mass (scalar)        : 0.15 | 0.25 | 0.04 | 0.07 | 0.05 | 0.02 | 0.03 |
+  // Mass (vector)        : 0.34 | 0.45 | 0.07 | 0.15 | 0.05 | 0.02 | 0.14 |
+  // Laplacian            : 0.13 | 0.37 | 0.03 | 0.05 | 0.04 | 0.02 | 0.04 |
+  // Homogeneous elas     : 0.37 | 1.28 | 0.13 | 0.19 | 0.04 | 0.01 | 0.14 |
+  // Non-homogeneous elast: 0.46 | 2.40 | 0.13 | 0.18 | 0.04 | 0.01 | 0.24 |
   test_new_assembly(3, 18, 2);  // ndofu = 151959 ndofp =  50653 ndofchi = 6553
-  // Mass (scalar)        : 0.25 | 0.30 | 0.11 | 0.16 | 0.05 | 0.02 | 0.04 |
-  // Mass (vector)        : 1.48 | 0.90 | 0.31 | 0.62 | 0.05 | 0.02 | 0.81 |
-  // Laplacian            : 0.22 | 0.57 | 0.09 | 0.11 | 0.04 | 0.02 | 0.09 |
-  // Homogeneous elas     : 2.43 | 3.48 | 1.19 | 1.37 | 0.03 | 0.02 | 1.03 |
-  // Non-homogeneous elast: 2.52 | 9.25 | 1.21 | 1.37 | 0.03 | 0.02 | 1.12 |
+  // Mass (scalar)        : 0.23 | 0.30 | 0.10 | 0.15 | 0.04 | 0.02 | 0.04 |
+  // Mass (vector)        : 1.46 | 0.90 | 0.33 | 0.64 | 0.04 | 0.02 | 0.78 |
+  // Laplacian            : 0.18 | 0.55 | 0.08 | 0.10 | 0.03 | 0.02 | 0.05 |
+  // Homogeneous elas     : 2.37 | 3.47 | 1.20 | 1.37 | 0.03 | 0.02 | 0.97 |
+  // Non-homogeneous elast: 2.46 | 9.25 | 1.20 | 1.37 | 0.03 | 0.02 | 1.06 |
   test_new_assembly(3, 9, 4);   // ndofu = 151959 ndofp =  50653 ndofchi = 6553
-  // Mass (scalar)        : 0.75 | 0.38 | 0.31 | 0.39 | 0.01 | .005 | 0.35 |
-  // Mass (vector)        : 5.12 | 1.33 | 0.56 | 1.85 | 0.01 | .005 | 3.31 |
-  // Laplacian            : 0.62 | 0.79 | 0.25 | 0.34 | 0.01 | .005 | 0.27 |
-  // Homogeneous elas     : 10.9 | 5.52 | 2.23 | 3.14 | 0.01 | .005 | 7.75 |
-  // Non-homogeneous elast: 11.0 | 48.0 | 1.97 | 3.17 | 0.01 | .005 | 7.82 |
+  // Mass (scalar)        : 0.75 | 0.38 | 0.25 | 0.39 | 0.01 | .005 | 0.35 |
+  // Mass (vector)        : 5.10 | 1.33 | 0.42 | 1.86 | 0.01 | .005 | 3.23 |
+  // Laplacian            : 0.61 | 0.79 | 0.25 | 0.34 | 0.01 | .005 | 0.26 |
+  // Homogeneous elas     : 10.8 | 5.52 | 2.40 | 3.29 | 0.01 | .005 | 7.50 |
+  // Non-homogeneous elast: 11.0 | 48.0 | 2.35 | 3.19 | 0.01 | .005 | 7.80 |
 
   // Conclusions :
   // - Desactivation of debug test has no sensible effect.
@@ -517,8 +517,8 @@
   // Mass (scalar)        : 0.78 | 0.38 |
   // Mass (vector)        : 6.90 | 1.60 | 0.65 | 1.77 | 0.01 | .005 | 0.21 |
   // Laplacian            : 0.90 | 0.89 | 0.32 | 0.43 | 0.01 | .005 | 0.07 |
-  // Homogeneous elas     : 13.1 | 6.65 | 0.90 | 1.69 | 0.01 | .005 | 2.50 |
-  // Non-homogeneous elast: 13.0 | 49.1 | 0.95 | 1.48 | 0.01 | .005 | 2.45 |
+  // Homogeneous elas     : 12.1 | 6.65 | 0.90 | 1.69 | 0.01 | .005 | 2.50 |
+  // Non-homogeneous elast: 12.0 | 49.1 | 0.95 | 1.48 | 0.01 | .005 | 2.45 |
 #endif
 
   return 0; 

Modified: trunk/getfem/src/getfem/getfem_mesh.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_mesh.h?rev=5399&r1=5398&r2=5399&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_mesh.h       (original)
+++ trunk/getfem/src/getfem/getfem_mesh.h       Mon Oct 10 17:54:01 2016
@@ -684,6 +684,20 @@
 
 
   ///@}
+
+  inline void vectors_to_base_matrix(base_matrix &G, const mesh &m,
+                                    size_type cv) {
+    const bgeot::mesh_structure::ind_cv_ct &ct = m.ind_points_of_convex(cv);
+    size_type N = m.dim(), Np = ct.size();
+    G.base_resize(N, Np);
+    auto it = G.begin();
+    for (size_type i = 0; i < Np; ++i, it += N) {
+      const base_node &P = m.points()[ct[i]];
+      std::copy(P.begin(),  P.end(), it);
+    }
+  }
+
+
 }  /* end of namespace getfem.                                             */
 
 

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=5399&r1=5398&r2=5399&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_mesh_fem.h   (original)
+++ trunk/getfem/src/getfem/getfem_mesh_fem.h   Mon Oct 10 17:54:01 2016
@@ -631,7 +631,7 @@
   /** Given a mesh_fem @param mf and a vector @param vec of size equal to
    *  mf.nb_basic_dof(), the output vector @param coeff will contain the
    *  values of @param vec corresponding to the basic dofs of element
-   *  @param cv . The size of @param coeff is adjusted if necessary.
+   *  @param cv. The size of @param coeff is adjusted if necessary.
    */
   template <typename VEC1, typename VEC2>
   void slice_vector_on_basic_dof_of_element(const mesh_fem &mf,
@@ -640,17 +640,22 @@
     size_type nbdof = mf.nb_basic_dof();
     size_type qmult = gmm::vect_size(vec) / nbdof;
     GMM_ASSERT1(gmm::vect_size(vec) == qmult * nbdof, "Bad dof vector size");
-    size_type cvnbdof = mf.nb_basic_dof_of_element(cv);
-    gmm::resize(coeff, cvnbdof*qmult);
-    mesh_fem::ind_dof_ct::const_iterator
-      itdof = mf.ind_basic_dof_of_element(cv).begin();
-    if (qmult == 1) {
-      for (size_type k = 0; k < cvnbdof; ++k, ++itdof)
-          coeff[k] = vec[*itdof];
+
+    auto &ct = mf.ind_scalar_basic_dof_of_element(cv);
+    size_type qmult2 = mf.get_qdim();
+    if (qmult2 > 1) qmult2 /= mf.fem_of_element(cv)->target_dim();
+    size_type qmultot = qmult*qmult2;
+    gmm::resize(coeff, ct.size()*qmultot);
+ 
+    auto it = ct.begin();
+    auto itc = coeff.begin();
+    if (qmultot == 1) {
+      for (; it != ct.end(); ++it) *itc++ = vec[*it];
     } else {
-      for (size_type k = 0; k < cvnbdof; ++k, ++itdof)
-        for (size_type l = 0; l < qmult; ++l)
-          coeff[l+k*qmult] = vec[l+(*itdof)*qmult];
+      for (; it != ct.end(); ++it) {
+       auto itv = vec.begin()+(*it)*qmult;
+       for (size_type m = 0; m < qmultot; ++m) *itc++ = *itv++;
+      }
     }
   }
 

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=5399&r1=5398&r2=5399&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Mon Oct 10 17:54:01 2016
@@ -5206,7 +5206,7 @@
                                       inin.derivatives, compute_der);
       if (inin.pt_type) {
         if (cv != size_type(-1)) {
-          bgeot::vectors_to_base_matrix(inin.G, 
(inin.m)->points_of_convex(cv));
+          vectors_to_base_matrix(inin.G, *(inin.m), cv);
           inin.ctx.change((inin.m)->trans_of_convex(cv),
                          0, P_ref, inin.G, cv, face_num);
           inin.has_ctx = true;
@@ -5301,7 +5301,7 @@
               const bgeot::stored_point_tab
                 &spt = *(pai->pintegration_points());
               base_matrix G;
-              bgeot::vectors_to_base_matrix(G, m.points_of_convex(cv));
+              vectors_to_base_matrix(G, m, cv);
               fem_interpolation_context ctx_x(gpc.pgt1, 0, spt[0], G, cv, f);
               std::vector<base_node> P_ref(nbpt);
 
@@ -5315,8 +5315,7 @@
               pspt = store_point_tab(P_ref);
               neighbour_corresp[gpc] = pspt;
             }
-            bgeot::vectors_to_base_matrix(inin.G,
-                                          m.points_of_convex(adj_face.cv));
+            vectors_to_base_matrix(inin.G, m, adj_face.cv);
             bgeot::pgeotrans_precomp pgp = gp_pool(gpc.pgt2, pspt);
             inin.ctx.change(pgp, 0, 0, inin.G, adj_face.cv, adj_face.f);
           }
@@ -5341,8 +5340,7 @@
                                         inin.derivatives, false);
         if (inin.pt_type) {
           if (cv != size_type(-1)) {
-            bgeot::vectors_to_base_matrix(inin.G,
-                                          (inin.m)->points_of_convex(cv));
+            vectors_to_base_matrix(inin.G, *(inin.m), cv);
             inin.ctx.change((inin.m)->trans_of_convex(cv),
                            0, P_ref, inin.G, cv, face_num);
             inin.has_ctx = true;
@@ -5422,8 +5420,8 @@
        GA_DEBUG_ASSERT(V.size() >= I.first() + mf.nb_basic_dof(),
                        "Bad assembly vector size");
        auto &ct = mf.ind_scalar_basic_dof_of_element(cv_1);
-       size_type qmult = size_type(mf.get_qdim()
-                                   / mf.fem_of_element(cv_1)->target_dim());
+       size_type qmult = mf.get_qdim();
+       if (qmult > 1) qmult /= mf.fem_of_element(cv_1)->target_dim();
        size_type ifirst = I.first();
        auto ite = elem.begin();
        for (auto itc = ct.begin(); itc != ct.end(); ++itc)
@@ -5503,8 +5501,8 @@
          if (!(ctx1.is_convex_num_valid())) return 0;
          size_type cv1 = ctx1.convex_num();
          auto &ct1 = pmf1->ind_scalar_basic_dof_of_element(cv1);
-         size_type qmult1 = size_type(pmf1->get_qdim()
-                                  / pmf1->fem_of_element(cv1)->target_dim());
+         size_type qmult1 = pmf1->get_qdim();
+         if (qmult1 > 1) qmult1 /= pmf1->fem_of_element(cv1)->target_dim();
          auto itd = dofs1.begin();
          if (qmult1 == 1) {
            for (auto itt = ct1.begin(); itt != ct1.end(); ++itt)
@@ -5522,8 +5520,8 @@
          if (!(ctx2.is_convex_num_valid())) return 0;
          size_type cv2 = ctx2.convex_num();
          auto &ct2 = pmf2->ind_scalar_basic_dof_of_element(cv2);
-         size_type qmult2 = size_type(pmf2->get_qdim()
-                                  / pmf2->fem_of_element(cv2)->target_dim());
+         size_type qmult2 = pmf2->get_qdim();
+         if (qmult2 > 1) qmult2 /= pmf2->fem_of_element(cv2)->target_dim();
          auto itd = dofs2.begin();
          if (qmult2 == 1) {
            for (auto itt = ct2.begin(); itt != ct2.end(); ++itt)
@@ -5541,7 +5539,7 @@
         for (const size_type &dof2 : dofs2)
           for (const size_type &dof1 : dofs1) {
             if (gmm::abs(*it) > threshold)
-              K(dof1, dof2) += *it;
+             K(dof1, dof2) += *it;
             ++it;
           }
       }
@@ -5607,7 +5605,7 @@
         for (const size_type &dof2 : ct2)
           for (const size_type &dof1 : ct1) {
             if (gmm::abs(*it) > threshold)
-              K(dof1+i1, dof2+i2) += *it;
+             K(dof1+i1, dof2+i2) += *it;
             ++it;
           }
       }
@@ -5658,20 +5656,19 @@
         size_type cv1 = ctx1.convex_num();
        if (cv1 == size_type(-1)) return 0;
        auto &ct1 = pmf1->ind_scalar_basic_dof_of_element(cv1);
-       size_type qmult1 = size_type(pmf1->get_qdim()
-                                    / pmf1->fem_of_element(cv1)->target_dim());
+       size_type qmult1 = pmf1->get_qdim();
+       if (qmult1 > 1) qmult1 /= pmf1->fem_of_element(cv1)->target_dim();
        dofs1.assign(s1, I1.first());
        auto itd = dofs1.begin();
        for (auto itt = ct1.begin(); itt != ct1.end(); ++itt)
          for (size_type q = 0; q < qmult1; ++q)
            *itd++ += *itt + q;
-
        
        size_type cv2 = ctx2.convex_num();
        if (cv2 == size_type(-1)) return 0;
-       auto &ct2 = pmf2->ind_scalar_basic_dof_of_element(cv2); 
-       size_type qmult2 = size_type(pmf2->get_qdim()
-                                    / pmf2->fem_of_element(cv2)->target_dim());
+       auto &ct2 = pmf2->ind_scalar_basic_dof_of_element(cv2);
+       size_type qmult2 = pmf2->get_qdim();
+       if (qmult2 > 1) qmult2 /= pmf2->fem_of_element(cv2)->target_dim();
        dofs2.assign(s2, I2.first());
        itd = dofs2.begin();
        for (auto itt = ct2.begin(); itt != ct2.end(); ++itt)
@@ -6484,7 +6481,6 @@
 
 
   void ga_workspace::assembly(size_type order) {
-
     GA_TIC;
     ga_instruction_set gis;
     ga_compile(*this, gis, order);
@@ -11646,7 +11642,7 @@
           = gic.ppoints_for_element(v.cv(), v.f(), ind);
 
         if (pspt.get() && ind.size() && pspt->size()) {
-          bgeot::vectors_to_base_matrix(G, m.points_of_convex(v.cv()));
+          vectors_to_base_matrix(G, m, v.cv());
           bgeot::pgeometric_trans pgt = m.trans_of_convex(v.cv());
           up.resize(G.nrows());
           un.resize(pgt->dim());
@@ -11751,15 +11747,7 @@
         if (mim.convex_index().is_in(v.cv())) {
           // cout << "proceed with element " << v.cv() << endl;
           if (v.cv() != old_cv) {
-
-           // const bgeot::mesh_structure::ind_cv_ct 
&ct=m.ind_points_of_convex(v.cv());
-           // G.base_resize(m.dim(), ct.size());
-           // auto it = G.begin();
-           // for (size_type i = 0; i < ct.size(); ++i, it += m.dim())
-           //   std::copy(m.points()[ct[i]].begin(),  m.points()[ct[i]].end(), 
it);
-
-           
-           bgeot::vectors_to_base_matrix(G, m.points_of_convex(v.cv()));
+           vectors_to_base_matrix(G, m, v.cv());
             pgt = m.trans_of_convex(v.cv());
             up.resize(G.nrows());
             un.resize(pgt->dim());