[Getfem-commits] r4935 - /trunk/getfem/src/getfem_contact_and_friction_common.cc

[logari81]

Author: logari81
Date: Sun Apr  5 23:35:05 2015
New Revision: 4935

URL: http://svn.gna.org/viewcvs/getfem?rev=4935&view=rev
Log:
fix bug in computation of influence boxes in contact and add caching of normals

Modified:
    trunk/getfem/src/getfem_contact_and_friction_common.cc

Modified: trunk/getfem/src/getfem_contact_and_friction_common.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_contact_and_friction_common.cc?rev=4935&r1=4934&r2=4935&view=diff
==============================================================================
--- trunk/getfem/src/getfem_contact_and_friction_common.cc      (original)
+++ trunk/getfem/src/getfem_contact_and_friction_common.cc      Sun Apr  5 
23:35:05 2015
@@ -605,6 +605,7 @@
 
         dal::bit_vector points_already_interpolated;
         std::vector<base_node> transformed_points(m.nb_max_points());
+        std::vector<base_node> normals(m.nb_max_points());
         for (getfem::mr_visitor v(region,m); !v.finished(); ++v) {
           size_type cv = v.cv();
           bgeot::pgeometric_trans pgt = m.trans_of_convex(cv);
@@ -638,8 +639,20 @@
                 transformed_points[ind] = ctx.xreal();
               }
               points_already_interpolated.add(ind);
+              // computation of unit normal vector if the vertex is on the face
+              if (points_on_face[ip]) {
+                compute_normal(ctx, v.f(), ref_conf, coeff, n0, n, grad);
+                n /= gmm::vect_norm2(n);
+                normals[ind] = n;
+                n_mean += n;
+                ++nb_pt_on_face;
+              }
             } else {
               val = transformed_points[ind];
+              if (points_on_face[ip]) {
+                n_mean += normals[ind];
+                ++nb_pt_on_face;
+              }
             }
 
             if (ip == 0) // computation of bounding box
@@ -649,14 +662,6 @@
                 bmin[k] = std::min(bmin[k], val[k]);
                 bmax[k] = std::max(bmax[k], val[k]);
               }
-            }
-
-            // computation of unit normal vector if the vertex is on the face
-            if (points_on_face[ip]) {
-              compute_normal(ctx, v.f(), ref_conf, coeff, n0, n, grad);
-              n /= gmm::vect_norm2(n);
-              n_mean += n;
-              ++nb_pt_on_face;
             }
 
           }
@@ -1389,6 +1394,7 @@
           
           dal::bit_vector points_already_interpolated;
           std::vector<base_node> transformed_points(m.nb_max_points());
+          std::vector<base_node> normals(m.nb_max_points());
           for (getfem::mr_visitor v(region,m); !v.finished(); ++v) {
             size_type cv = v.cv();
             bgeot::pgeometric_trans pgt = m.trans_of_convex(cv);
@@ -1397,7 +1403,7 @@
             slice_vector_on_basic_dof_of_element(mfu, U, cv, coeff);
             bgeot::vectors_to_base_matrix
               (G, mfu.linked_mesh().points_of_convex(cv));
-            fem_interpolation_context ctx(pgt,pfp,size_type(-1), G, cv);
+            fem_interpolation_context ctx(pgt, pfp, size_type(-1), G, cv);
             
             bgeot::pconvex_structure cvs = pgt->structure();
             size_type nb_pt_on_face = cvs->nb_points_of_face(v.f());
@@ -1416,8 +1422,14 @@
                 val += ctx.xreal();
                 transformed_points[ind] = val;
                 points_already_interpolated.add(ind);
+                // computation of unit normal vector if the vertex is on the 
face
+                compute_normal(ctx, v.f(), false, coeff, n0_x, n_x, grad);
+                n_x /= gmm::vect_norm2(n_x);
+                normals[ind] = n_x;
+                n_mean += n_x;
               } else {
                 val = transformed_points[ind];
+                n_mean += normals[ind];
               }
               
               if (k == 0) // computation of bounding box
@@ -1428,11 +1440,6 @@
                   bmax[l] = std::max(bmax[l], val[l]);
                 }
               }
-              
-              // computation of unit normal vector if the vertex is on the face
-              compute_normal(ctx, v.f(), false, coeff, n0_x, n_x, grad);
-              n_x /= gmm::vect_norm2(n_x);
-              n_mean += n_x;
             }
             
             // Security coefficient of 1.3 (for nonlinear transformations)