getfem-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Getfem-commits] r5442 - in /trunk/getfem: contrib/opt_assembly/opt_asse


From: Yves . Renard
Subject: [Getfem-commits] r5442 - in /trunk/getfem: contrib/opt_assembly/opt_assembly.cc src/bgeot_geometric_trans.cc src/getfem/bgeot_geometric_trans.h
Date: Wed, 26 Oct 2016 08:06:53 -0000

Author: renard
Date: Wed Oct 26 10:06:52 2016
New Revision: 5442

URL: http://svn.gna.org/viewcvs/getfem?rev=5442&view=rev
Log:
direct formula for 3x3 matrix invert

Modified:
    trunk/getfem/contrib/opt_assembly/opt_assembly.cc
    trunk/getfem/src/bgeot_geometric_trans.cc
    trunk/getfem/src/getfem/bgeot_geometric_trans.h

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=5442&r1=5441&r2=5442&view=diff
==============================================================================
--- trunk/getfem/contrib/opt_assembly/opt_assembly.cc   (original)
+++ trunk/getfem/contrib/opt_assembly/opt_assembly.cc   Wed Oct 26 10:06:52 2016
@@ -256,7 +256,7 @@
   
   bool all = false;
   bool select = true;
-  int only_one = 11;
+  int only_one = 2;
 
   if (all || select || only_one == 1) {
     VEC_TEST_1("Test for source term", ndofu, "u.Test_u", mim, size_type(-1),
@@ -462,13 +462,13 @@
   // Non-homogeneous elast: 0.36 | 2.26 | 0.09 | 0.16 | 0.06 | 0.14 |
   if (all || only_one == 2) // ndofu = 151959 ndofp =  50653 ndofchi = 6553
     test_new_assembly(3, 36, 1);
-  // Vector source term   : 0.33 | 1.00 |
-  // Nonlinear residual   : 0.84 |      |
-  // Mass (scalar)        : 0.23 | 0.75 | 0.05 | 0.08 | 0.09 | 0.06 |
-  // Mass (vector)        : 0.70 | 1.54 | 0.17 | 0.27 | 0.09 | 0.34 |
-  // Laplacian            : 0.23 | 1.37 | 0.03 | 0.07 | 0.09 | 0.07 |
-  // Homogeneous elas     : 0.85 | 4.58 | 0.26 | 0.33 | 0.09 | 0.43 |
-  // Non-homogeneous elast: 0.91 | 6.55 | 0.26 | 0.33 | 0.09 | 0.49 |
+  // Vector source term   : 0.31 | 0.82 |
+  // 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 |
   if (all || only_one == 3) // ndofu = 321602 ndofp = 160801 ndofchi = 1201
     test_new_assembly(2, 200, 2);
   // Vector source term   : 0.11 | 0.24 |
@@ -480,9 +480,9 @@
   // Non-homogeneous elast: 0.30 | 2.38 | 0.07 | 0.10 | 0.03 | 0.17 |
   if (all || only_one == 4) // ndofu = 151959 ndofp =  50653 ndofchi = 6553
     test_new_assembly(3, 18, 2);
-  // Vector source term   : 0.17 | 0.26 |
-  // Nonlinear residual   : 0.50 |      |
-  // Mass (scalar)        : 0.11 | 0.28 | 0.05 | 0.05 | 0.03 | 0.03 |
+  // Vector source term   : 0.16 | 0.24 |
+  // Nonlinear residual   : 0.40 |      |
+  // 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 |
@@ -490,7 +490,7 @@
   if (all || only_one == 5) // ndofu = 151959 ndofp =  50653 ndofchi = 6553
     test_new_assembly(3, 9, 4);
   // Vector source term   : 0.13 | 0.20 |
-  // Nonlinear residual   : 0.43 |      |
+  // Nonlinear residual   : 0.38 |      |
   // 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 |

Modified: trunk/getfem/src/bgeot_geometric_trans.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/bgeot_geometric_trans.cc?rev=5442&r1=5441&r2=5442&view=diff
==============================================================================
--- trunk/getfem/src/bgeot_geometric_trans.cc   (original)
+++ trunk/getfem/src/bgeot_geometric_trans.cc   Wed Oct 26 10:06:52 2016
@@ -115,6 +115,10 @@
       return *A;
     } else if (N == 2) {
       return (*A) * (A[3]) - (A[1]) * (A[2]);
+    } 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];
+      return A[0] * a0 + A[1] * a3 + A[2] * a6;
     } else {
       size_type NN = N*N;
       if (__aux1.size() < NN) __aux1.resize(N*N);
@@ -147,6 +151,18 @@
       scalar_type det = a * d - b * c;
       GMM_ASSERT1(det != scalar_type(0), "Non invertible matrix");
       *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;
+      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;
       return det;
     } else {
       size_type NN = N*N;
@@ -199,8 +215,6 @@
     return xreal_;
   }
 
