[Getfem-commits] r5443 - in /trunk/getfem: contrib/opt_assembly/opt_assembly.cc src/bgeot_geometric_trans.cc

[Yves . Renard]

Author: renard
Date: Wed Oct 26 11:00:19 2016
New Revision: 5443

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

Modified:
    trunk/getfem/contrib/opt_assembly/opt_assembly.cc
    trunk/getfem/src/bgeot_geometric_trans.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=5443&r1=5442&r2=5443&view=diff
==============================================================================
--- trunk/getfem/contrib/opt_assembly/opt_assembly.cc   (original)
+++ trunk/getfem/contrib/opt_assembly/opt_assembly.cc   Wed Oct 26 11:00:19 2016
@@ -454,7 +454,7 @@
   if (all || only_one == 1) // ndofu = 321602 ndofp = 160801 ndofchi = 1201
     test_new_assembly(2, 400, 1);
   // Vector source term   : 0.25 | 0.68 |
-  // Nonlinear residual   : 0.35 |      |
+  // Nonlinear residual   : 0.34 |      |
   // Mass (scalar)        : 0.18 | 0.58 | 0.04 | 0.06 | 0.06 | 0.06 |
   // Mass (vector)        : 0.30 | 0.82 | 0.09 | 0.15 | 0.06 | 0.09 |
   // Laplacian            : 0.16 | 0.80 | 0.04 | 0.05 | 0.06 | 0.05 |
@@ -466,9 +466,9 @@
   // Nonlinear residual   : 0.59 |      |
   // Mass (scalar)        : 0.21 | 0.58 | 0.05 | 0.08 | 0.08 | 0.06 |
   // Mass (vector)        : 0.69 | 1.37 | 0.17 | 0.27 | 0.08 | 0.34 |
-  // Laplacian            : 0.18 | 1.37 | 0.03 | 0.07 | 0.08 | 0.03 |
-  // Homogeneous elas     : 0.84 | 4.58 | 0.26 | 0.33 | 0.08 | 0.43 |
-  // Non-homogeneous elast: 0.90 | 6.55 | 0.26 | 0.33 | 0.08 | 0.49 |
+  // Laplacian            : 0.18 | 1.15 | 0.03 | 0.07 | 0.08 | 0.03 |
+  // Homogeneous elas     : 0.79 | 4.29 | 0.26 | 0.33 | 0.08 | 0.38 |
+  // Non-homogeneous elast: 0.87 | 6.35 | 0.26 | 0.33 | 0.08 | 0.46 |
   if (all || only_one == 3) // ndofu = 321602 ndofp = 160801 ndofchi = 1201
     test_new_assembly(2, 200, 2);
   // Vector source term   : 0.11 | 0.24 |
@@ -481,11 +481,11 @@
   if (all || only_one == 4) // ndofu = 151959 ndofp =  50653 ndofchi = 6553
     test_new_assembly(3, 18, 2);
   // Vector source term   : 0.16 | 0.24 |
-  // Nonlinear residual   : 0.40 |      |
+  // Nonlinear residual   : 0.39 |      |
   // Mass (scalar)        : 0.11 | 0.25 | 0.05 | 0.05 | 0.03 | 0.03 |
-  // Mass (vector)        : 1.14 | 0.90 | 0.21 | 0.35 | 0.03 | 0.76 |
-  // Laplacian            : 0.10 | 0.55 | 0.03 | 0.04 | 0.03 | 0.03 |
-  // Homogeneous elas     : 1.66 | 3.41 | 0.59 | 0.73 | 0.03 | 0.90 |
+  // Mass (vector)        : 1.14 | 0.89 | 0.21 | 0.35 | 0.03 | 0.76 |
+  // Laplacian            : 0.10 | 0.53 | 0.03 | 0.04 | 0.03 | 0.03 |
+  // Homogeneous elas     : 1.66 | 3.35 | 0.59 | 0.73 | 0.03 | 0.90 |
   // Non-homogeneous elast: 1.72 | 9.15 | 0.59 | 0.73 | 0.03 | 0.96 |
   if (all || only_one == 5) // ndofu = 151959 ndofp =  50653 ndofchi = 6553
     test_new_assembly(3, 9, 4);
@@ -494,7 +494,7 @@
   // Mass (scalar)        : 0.51 | 0.34 | 0.09 | 0.16 | 0.01 | 0.33 |
   // Mass (vector)        : 4.37 | 1.31 | 0.41 | 1.27 | 0.01 | 3.10 |
   // Laplacian            : 0.40 | 0.77 | 0.09 | 0.14 | 0.01 | 0.25 |
-  // Homogeneous elas     : 8.93 | 5.24 | 0.88 | 1.43 | 0.01 | 7.49 |
+  // Homogeneous elas     : 8.93 | 5.23 | 0.88 | 1.43 | 0.01 | 7.49 |
   // Non-homogeneous elast: 9.01 | 47.4 | 0.76 | 1.35 | 0.01 | 7.65 |
 
   // Conclusions :

