[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;
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] r5442 - in /trunk/getfem: contrib/opt_assembly/opt_assembly.cc src/bgeot_geometric_trans.cc src/getfem/bgeot_geometric_trans.h,
Yves . Renard <=