getfem-commits
[Top][All Lists]
Advanced

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

[Getfem-commits] r5410 - in /trunk/getfem: contrib/opt_assembly/ src/ sr


From: Yves . Renard
Subject: [Getfem-commits] r5410 - in /trunk/getfem: contrib/opt_assembly/ src/ src/gmm/
Date: Wed, 12 Oct 2016 12:21:03 -0000

Author: renard
Date: Wed Oct 12 14:21:01 2016
New Revision: 5410

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

Modified:
    trunk/getfem/contrib/opt_assembly/opt_assembly.cc
    trunk/getfem/src/bgeot_geometric_trans.cc
    trunk/getfem/src/getfem_generic_assembly.cc
    trunk/getfem/src/gmm/gmm_vector.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=5410&r1=5409&r2=5410&view=diff
==============================================================================
--- trunk/getfem/contrib/opt_assembly/opt_assembly.cc   (original)
+++ trunk/getfem/contrib/opt_assembly/opt_assembly.cc   Wed Oct 12 14:21:01 2016
@@ -433,29 +433,29 @@
   // - Instructions execution except for assembly ones
   //                        new  | old  | sto  | asse | exec |  J   | Ins  |
   test_new_assembly(2, 400, 1); // ndofu = 321602 ndofp = 160801 ndofchi = 1201
-  // Mass (scalar)        : 0.37 | 0.61 | 0.07 | 0.10 | 0.21 | 0.06 | 0.06 |
+  // Mass (scalar)        : 0.36 | 0.61 | 0.07 | 0.10 | 0.20 | 0.06 | 0.06 |
   // Mass (vector)        : 0.46 | 0.82 | 0.12 | 0.22 | 0.15 | 0.05 | 0.09 |
   // Laplacian            : 0.26 | 0.83 | 0.05 | 0.07 | 0.14 | 0.06 | 0.05 |
-  // Homogeneous elas     : 0.47 | 1.88 | 0.14 | 0.22 | 0.14 | 0.05 | 0.11 |
-  // Non-homogeneous elast: 0.58 | 2.26 | 0.14 | 0.22 | 0.14 | 0.05 | 0.21 |
+  // Homogeneous elas     : 0.46 | 1.88 | 0.14 | 0.22 | 0.13 | 0.05 | 0.11 |
+  // Non-homogeneous elast: 0.57 | 2.26 | 0.14 | 0.22 | 0.13 | 0.05 | 0.21 |
   test_new_assembly(3, 36, 1);  // ndofu = 151959 ndofp =  50653 ndofchi = 6553
   // Mass (scalar)        : 0.44 | 0.77 | 0.08 | 0.12 | 0.24 | 0.11 | 0.08 |
   // Mass (vector)        : 0.94 | 1.57 | 0.18 | 0.39 | 0.20 | 0.11 | 0.35 |
-  // Laplacian            : 0.42 | 1.38 | 0.06 | 0.09 | 0.19 | 0.11 | 0.14 |
+  // Laplacian            : 0.41 | 1.38 | 0.06 | 0.09 | 0.19 | 0.11 | 0.14 |
   // Homogeneous elas     : 1.28 | 4.58 | 0.48 | 0.59 | 0.19 | 0.11 | 0.51 |
-  // Non-homogeneous elast: 1.44 | 6.81 | 0.49 | 0.63 | 0.19 | 0.11 | 0.62 |
+  // Non-homogeneous elast: 1.42 | 6.81 | 0.49 | 0.63 | 0.18 | 0.11 | 0.61 |
   test_new_assembly(2, 200, 2); // ndofu = 321602 ndofp = 160801 ndofchi = 1201
   // Mass (scalar)        : 0.15 | 0.25 | 0.04 | 0.07 | 0.05 | 0.02 | 0.03 |
   // Mass (vector)        : 0.34 | 0.45 | 0.07 | 0.15 | 0.05 | 0.02 | 0.14 |
   // Laplacian            : 0.13 | 0.37 | 0.03 | 0.05 | 0.04 | 0.02 | 0.04 |
   // Homogeneous elas     : 0.37 | 1.28 | 0.13 | 0.19 | 0.04 | 0.01 | 0.14 |