-
-
   void geotrans_interpolation_context::compute_J(void) const {
     GMM_ASSERT1(have_G() && have_pgt(), "Unable to compute J\n");
     size_type P = pgt_->structure()->dim();
@@ -211,35 +225,28 @@
       ipvt.resize(P);
       bgeot::lu_factor(&(*(B_factors.begin())), ipvt, P);
       // gmm::abs below because on flat convexes determinant could be -1e-27.
-      J_ = ::sqrt(gmm::abs(bgeot::lu_det(&(*(B_factors.begin())), ipvt, P)));
-    }
-    else {
+      J__=J_=::sqrt(gmm::abs(bgeot::lu_det(&(*(B_factors.begin())), ipvt, P)));
+    } else {
       // J_ = gmm::abs(gmm::lu_det(KK));
       if (P <= 2) {
        auto it = KK.begin();
-       if (P == 1) J_ = gmm::abs(*it);
-       else J_ = gmm::abs((*it) * (it[3]) - (it[1]) * (it[2]));
-      } else if (P == 3 && false) {
-       B_factors.base_resize(P, P); // co-factors
-       auto it = KK.begin();
-       auto itB = B_factors.begin();
-       *itB++ = it[4]*it[8]-it[5]*it[7];
-       *itB++ = it[2]*it[7]-it[1]*it[8];
-       *itB++ = it[1]*it[5]-it[2]*it[4];
-       *itB++ = it[5]*it[6]-it[3]*it[8];
-       *itB++ = it[0]*it[8]-it[2]*it[6];
-       *itB++ = it[2]*it[3]-it[0]*it[5];
-       *itB++ = it[3]*it[7]-it[4]*it[6];
-       *itB++ = it[1]*it[6]-it[0]*it[7];
-       *itB++ = it[0]*it[4]-it[1]*it[3];
-       itB = B_factors.begin();
-       J_ = it[0] * itB[0] + it[1] * itB[3] + it[2] * itB[6];
+       if (P == 1) { J__ = *it; J_ = gmm::abs(J__); } 
+       else { J__ = (*it) * (it[3]) - (it[1]) * (it[2]); J_ = gmm::abs(J__); } 
+      } 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];
+       J_ = gmm::abs(J__);
       } else {
        B_factors.base_resize(P, P); // store factorization for B computation
        gmm::copy(gmm::transposed(KK), B_factors);
        ipvt.resize(P);
        bgeot::lu_factor(&(*(B_factors.begin())), ipvt, P);
-       J_ = gmm::abs(bgeot::lu_det(&(*(B_factors.begin())), ipvt, P));
+       J__ = bgeot::lu_det(&(*(B_factors.begin())), ipvt, P);
+       J_ = gmm::abs(J__);
       }
     }
     have_J_ = true;
@@ -267,27 +274,28 @@
       const base_matrix &KK = K();
       size_type P = pgt_->structure()->dim(), N_ = gmm::mat_nrows(KK);
       B_.base_resize(N_, P);
+      if (!have_J_) compute_J();
+      GMM_ASSERT1(J__ != scalar_type(0), "Non invertible matrix");
       if (P != N_) {
-       if (!have_J_) compute_J();
        PC.base_resize(P, P);
         gmm::lu_inverse(B_factors, ipvt, PC);
         gmm::mult(KK, PC, B_);
       } else if (P == 1) {
-        scalar_type det = KK(0, 0);
-        GMM_ASSERT1(det != scalar_type(0), "Non invertible matrix");
-        B_(0, 0) = scalar_type(1)/det;
-        J_ = gmm::abs(det);
+        B_(0, 0) = scalar_type(1) / J__;
       } else if (P == 2) {
-        const scalar_type *p = &(KK(0,0));
-        scalar_type det = (*p) * (*(p+3)) - (*(p+1)) * (*(p+2));
-        GMM_ASSERT1(det != scalar_type(0), "Non invertible matrix");
-        scalar_type *q = &(B_(0,0));
-        *q++ =  (*(p+3)) / det;  *q++ = -(*(p+2)) / det;
-        *q++ = -(*(p+1)) / det;  *q++ =  (*p) / det;
-        J_ = gmm::abs(det);
+       auto it = KK.begin(); auto itB = B_.begin();
+       *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[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[8] = (it[0]*it[4] - it[1]*it[3]) / J__;
       } else {
-        scalar_type det = J();
-        GMM_ASSERT1(det != scalar_type(0), "Non invertible matrix");
        bgeot::lu_inverse(&(*(B_factors.begin())), ipvt, &(*(B_.begin())), P);
       }
       have_B_ = true;

Modified: trunk/getfem/src/getfem/bgeot_geometric_trans.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/bgeot_geometric_trans.h?rev=5442&r1=5441&r2=5442&view=diff
==============================================================================
--- trunk/getfem/src/getfem/bgeot_geometric_trans.h     (original)
+++ trunk/getfem/src/getfem/bgeot_geometric_trans.h     Wed Oct 26 10:06:52 2016
@@ -395,15 +395,15 @@
   */
   class APIDECL geotrans_interpolation_context {
   protected:
-    mutable base_node xref_; /** reference point */
+    mutable base_node xref_;  /** reference point */
     mutable base_node xreal_; /** transformed point */
-    const base_matrix *G_; /** pointer to the matrix of real nodes of the 
convex */
-    mutable base_matrix K_,B_, B3_, B32_; /** see documentation (getfem kernel 
doc) for more details */
+    const base_matrix *G_;    /** pointer to the matrix of real nodes of the 
convex */
+    mutable base_matrix K_, B_, B3_, B32_; /** see documentation (getfem 
kernel doc) for more details */
     pgeometric_trans pgt_;
     pgeotrans_precomp pgp_;
     pstored_point_tab pspt_; /** if pgp != 0, it is the same as pgp's one */
-    size_type ii_; /** index of current point in the pgp */
-    mutable scalar_type J_; /** Jacobian */
+    size_type ii_;           /** index of current point in the pgp        */
+    mutable scalar_type J_, J__; /** Jacobian */
     mutable base_matrix PC, B_factors;
     mutable base_vector aux1, aux2;
     mutable std::vector<int> ipvt;




reply via email to

[Prev in Thread] Current Thread [Next in Thread]