[Getfem-commits] r5467 - in /trunk/getfem: src/getfem/getfem_assembling.h src/getfem_generic_assembly.cc tests/crack.cc

[Yves . Renard]

Author: renard
Date: Sun Nov  6 19:02:25 2016
New Revision: 5467

URL: http://svn.gna.org/viewcvs/getfem?rev=5467&view=rev
Log:
a fix

Modified:
    trunk/getfem/src/getfem/getfem_assembling.h
    trunk/getfem/src/getfem_generic_assembly.cc
    trunk/getfem/tests/crack.cc

Modified: trunk/getfem/src/getfem/getfem_assembling.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_assembling.h?rev=5467&r1=5466&r2=5467&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_assembling.h (original)
+++ trunk/getfem/src/getfem/getfem_assembling.h Sun Nov  6 19:02:25 2016
@@ -557,41 +557,6 @@
  
 
   // -------- Before this : cleaned ----------
-
-  template<typename MAT>
-  void asm_mass_matrix(const MAT &M, const mesh_im &mim,
-                      const mesh_fem &mf_u1,
-                      const mesh_region &rg = mesh_region::all_convexes()) {
-    generic_assembly assem;
-    if (mf_u1.get_qdim() == 1)
-      assem.set("M(#1,#1)+=sym(comp(Base(#1).Base(#1)))");
-    else
-      assem.set("M(#1,#1)+=sym(comp(vBase(#1).vBase(#1))(:,i,:,i));");
-    assem.push_mi(mim);
-    assem.push_mf(mf_u1);
-    assem.push_mat(const_cast<MAT &>(M));
-    assem.assembly(rg);
-  }
-
-  template<typename MAT>
-  void asm_mass_matrix(const MAT &M, const mesh_im &mim, const mesh_fem &mf_u1,
-                      const mesh_fem &mf_u2,
-                      const mesh_region &rg = mesh_region::all_convexes()) {
-    generic_assembly assem;
-    if (mf_u1.get_qdim() == 1 && mf_u2.get_qdim() == 1)
-      assem.set("M(#1,#2)+=comp(Base(#1).Base(#2))");
-    else if (mf_u1.get_qdim() == 1)
-      assem.set("M(#1,#2)+=comp(Base(#1).vBase(#2))(:,:,1);"); // could be i 
in place of 1
-    else if (mf_u2.get_qdim() == 1)
-      assem.set("M(#1,#2)+=comp(vBase(#1).Base(#2))(:,1,:);");
-    else
-      assem.set("M(#1,#2)+=comp(vBase(#1).vBase(#2))(:,i,:,i);");
-    assem.push_mi(mim);
-    assem.push_mf(mf_u1);
-    assem.push_mf(mf_u2);
-    assem.push_mat(const_cast<MAT &>(M));
-    assem.assembly(rg);
-  }
 
   
   /*
@@ -661,7 +626,7 @@
       @ingroup asm
   */
   template<typename MAT>
-  inline void new_asm_mass_matrix
+  inline void asm_mass_matrix
   (const MAT &M, const mesh_im &mim, const mesh_fem &mf1,
    const mesh_region &rg = mesh_region::all_convexes()) {
 
@@ -674,7 +639,7 @@
     gmm::add(workspace.assembled_matrix(), const_cast<MAT &>(M));
   }
 
-  inline void new_asm_mass_matrix
+  inline void asm_mass_matrix
   (model_real_sparse_matrix &M, const mesh_im &mim,
    const mesh_fem &mf1,
    const mesh_region &rg = mesh_region::all_convexes()) {
@@ -693,7 +658,7 @@
    */
 
   template<typename MAT>
-  inline void new_asm_mass_matrix
+  inline void asm_mass_matrix
   (const MAT &M, const mesh_im &mim, const mesh_fem &mf1, const mesh_fem &mf2,
    const mesh_region &rg = mesh_region::all_convexes()) {
     ga_workspace workspace;
@@ -707,7 +672,7 @@
             const_cast<MAT &>(M));
   }
   
