[Getfem-commits] r5486 - in /trunk/getfem/src: ./ getfem/

[logari81]

Author: logari81
Date: Mon Nov 21 23:13:27 2016
New Revision: 5486

URL: http://svn.gna.org/viewcvs/getfem?rev=5486&view=rev
Log:
return compatible Hessian tensor sizes from all implementations of 
real_hess_base_value() (last two dimensions are flattened)

Modified:
    trunk/getfem/src/getfem/getfem_fem_global_function.h
    trunk/getfem/src/getfem_fem.cc
    trunk/getfem/src/getfem_fem_global_function.cc
    trunk/getfem/src/getfem_fem_level_set.cc
    trunk/getfem/src/getfem_generic_assembly.cc
    trunk/getfem/src/getfem_mesh_fem_product.cc
    trunk/getfem/src/getfem_mesh_fem_sum.cc

Modified: trunk/getfem/src/getfem/getfem_fem_global_function.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_fem_global_function.h?rev=5486&r1=5485&r2=5486&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_fem_global_function.h        (original)
+++ trunk/getfem/src/getfem/getfem_fem_global_function.h        Mon Nov 21 
23:13:27 2016
@@ -52,7 +52,6 @@
     const mesh_im &mim;
     const bool has_mesh_im;
 
-    mutable bgeot::multi_index mib,mig,mih;
     mutable std::vector<std::vector<size_type> > index_of_global_dof_;
     mutable bgeot::pstored_point_tab pspt_override;
 