-  // Non-homogeneous elast: 0.46 | 2.40 | 0.13 | 0.18 | 0.04 | 0.01 | 0.24 |
+  // Non-homogeneous elast: 0.45 | 2.40 | 0.13 | 0.18 | 0.04 | 0.01 | 0.23 |
   test_new_assembly(3, 18, 2);  // ndofu = 151959 ndofp =  50653 ndofchi = 6553
   // Mass (scalar)        : 0.23 | 0.30 | 0.10 | 0.15 | 0.04 | 0.02 | 0.04 |
   // Mass (vector)        : 1.46 | 0.90 | 0.33 | 0.64 | 0.04 | 0.02 | 0.78 |
   // Laplacian            : 0.18 | 0.55 | 0.08 | 0.10 | 0.03 | 0.02 | 0.05 |
   // Homogeneous elas     : 2.37 | 3.47 | 1.20 | 1.37 | 0.03 | 0.02 | 0.97 |
-  // Non-homogeneous elast: 2.46 | 9.25 | 1.20 | 1.37 | 0.03 | 0.02 | 1.06 |
+  // Non-homogeneous elast: 2.45 | 9.25 | 1.20 | 1.37 | 0.03 | 0.02 | 1.05 |
   test_new_assembly(3, 9, 4);   // ndofu = 151959 ndofp =  50653 ndofchi = 6553
   // Mass (scalar)        : 0.75 | 0.38 | 0.25 | 0.39 | 0.01 | .005 | 0.35 |
   // Mass (vector)        : 5.10 | 1.33 | 0.42 | 1.86 | 0.01 | .005 | 3.23 |

Modified: trunk/getfem/src/bgeot_geometric_trans.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/bgeot_geometric_trans.cc?rev=5410&r1=5409&r2=5410&view=diff
==============================================================================
--- trunk/getfem/src/bgeot_geometric_trans.cc   (original)
+++ trunk/getfem/src/bgeot_geometric_trans.cc   Wed Oct 12 14:21:01 2016
@@ -36,7 +36,7 @@
   const base_node& geotrans_interpolation_context::xref() const {
     if (!have_xref()) {
       if (pspt_) xref_ = (*pspt_)[ii_];
-      else GMM_ASSERT1(false, "missing xref");
+      else GMM_ASSERT1(false, "Missing xref");
     }
     return xref_;
   }