-  inline void new_asm_mass_matrix
+  inline void asm_mass_matrix
   (model_real_sparse_matrix &M, const mesh_im &mim,
    const mesh_fem &mf1, const mesh_fem &mf2,
    const mesh_region &rg = mesh_region::all_convexes()) {

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=5467&r1=5466&r2=5467&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Sun Nov  6 19:02:25 2016
@@ -4771,21 +4771,21 @@
     virtual int exec() {
       GA_DEBUG_INFO("Instruction: reduction operation of size " << nn <<
                    " optimized for vectorized tensor");
-      size_type s1 = tc1.size()/nn, s2 = tc2.size()/nn, s1_1 = s1+1;
-      GA_DEBUG_ASSERT(t.size() == s1*s2, "Internal error");
+      size_type s1 = tc1.size()/nn, s2 = tc2.size()/nn, s2_1 = s2+1;
+      GA_DEBUG_ASSERT(t.size() == s2*s1, "Internal error");
       size_type ss1 = s1/nn, ss2 = s2/nn;
 
       // std::fill(t.begin(), t.end(), scalar_type(0)); // Factorized
       auto it2 = tc2.begin();
       for (size_type j = 0; j < ss2; ++j, it2 += nn) {
-       auto it1 = tc1.begin(), it = t.begin() + j*tc1.size();
-       for (size_type i = 0; i < ss1; ++i, it1 += nn, it += nn) {
+       auto it1 = tc1.begin(), it = t.begin() + j*nn;
+       for (size_type i = 0; i < ss1; ++i, it1 += nn, it += s2*nn) {
                  scalar_type a = (*it1) * (*it2);
          auto itt = it;
-         *itt = a; itt += s1_1; *itt = a;
-                 for (size_type k = 2; k < nn; ++k) { itt += s1_1; *itt = a; }
+         *itt = a; itt += s2_1; *itt = a;
+                 for (size_type k = 2; k < nn; ++k) { itt += s2_1; *itt = a; }
                }
-      }
+      }      
       return 0;
     }
     ga_instruction_reduction_opt1(base_tensor &t_, base_tensor &tc1_,
@@ -4919,7 +4919,7 @@
       to_clear = true;
       t.set_sparsity(10, tc1.qdim());
       return std::make_shared<ga_instruction_reduction_opt1>
-       (t.tensor(), tc1.tensor(), tc2.tensor(), n);
+       (t.tensor(), tc1.tensor(), tc2.tensor(), n);
     }
     
     switch(n) {

Modified: trunk/getfem/tests/crack.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/tests/crack.cc?rev=5467&r1=5466&r2=5467&view=diff
==============================================================================
--- trunk/getfem/tests/crack.cc (original)
+++ trunk/getfem/tests/crack.cc Sun Nov  6 19:02:25 2016
@@ -549,7 +549,6 @@
 bool crack_problem::solve(plain_vector &U) {
   size_type N = mesh.dim();
   ls.reinit();  
-  cout << "ls.get_mesh_fem().nb_dof() = " << ls.get_mesh_fem().nb_dof() << 
"\n";
   for (size_type d = 0; d < ls.get_mesh_fem().nb_basic_dof(); ++d) {
     ls.values(0)[d] = ls_function(ls.get_mesh_fem().point_of_basic_dof(d), 
0)[0];
     ls.values(1)[d] = ls_function(ls.get_mesh_fem().point_of_basic_dof(d), 
0)[1];
@@ -815,20 +814,8 @@
 
     gmm::clear(H0);
     getfem::asm_mass_matrix(H0, mim, mf_mortar, mf_u(), 
-                                MORTAR_BOUNDARY_IN);
-    // cout << "H0 = " << H0 << "\n H0 " << endl; getchar();
-
-    // sparse_matrix H8(mf_mortar.nb_dof(), mf_u().nb_dof());
-    // gmm::clear(H8);
-    // getfem::asm_mass_matrix(H8, mim, mf_mortar, mf_u(), MORTAR_BOUNDARY_IN);
-    // // gmm::add(gmm::scaled(H0, -1.0), H8);
-    // // cout << "H0 = " << H0 << endl;
-    // cout << "H8 = " << H8 << endl;
-    
-
-
-
-    gmm::add(gmm::scaled(gmm::sub_matrix(H0, sub_i, sub_j), -1), H);
+                           MORTAR_BOUNDARY_IN);
+    gmm::add(gmm::scaled(gmm::sub_matrix(H0, sub_i, sub_j), -1.0), H);
 
 
     /* because of the discontinuous partition of mf_u(), some levelset