[Getfem-commits] (no subject)

[Konstantinos Poulios]

branch: master
commit 7bb6268049999295c04f65cf587428792dc879ec
Author: Konstantinos Poulios <address@hidden>
Date:   Sat Sep 22 17:02:55 2018 +0200

    Final fix in second derivative of Det()
---
 ...fem_generic_assembly_functions_and_operators.cc | 23 +++++++++++-----------
 1 file changed, 11 insertions(+), 12 deletions(-)

diff --git a/src/getfem_generic_assembly_functions_and_operators.cc 
b/src/getfem_generic_assembly_functions_and_operators.cc
index 0c73844..796e2c7 100644
--- a/src/getfem_generic_assembly_functions_and_operators.cc
+++ b/src/getfem_generic_assembly_functions_and_operators.cc
@@ -297,7 +297,8 @@ namespace getfem {
       }
     }
 
-    // Second derivative : det(M)(address@hidden - M^{-T}_{jk}M^{-T}_{li})
+    // Second derivative : det(M)(address@hidden - M^{-T}_{kj}M^{-T}_{il})
+    //                   = det(M)(address@hidden - M^{-1}_{jk}M^{-1}_{li})
     void second_derivative(const arg_list &args, size_type, size_type,
                            base_tensor &result) const {
       size_type N = args[0]->sizes()[0];
@@ -308,23 +309,21 @@ namespace getfem {
         gmm::clear(result.as_vector());
       else {
         auto it = result.begin();
-        auto ita = __mat_aux1().begin(), ita_l = ita, ita_0l = ita;
-        for (size_type l = 0; l < N; ++l, ++ita_l, ita_0l += N) {
-          auto ita_lk = ita_l;
-          for (size_type k = 0; k < N; ++k, ita_lk += N) {
-            auto ita_j = ita, ita_kj = ita + k;
-            for (size_type j = 0; j < N; ++j, ++ita_j, ita_kj += N) {
-              auto ita_ji = ita_j, ita_il = ita_0l;
-              for (size_type i = 0; i < N; ++i, ++it, ita_ji += N, ita_il++)
-                *it = ((*ita_ji) * (*ita_lk) - (*ita_kj) * (*ita_il)) * det;
+        auto ita = __mat_aux1().begin();
+        for (size_type l = 0; l < N; ++l) {
+          auto ita_lk = ita + l, ita_0k = ita;
+          for (size_type k = 0; k < N; ++k, ita_lk += N, ita_0k += N) {
+            auto ita_jk = ita_0k;
+            for (size_type j = 0; j < N; ++j, ++ita_jk) {
+              auto ita_ji = ita + j, ita_li = ita + l;
+              for (size_type i = 0; i < N; ++i, ++it, ita_ji += N, ita_li += N)
+                *it = ((*ita_ji) * (*ita_lk) - (*ita_jk) * (*ita_li)) * det;
             }
           }
         }
         GA_DEBUG_ASSERT(it == result.end(), "Internal error");
       }
     }
-
-    
   };
 
   // Inverse Operator (for square matrices)