@@ -51,7 +51,7 @@
   }
 
   void geotrans_interpolation_context::compute_J(void) const {
-    GMM_ASSERT1(have_G() && have_pgt(), "unable to compute J\n");
+    GMM_ASSERT1(have_G() && have_pgt(), "Unable to compute J\n");
     size_type P = pgt_->structure()->dim();
     if (P != N()) {
       B_factors.base_resize(P, P);
@@ -80,15 +80,15 @@
 
   const base_matrix& geotrans_interpolation_context::K() const {
     if (!have_K()) {
-      GMM_ASSERT1(have_G() && have_pgt(), "unable to compute K\n");
+      GMM_ASSERT1(have_G() && have_pgt(), "Unable to compute K\n");
       size_type P = pgt_->structure()->dim();
       K_.base_resize(N(), P);
       if (have_pgp()) {
-             pgt_->compute_K_matrix(G(), pgp_->grad(ii_), K_);
+       pgt_->compute_K_matrix(G(), pgp_->grad(ii_), K_);
       } else {
        PC.base_resize(pgt()->nb_points(), P);
         pgt()->poly_vector_grad(xref(), PC);
-        pgt_->compute_K_matrix(G(), PC, K_);
+       pgt_->compute_K_matrix(G(), PC, K_);
       }
       have_K_ = true;
     }
@@ -97,7 +97,7 @@
 
   const base_matrix& geotrans_interpolation_context::B() const {
     if (!have_B()) {
-      GMM_ASSERT1(have_G() && have_pgt(), "unable to compute B\n");
+      GMM_ASSERT1(have_G() && have_pgt(), "Unable to compute B\n");
       size_type P = pgt_->structure()->dim();
       const base_matrix &KK = K();
       B_.base_resize(N(), P);
@@ -108,20 +108,20 @@
         gmm::mult(KK, PC, B_);
       } else if (P == 1) {
         scalar_type det = KK(0, 0);
-        GMM_ASSERT1(det != scalar_type(0), "non invertible matrix");
+        GMM_ASSERT1(det != scalar_type(0), "Non invertible matrix");
         B_(0, 0) = scalar_type(1)/det;
         J_ = gmm::abs(det);
       } 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");
+        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);
       } else {
         scalar_type det = J();
-        GMM_ASSERT1(det != scalar_type(0), "non invertible matrix");
+        GMM_ASSERT1(det != scalar_type(0), "Non invertible matrix");
         gmm::lu_inverse(B_factors, ipvt, B_);
       }
       have_B_ = true;

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=5410&r1=5409&r2=5410&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Wed Oct 12 14:21:01 2016
@@ -2749,7 +2749,7 @@
       size_type cv_1 = ctx.is_convex_num_valid()
                      ? ctx.convex_num() : mf.convex_index().first_true();
       pfem pf = mf.fem_of_element(cv_1);
-      GMM_ASSERT1(pf, "An element without finite element methode defined");
+      GMM_ASSERT1(pf, "An element without finite element method defined");
       size_type Qmult = qdim / pf->target_dim();
       size_type s = pf->nb_dof(cv_1) * Qmult;
       if (t.sizes()[0] != s)
@@ -6481,10 +6481,16 @@
 
 
   void ga_workspace::assembly(size_type order) {
+    size_type ndof;
+    const ga_workspace *w = this;
+    while (w->parent_workspace) w = w->parent_workspace;
+    if (w->md) ndof = w->md->nb_dof(); // To eventually call actualize_sizes()
+
     GA_TIC;
     ga_instruction_set gis;
     ga_compile(*this, gis, order);
-    size_type ndof = gis.nb_dof, max_dof =  gis.max_dof;
+    ndof = gis.nb_dof;
+    size_type max_dof =  gis.max_dof;
     GA_TOCTIC("Compile time");
 
     if (order == 2) {
@@ -11721,7 +11727,7 @@
 
   static void ga_exec(ga_instruction_set &gis, ga_workspace &workspace) {
     base_matrix G;
-    base_small_vector un, up;
+    base_small_vector un;
     scalar_type J(0);
 
     for (const std::string &t : gis.transformations)
@@ -11739,10 +11745,11 @@
       mesh_region rg(region);
       m.intersect_with_mpi_region(rg);
       size_type old_cv = size_type(-1);
-      bgeot::pgeometric_trans pgt = 0;
+      bgeot::pgeometric_trans pgt = 0, pgt_old = 0;
       pintegration_method pim = 0;
       papprox_integration pai = 0;
       bgeot::pstored_point_tab pspt = 0, old_pspt = 0;
+      bgeot::pgeotrans_precomp pgp = 0;
       bool first_gp = true;
       for (getfem::mr_visitor v(rg, m); !v.finished(); ++v) {
         if (mim.convex_index().is_in(v.cv())) {
@@ -11750,27 +11757,27 @@
           if (v.cv() != old_cv) {
             m.points_of_convex(v.cv(), G);
             pgt = m.trans_of_convex(v.cv());
-            up.resize(G.nrows());
-            un.resize(pgt->dim());
-            pim = mim.int_method_of_element(v.cv());
+           pim = mim.int_method_of_element(v.cv());
             if (pim->type() == IM_NONE) continue;
             GMM_ASSERT1(pim->type() == IM_APPROX, "Sorry, exact methods cannot 
"
                         "be used in high level generic assembly");
             pai = pim->approx_method();
             pspt = pai->pintegration_points();
             if (pspt->size()) {
-              if (gis.ctx.have_pgp() && gis.pai == pai && gis.ctx.pgt()==pgt) {
-                gis.ctx.change(gis.ctx.pgp(), 0, 0, G, v.cv(), v.f());
+              if (pgp && gis.pai == pai && pgt_old == pgt) {
+                gis.ctx.change(pgp, 0, 0, G, v.cv(), v.f());
               } else {
                 if (pai->is_built_on_the_fly()) {
                   gis.ctx.change(pgt, 0, (*pspt)[0], G, v.cv(), v.f());
+                 pgp = 0;
                 } else {
-                  gis.ctx.change(gis.gp_pool(pgt, pspt), 0,0, G, v.cv(), 
v.f());
+                 pgp = gis.gp_pool(pgt, pspt);
+                  gis.ctx.change(pgp, 0, 0, G, v.cv(), v.f());
                 }
+               pgt_old = pgt; gis.pai = pai;
               }
-              gis.pai = pai;
-              if (gis.need_elt_size)
-                gis.elt_size = m.convex_radius_estimate(v.cv())*scalar_type(2);
+              if (gis.need_elt_size) 
+                gis.elt_size = convex_radius_estimate(pgt, G)*scalar_type(2);
             }
             old_cv = v.cv();
           } else {
@@ -11787,19 +11794,20 @@
               first_ind = pai->ind_first_point_on_face(v.f());
             }
             for (gis.ipt = 0; gis.ipt < gis.nbpt; ++(gis.ipt)) {
-              if (gis.ctx.have_pgp()) gis.ctx.set_ii(first_ind+gis.ipt);
+              if (pgp) gis.ctx.set_ii(first_ind+gis.ipt);
               else gis.ctx.set_xref((*pspt)[first_ind+gis.ipt]);
               if (gis.ipt == 0 || !(pgt->is_linear())) {
                 J = gis.ctx.J();
                 // Computation of unit normal vector in case of a boundary
                 if (v.f() != short_type(-1)) {
-                  gmm::copy(pgt->normals()[v.f()], un);
-                  gmm::mult(gis.ctx.B(), un, up);
-                  scalar_type nup = gmm::vect_norm2(up);
+                 gis.Normal.resize(G.nrows());
+                 un.resize(pgt->dim());
+                 gmm::copy(pgt->normals()[v.f()], un);
+                  gmm::mult(gis.ctx.B(), un, gis.Normal);
+                  scalar_type nup = gmm::vect_norm2(gis.Normal);
                   J *= nup;
-                  gmm::scale(up,1.0/nup);
-                  gmm::clean(up, 1e-13);
-                  gis.Normal = up;
+                  gmm::scale(gis.Normal,1.0/nup);
+                  gmm::clean(gis.Normal, 1e-13);
                 } else gis.Normal.resize(0);
               }
               gis.coeff = J * pai->coeff(first_ind+gis.ipt);

Modified: trunk/getfem/src/gmm/gmm_vector.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/gmm/gmm_vector.h?rev=5410&r1=5409&r2=5410&view=diff
==============================================================================
--- trunk/getfem/src/gmm/gmm_vector.h   (original)
+++ trunk/getfem/src/gmm/gmm_vector.h   Wed Oct 12 14:21:01 2016
@@ -1065,8 +1065,8 @@
          base_type_::resize(this->nb_stored()+1, ev);
          if (ind != this->nb_stored() - 1) {
            it = this->begin() + ind;
-           iterator ite = this->end(); --ite; iterator itee = ite; --itee; 
-           for (; ite != it; --ite, --itee) *ite = *itee;
+           iterator ite = this->end(); --ite; iterator itee = ite; 
+           for (; ite != it; --ite) { --itee; *ite = *itee; }
            *it = ev;
          }
        }
@@ -1092,8 +1092,8 @@
          base_type_::resize(this->nb_stored()+1, ev);
          if (ind != this->nb_stored() - 1) {
            it = this->begin() + ind;
-           iterator ite = this->end(); --ite; iterator itee = ite; --itee; 
-           for (; itee != it; --ite, --itee) *ite = *itee;
+           iterator ite = this->end(); --ite; iterator itee = ite; 
+           for (; ite != it; --ite) { --itee; *ite = *itee; }
            *it = ev;
          }
        }




reply via email to

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