Modified: trunk/getfem/src/getfem_fem.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_fem.cc?rev=5486&r1=5485&r2=5486&view=diff
==============================================================================
--- trunk/getfem/src/getfem_fem.cc      (original)
+++ trunk/getfem/src/getfem_fem.cc      Mon Nov 21 23:13:27 2016
@@ -269,10 +269,7 @@
       }
 
       if (tt.size()) { /* only if the FEM can provide hess_base_value */
-        bgeot::multi_index mim(3);
-        mim[2] = gmm::sqr(tt.sizes()[2]); mim[1] = tt.sizes()[1];
-        mim[0] = tt.sizes()[0];
-        tt.adjust_sizes(mim);
+        tt.adjust_sizes(tt.sizes()[0], tt.sizes()[1], gmm::sqr(tt.sizes()[2]));
         t.mat_transp_reduction(tt, B3(), 2);
         if (!pgt()->is_linear()) {
           if (have_pfp()) {

Modified: trunk/getfem/src/getfem_fem_global_function.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_fem_global_function.cc?rev=5486&r1=5485&r2=5486&view=diff
==============================================================================
--- trunk/getfem/src/getfem_fem_global_function.cc      (original)
+++ trunk/getfem/src/getfem_fem_global_function.cc      Mon Nov 21 23:13:27 2016
@@ -79,16 +79,6 @@
         = std::make_shared<precomp_pool_key>(precomps);
       dal::add_stored_object(pkey, precomps);
     }
-
-    mib.resize(2);
-    mib[0] = short_type(1);
-    mib[1] = target_dim();
-    mig.resize(3);
-    mig[0] = short_type(1);
-    mig[1] = target_dim(); mig[2] = dim();
-    mih.resize(4);
-    mih[0] = short_type(1);
-    mih[1] = target_dim(); mih[3] = mih[2] = dim();
 
     size_type nb_total_dof(functions.size());
     base_node bmin(dim()), bmax(dim());
@@ -193,8 +183,7 @@
     assert(target_dim() == 1);
     size_type cv = c.convex_num();
     size_type nbdof = nb_dof(cv);
-    mib[0] = short_type(nbdof);
-    t.adjust_sizes(mib);
+    t.adjust_sizes(nbdof, target_dim());
     if (c.have_pfp() && c.ii() != size_type(-1)) {
       GMM_ASSERT1(precomps, "Internal error");
       if (precomps->size() == 0)
@@ -235,8 +224,7 @@
     assert(target_dim() == 1);
     size_type cv = c.convex_num();
     size_type nbdof = nb_dof(cv);
-    mig[0] = short_type(nbdof);
-    t.adjust_sizes(mig);
+    t.adjust_sizes(nbdof, target_dim(), dim());
     if (c.have_pfp() && c.ii() != size_type(-1)) {
       GMM_ASSERT1(precomps, "Internal error");
       if (precomps->size() == 0)
@@ -274,8 +262,7 @@
     assert(target_dim() == 1);
     size_type cv = c.convex_num();
     size_type nbdof = nb_dof(cv);
-    mih[0] = short_type(nbdof);
-    t.adjust_sizes(mih);
+    t.adjust_sizes(nbdof, target_dim(), gmm::sqr(dim()));
     if (c.have_pfp() && c.ii() != size_type(-1)) {
       GMM_ASSERT1(precomps, "Internal error");
       if (precomps->size() == 0)

Modified: trunk/getfem/src/getfem_fem_level_set.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_fem_level_set.cc?rev=5486&r1=5485&r2=5486&view=diff
==============================================================================
--- trunk/getfem/src/getfem_fem_level_set.cc    (original)
+++ trunk/getfem/src/getfem_fem_level_set.cc    Mon Nov 21 23:13:27 2016
@@ -201,10 +201,7 @@
   
   void fem_level_set::real_hess_base_value(const fem_interpolation_context &c,
                                            base_tensor &t, bool) const {
-    // bgeot::multi_index mi(4);
-    // mi[3] = mi[2] = short_type(c.N()); mi[1] = target_dim();
-    // mi[0] = short_type(nb_base(0));
-    t.adjust_sizes(nb_base(0), target_dim(), c.N(), c.N());
+    t.adjust_sizes(nb_base(0), target_dim(), gmm::sqr(c.N()));
     fem_interpolation_context c0 = c;
     if (c0.have_pfp())
       c0.set_pfp(fem_precomp(bfem, c0.pfp()->get_ppoint_tab(), c0.pfp()));

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=5486&r1=5485&r2=5486&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Mon Nov 21 23:13:27 2016
@@ -3194,7 +3194,7 @@
 
   struct ga_instruction_hess_base : public ga_instruction_val_base {
 
-    virtual int exec() { // --> t(ndof,target_dim,N,N)
+    virtual int exec() { // --> t(ndof,target_dim,N*N)
       GA_DEBUG_INFO("Instruction: compute Hessian of base functions");
       if (ctx.have_pgp()) ctx.set_pfp(pfp);
       else ctx.set_pf(mf.fem_of_element(ctx.convex_num()));
@@ -3211,7 +3211,7 @@
 
   struct ga_instruction_xfem_plus_hess_base : public ga_instruction_val_base {
 
-    virtual int exec() { // --> t(ndof,target_dim,N,N)
+    virtual int exec() { // --> t(ndof,target_dim,N*N)
       GA_DEBUG_INFO("Instruction: compute Hessian of base functions");
       if (ctx.have_pgp()) ctx.set_pfp(pfp);
       else ctx.set_pf(mf.fem_of_element(ctx.convex_num()));
@@ -3232,7 +3232,7 @@
 
   struct ga_instruction_xfem_minus_hess_base : public ga_instruction_val_base {
 
-    virtual int exec() { // --> t(ndof,target_dim,N,N)
+    virtual int exec() { // --> t(ndof,target_dim,N*N)
       GA_DEBUG_INFO("Instruction: compute Hessian of base functions");
       if (ctx.have_pgp()) ctx.set_pfp(pfp);
       else ctx.set_pf(mf.fem_of_element(ctx.convex_num()));
@@ -3364,55 +3364,50 @@
   };
 
   struct ga_instruction_hess : public ga_instruction_val {
-    // Z(ndof,target_dim,N,N), coeff(Qmult,ndof) --> t(target_dim*Qmult,N,N)
+    // Z(ndof,target_dim,N*N), coeff(Qmult,ndof) --> t(target_dim*Qmult,N,N)
     virtual int exec() {
       GA_DEBUG_INFO("Instruction: Hessian");
       size_type ndof = Z.sizes()[0];
       if (!ndof) { gmm::clear(t.as_vector()); return 0; }
-      size_type N = t.sizes().back();
-      GA_DEBUG_ASSERT(N == Z.sizes()[2], "Internal error");
-      GA_DEBUG_ASSERT(N == Z.sizes()[3], "Internal error");
+      size_type NN = gmm::sqr(t.sizes().back());
+      GA_DEBUG_ASSERT(NN == Z.sizes()[2], "Internal error");
       if (qdim == 1) {
-       GA_DEBUG_ASSERT(gmm::vect_size(coeff) == ndof,
-                       "Wrong size for coeff vector");
-       auto it = Z.begin(); auto itt = t.begin();
-       for (size_type k = 0; k < N; ++k)
-         for (size_type l = 0; l < N; ++l, ++itt) {
-           *itt = scalar_type(0);
-           for (auto itc = coeff.begin(); itc != coeff.end(); ++itc, ++it)
-             *itt += (*itc) * (*it);
-         }
-       GMM_ASSERT1(itt == t.end(),  "dimensions mismatch");
+        GA_DEBUG_ASSERT(gmm::vect_size(coeff) == ndof,
+                        "Wrong size for coeff vector");
+        auto it = Z.begin(); auto itt = t.begin();
+        for (size_type kl = 0; kl < NN; ++kl, ++itt) {
+          *itt = scalar_type(0);
+          for (auto itc = coeff.begin(); itc != coeff.end(); ++itc, ++it)
+            *itt += (*itc) * (*it);
+        }
+        GMM_ASSERT1(itt == t.end(),  "dimensions mismatch");
       } else {
-       size_type target_dim = Z.sizes()[1];
-       if (target_dim == 1) {
-         GA_DEBUG_ASSERT(t.size() == N*N*qdim, "dimensions mismatch");
-         GA_DEBUG_ASSERT(gmm::vect_size(coeff) == ndof*qdim,
-                         "Wrong size for coeff vector");
-         gmm::clear(t.as_vector());
-         for (size_type q = 0; q < qdim; ++q) {
-           base_tensor::const_iterator it = Z.begin();
-           for (size_type k = 0; k < N; ++k)
-             for (size_type l = 0; l < N; ++l)
-               for (size_type j = 0; j < ndof; ++j, ++it)
-                 t[q + k*qdim + l*qdim*N] += coeff[j*qdim+q] * (*it);
-         }
-       } else {
-         size_type Qmult = qdim / target_dim;
-         GA_DEBUG_ASSERT(t.size() == N*N*qdim, "dimensions mismatch");
-         GA_DEBUG_ASSERT(gmm::vect_size(coeff) == ndof*Qmult,
-                         "Wrong size for coeff vector");
-         gmm::clear(t.as_vector());
-         for (size_type q = 0; q < Qmult; ++q) {
-           base_tensor::const_iterator it = Z.begin();
-           for (size_type k = 0; k < N; ++k)
-             for (size_type l = 0; l < N; ++l)
-               for (size_type r = 0; r < target_dim; ++r)
-                 for (size_type j = 0; j < ndof; ++j, ++it)
-                   t[r + q*target_dim + k*qdim + l*qdim*N]
-                     += coeff[j*Qmult+q] * (*it);
-         } 
-       }
+        size_type target_dim = Z.sizes()[1];
+        if (target_dim == 1) {
+          GA_DEBUG_ASSERT(t.size() == NN*qdim, "dimensions mismatch");
+          GA_DEBUG_ASSERT(gmm::vect_size(coeff) == ndof*qdim,
+                          "Wrong size for coeff vector");
+          gmm::clear(t.as_vector());
+          for (size_type q = 0; q < qdim; ++q) {
+            base_tensor::const_iterator it = Z.begin();
+            for (size_type kl = 0; kl < NN; ++kl)
+              for (size_type j = 0; j < ndof; ++j, ++it)
+                t[q + kl*qdim] += coeff[j*qdim+q] * (*it);
+          }
+        } else {
+          size_type Qmult = qdim / target_dim;
+          GA_DEBUG_ASSERT(t.size() == NN*qdim, "dimensions mismatch");
+          GA_DEBUG_ASSERT(gmm::vect_size(coeff) == ndof*Qmult,
+                          "Wrong size for coeff vector");
+          gmm::clear(t.as_vector());
+          for (size_type q = 0; q < Qmult; ++q) {
+            base_tensor::const_iterator it = Z.begin();
+            for (size_type kl = 0; kl < NN; ++kl)
+              for (size_type r = 0; r < target_dim; ++r)
+                for (size_type j = 0; j < ndof; ++j, ++it)
+                  t[r + q*target_dim + kl*qdim] += coeff[j*Qmult+q] * (*it);
+          } 
+        }
       }
       return 0;
     }
@@ -3643,7 +3638,7 @@
   };
 
   struct ga_instruction_copy_hess_base : public ga_instruction_copy_val_base {
-    // Z(ndof,target_dim,N,N) --> t(Qmult*ndof,Qmult*target_dim,N,N)
+    // Z(ndof,target_dim,N*N) --> t(Qmult*ndof,Qmult*target_dim,N,N)
     virtual int exec() {
       GA_DEBUG_INFO("Instruction: Hessian of test functions");
       size_type target_dim = Z.sizes()[1];
@@ -3652,24 +3647,22 @@
         gmm::copy(Z.as_vector(), t.as_vector());
       } else {
         size_type ndof = Z.sizes()[0];
-        size_type N2 = Z.sizes()[2];
         GA_DEBUG_ASSERT(t.size() == Z.size() * Qmult * Qmult,
                       "Wrong size for Hessian vector");
         gmm::clear(t.as_vector());
         base_tensor::const_iterator itZ = Z.begin();
         size_type s = t.sizes()[0], ss = s * Qmult, sss = s+1;
-        size_type ssss=ss*target_dim;
 
         // Performs t(i*Qmult+j, k*Qmult + j, l, m) = Z(i,k,l*N+m)
-        for (size_type l = 0; l < N2; ++l)
-          for (size_type k = 0; k < target_dim; ++k) {
-            base_tensor::iterator it = t.begin() + (ss * k + ssss*l);
-            for (size_type i = 0; i < ndof; ++i, ++itZ, it += Qmult) {
-              base_tensor::iterator it2 = it;
-              *it2 = *itZ;
-              for (size_type j = 1; j < Qmult; ++j) { it2 += sss; *it2 = *itZ; 
}
-            }
+        size_type NNdim = Z.sizes()[2]*target_dim;
+        for (size_type klm = 0; klm < NNdim; ++klm) {
+          base_tensor::iterator it = t.begin() + (ss * klm);
+          for (size_type i = 0; i < ndof; ++i, ++itZ, it += Qmult) {
+            base_tensor::iterator it2 = it;
+            *it2 = *itZ;
+            for (size_type j = 1; j < Qmult; ++j) { it2 += sss; *it2 = *itZ; }
           }
+        }
       }
       return 0;
     }
@@ -3976,7 +3969,7 @@
     virtual int exec() {
       GA_DEBUG_INFO("Instruction: interpolated variable hessian");
       ga_instruction_interpolate::exec();
-      base_matrix v(qdim, ctx.N()*ctx.N());
+      base_matrix v(qdim, ctx.N()*ctx.N()); // To be optimized
       ctx.pf()->interpolation_hess(ctx, coeff, v, dim_type(qdim));
       gmm::copy(v.as_vector(), t.as_vector());
       return 0;
@@ -4087,7 +4080,7 @@
 
   struct ga_instruction_interpolate_hess_base
     : public ga_instruction_copy_hess_base, ga_instruction_interpolate_base {
-    // ctx --> Z(ndof,target_dim,N,N) --> t(Qmult*ndof,Qmult*target_dim,N,N)
+    // ctx --> Z(ndof,target_dim,N*N) --> t(Qmult*ndof,Qmult*target_dim,N,N)
     virtual int exec() {
       GA_DEBUG_INFO("Instruction: interpolated base hessian");
       ga_instruction_interpolate_base::exec();
@@ -4204,7 +4197,7 @@
   struct ga_instruction_elementary_transformation_hess_base
     : public ga_instruction_copy_hess_base,
              ga_instruction_elementary_transformation_base {
-    // Z(ndof,target_dim,N,N) --> t_out(Qmult*ndof,Qmult*target_dim,N,N)
+    // Z(ndof,target_dim,N*N) --> t_out(Qmult*ndof,Qmult*target_dim,N,N)
     virtual int exec() {
       GA_DEBUG_INFO("Instruction: Hessian of test functions with elementary "
                     "transformation");

Modified: trunk/getfem/src/getfem_mesh_fem_product.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_mesh_fem_product.cc?rev=5486&r1=5485&r2=5486&view=diff
==============================================================================
--- trunk/getfem/src/getfem_mesh_fem_product.cc (original)
+++ trunk/getfem/src/getfem_mesh_fem_product.cc Mon Nov 21 23:13:27 2016
@@ -119,11 +119,8 @@
   }
   
   void fem_product::real_hess_base_value(const fem_interpolation_context &c,
-                                 base_tensor &t, bool) const {
-    bgeot::multi_index mi(4);
-    mi[3] = mi[2] = short_type(c.N()); mi[1] = target_dim();
-    mi[0] = short_type(nb_dof(0));
-    t.adjust_sizes(mi);
+                                         base_tensor &t, bool) const {
+    t.adjust_sizes(nb_dof(0), target_dim(), gmm::sqr(c.N()));
     base_tensor::iterator it = t.begin();
     
     fem_interpolation_context c0 = c;

Modified: trunk/getfem/src/getfem_mesh_fem_sum.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_mesh_fem_sum.cc?rev=5486&r1=5485&r2=5486&view=diff
==============================================================================
--- trunk/getfem/src/getfem_mesh_fem_sum.cc     (original)
+++ trunk/getfem/src/getfem_mesh_fem_sum.cc     Mon Nov 21 23:13:27 2016
@@ -215,10 +215,7 @@
   
   void fem_sum::real_hess_base_value(const fem_interpolation_context &c,
                                      base_tensor &t, bool withM) const {
-    bgeot::multi_index mi(4);
-    mi[3] = mi[2] = short_type(c.N()); mi[1] = target_dim();
-    mi[0] = short_type(nb_dof(0));
-    t.adjust_sizes(mi);
+    t.adjust_sizes(nb_dof(0), target_dim(), gmm::sqr(c.N()));
     base_tensor::iterator it = t.begin(), itf;
     
     fem_interpolation_context c0 = c;