Modified: trunk/getfem/src/bgeot_geometric_trans.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/bgeot_geometric_trans.cc?rev=5443&r1=5442&r2=5443&view=diff
==============================================================================
--- trunk/getfem/src/bgeot_geometric_trans.cc   (original)
+++ trunk/getfem/src/bgeot_geometric_trans.cc   Wed Oct 26 11:00:19 2016
@@ -153,16 +153,16 @@
       *A++ =  d/det;  *A++ /= -det; *A++ /= -det;  *A =  a/det;
       return det;
     } else if (N == 3) {
-      scalar_type a0 = A[4]*A[8] - A[5]*A[7], a3 = A[5]*A[6] - A[3]*A[8];
-      scalar_type a6 = A[3]*A[7] - A[4]*A[6];
-      scalar_type det =  A[0] * a0 + A[1] * a3 + A[2] * a6;
+      scalar_type a0 = A[4]*A[8] - A[5]*A[7], a1 = A[5]*A[6] - A[3]*A[8];
+      scalar_type a2 = A[3]*A[7] - A[4]*A[6];
+      scalar_type det =  A[0] * a0 + A[1] * a1 + A[2] * a2;
       GMM_ASSERT1(det != scalar_type(0), "Non invertible matrix");
-      scalar_type a1 = (A[2]*A[7] - A[1]*A[8]), a2 = (A[1]*A[5] - A[2]*A[4]);
-      scalar_type a4 = (A[0]*A[8] - A[2]*A[6]), a5 = (A[2]*A[3] - A[0]*A[5]);
-      scalar_type a7 = (A[1]*A[6] - A[0]*A[7]), a8 = (A[0]*A[4] - A[1]*A[3]);
-      *A++ = a0 / det; *A++ = a1 / det; *A++ = a2 / det;
-      *A++ = a3 / det; *A++ = a4 / det; *A++ = a5 / det;
-      *A++ = a6 / det; *A++ = a7 / det; *A++ = a8 / det;
+      scalar_type a3 = (A[2]*A[7] - A[1]*A[8]), a6 = (A[1]*A[5] - A[2]*A[4]);
+      scalar_type a4 = (A[0]*A[8] - A[2]*A[6]), a7 = (A[2]*A[3] - A[0]*A[5]);
+      scalar_type a5 = (A[1]*A[6] - A[0]*A[7]), a8 = (A[0]*A[4] - A[1]*A[3]);
+      *A++ = a0 / det; *A++ = a3 / det; *A++ = a6 / det;
+      *A++ = a1 / det; *A++ = a4 / det; *A++ = a7 / det;
+      *A++ = a2 / det; *A++ = a5 / det; *A++ = a8 / det;
       return det;
     } else {
       size_type NN = N*N;
@@ -235,10 +235,10 @@
       } else if (P == 3) {
        B_.base_resize(P, P); // co-factors
        auto it = KK.begin(); auto itB = B_.begin();
-       itB[0] = it[4]*it[8] - it[5]*it[7];
-       itB[3] = it[5]*it[6] - it[3]*it[8];
-       itB[6] = it[3]*it[7] - it[4]*it[6];
-       J__ = it[0] * itB[0] + it[1] * itB[3] + it[2] * itB[6];
+       scalar_type a0 = itB[0] = it[4]*it[8] - it[5]*it[7];
+       scalar_type a1 = itB[1] = it[5]*it[6] - it[3]*it[8];
+       scalar_type a2 = itB[2] = it[3]*it[7] - it[4]*it[6];
+       J__ = it[0] * a0 + it[1] * a1 + it[2] * a2;
        J_ = gmm::abs(J__);
       } else {
        B_factors.base_resize(P, P); // store factorization for B computation
@@ -284,16 +284,16 @@
         B_(0, 0) = scalar_type(1) / J__;
       } else if (P == 2) {
        auto it = KK.begin(); auto itB = B_.begin();
-       *itB++ = it[3] / J__; *itB++ = it[2] / J__; 
-       *itB++ = it[1] / J__; *itB = (*it) / J__;
+       *itB++ = it[3] / J__; *itB++ = -it[2] / J__; 
+       *itB++ = -it[1] / J__; *itB = (*it) / J__;
       } else if (P == 3) {
        auto it = KK.begin(); auto itB = B_.begin();
-       itB[0] /= J__; itB[3] /= J__; itB[6] /= J__; 
-       itB[1] = (it[2]*it[7] - it[1]*it[8]) / J__;
-       itB[2] = (it[1]*it[5] - it[2]*it[4]) / J__;
+       itB[0] /= J__; itB[1] /= J__; itB[2] /= J__; 
+       itB[3] = (it[2]*it[7] - it[1]*it[8]) / J__;
+       itB[6] = (it[1]*it[5] - it[2]*it[4]) / J__;
        itB[4] = (it[0]*it[8] - it[2]*it[6]) / J__;
-       itB[5] = (it[2]*it[3] - it[0]*it[5]) / J__;
-       itB[7] = (it[1]*it[6] - it[0]*it[7]) / J__;
+       itB[7] = (it[2]*it[3] - it[0]*it[5]) / J__;
+       itB[5] = (it[1]*it[6] - it[0]*it[7]) / J__;
        itB[8] = (it[0]*it[4] - it[1]*it[3]) / J__;
       } else {
        bgeot::lu_inverse(&(*(B_factors.begin())), ipvt, &(*(B_.begin())), P);