[Getfem-commits] r4688 - in /trunk/getfem/src: ./ getfem/

[logari81]

Author: logari81
Date: Thu Jun 19 09:40:14 2014
New Revision: 4688

URL: http://svn.gna.org/viewcvs/getfem?rev=4688&view=rev
Log:
add optional support for frame indifferent relative velocity in large sliding 
contact

Modified:
    trunk/getfem/src/getfem/getfem_contact_and_friction_common.h
    trunk/getfem/src/getfem_contact_and_friction_common.cc
    trunk/getfem/src/getfem_contact_and_friction_integral.cc
    trunk/getfem/src/getfem_contact_and_friction_large_sliding.cc
    trunk/getfem/src/getfem_generic_assembly.cc

Modified: trunk/getfem/src/getfem/getfem_contact_and_friction_common.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_contact_and_friction_common.h?rev=4688&r1=4687&r2=4688&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_contact_and_friction_common.h        
(original)
+++ trunk/getfem/src/getfem/getfem_contact_and_friction_common.h        Thu Jun 
19 09:40:14 2014
@@ -401,6 +401,10 @@
                                   size_type ndof, size_type qdim, size_type N);
 
 
+  void compute_normal(const fem_interpolation_context &ctx, size_type face,
+                      bool in_reference_conf, const model_real_plain_vector 
&coeff,
+                      base_node &n0, base_node &n, base_matrix &grad);
+
   //=========================================================================
   //
   //  Structure which stores the contact boundaries, rigid obstacles and

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=4688&r1=4687&r2=4688&view=diff
==============================================================================
--- trunk/getfem/src/getfem_contact_and_friction_common.cc      (original)
+++ trunk/getfem/src/getfem_contact_and_friction_common.cc      Thu Jun 19 
09:40:14 2014
@@ -34,8 +34,8 @@
 
   void compute_normal(const fem_interpolation_context &ctx,
                       size_type face, bool in_reference_conf,
+                      const model_real_plain_vector &coeff,
                       base_node &n0, base_node &n,
-                      model_real_plain_vector &coeff,
                       base_matrix &grad) {
       n0 = bgeot::compute_normal(ctx, face);
       if (in_reference_conf) {
@@ -457,8 +457,7 @@
             }
 
             // unit normal vector computation
-            compute_normal(ctx, v.f(), ref_conf,
-                           n0, n, coeff, grad);
+            compute_normal(ctx, v.f(), ref_conf, coeff, n0, n, grad);
             n /= gmm::vect_norm2(n);
 
             if (on_fem_nodes && dof_already_interpolated[ind]) {
@@ -653,8 +652,7 @@
 
             // computation of unit normal vector if the vertex is on the face
             if (points_on_face[ip]) {
-              compute_normal(ctx, v.f(), ref_conf,
-                             n0, n, coeff, grad);
+              compute_normal(ctx, v.f(), ref_conf, coeff, n0, n, grad);
               n /= gmm::vect_norm2(n);
               n_mean += n;
               ++nb_pt_on_face;
@@ -1230,7 +1228,7 @@
 
         // compute the unit normal vector at y and the signed distance.
         base_small_vector ny0(N);
-        compute_normal(ctx, iff, ref_conf, ny0, ny, coeff, grad);
+        compute_normal(ctx, iff, ref_conf, coeff, ny0, ny, grad);
         // ny /= gmm::vect_norm2(ny); // Useful only if the unit normal is kept
         signed_dist *= gmm::sgn(gmm::vect_sp(x - y, ny));
 

Modified: trunk/getfem/src/getfem_contact_and_friction_integral.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_contact_and_friction_integral.cc?rev=4688&r1=4687&r2=4688&view=diff
==============================================================================
--- trunk/getfem/src/getfem_contact_and_friction_integral.cc    (original)
+++ trunk/getfem/src/getfem_contact_and_friction_integral.cc    Thu Jun 19 
09:40:14 2014
@@ -490,7 +490,7 @@
             zt = (V - un1 * no) * (r * alpha) - zt;          // = zt1 - zt2 , 
with zt = r*alpha*u_T
           }
         }
-        un = un1 - un; // = un1 - unp1
+        un = un1 - un; // = un1 - un2
       }
       break;
 

Modified: trunk/getfem/src/getfem_contact_and_friction_large_sliding.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_contact_and_friction_large_sliding.cc?rev=4688&r1=4687&r2=4688&view=diff
==============================================================================
--- trunk/getfem/src/getfem_contact_and_friction_large_sliding.cc       
(original)
+++ trunk/getfem/src/getfem_contact_and_friction_large_sliding.cc       Thu Jun 
19 09:40:14 2014
@@ -246,6 +246,7 @@
   // For the moment, with raytrace detection and integral unsymmetric
   // Alart-Curnier augmented Lagrangian
 
+#undef CONSIDER_FRAME_INDIFFERENCE
 
   struct integral_large_sliding_contact_brick : public virtual_brick {
 
@@ -397,7 +398,7 @@
       GMM_ASSERT1(!isrigid(), "Rigid obstacle master node: no fem defined");
       if (!ctx_uy_init) {
         bgeot::vectors_to_base_matrix(Gy, meshy().points_of_convex(cvy_));
-        ctx_uy_ = fem_interpolation_context(pgty, pf_uy, y_ref(), Gy, cvy_,fy);
+        ctx_uy_ = fem_interpolation_context(pgty, pf_uy, y_ref(), Gy, cvy_, 
fy);
         ctx_uy_init = true;
       }
       return ctx_uy_;
@@ -544,8 +545,8 @@
     }
 
     scalar_type alpha;
-    base_small_vector x0_, y0_, Vs_;
-    bool x0_init, y0_init, Vs_init;
+    base_small_vector x0_, y0_, nx0_, Vs_;
+    bool x0_init, y0_init, nx0_init, Vs_init;
     base_matrix grad_phiy0_;
     bool grad_phiy0_init;
 
@@ -556,7 +557,7 @@
           pfem pf = ctx_ux().pf();
           slice_vector_on_basic_dof_of_element(*mf_ux_, all_x0, cvx_, coeff);
           pf->interpolation(ctx_ux(), coeff, x0_, dim_type(N));
-        }  else gmm::clear(x0_);
+        } else gmm::clear(x0_);
         gmm::add(ctx_ux().xreal(), x0_);
         x0_init = true;
       }
@@ -571,7 +572,7 @@
             pfem pf = ctx_uy().pf();
             slice_vector_on_basic_dof_of_element(*mf_uy_, all_y0, cvy_, coeff);
             pf->interpolation(ctx_uy(), coeff, y0_, dim_type(N));
-          }  else gmm::clear(y0_);
+          } else gmm::clear(y0_);
           gmm::add(ctx_uy().xreal(), y0_);
         } else gmm::copy(y(), y0_);
         y0_init = true;
@@ -579,12 +580,33 @@
       return y0_;
     }
 
+    const base_small_vector &nx0(void) {
+      if (!nx0_init) {
+        const model_real_plain_vector &all_x0 = mcf.disp0_of_boundary(ibx);
+        if (all_x0.size()) {
+          pfem pf = ctx_ux().pf();
+          slice_vector_on_basic_dof_of_element(*mf_ux_, all_x0, cvx_, coeff);
+          base_small_vector nx00_(N);
+          base_matrix grad_phix0_(N,N);
+          compute_normal(ctx_ux(), fx, false, coeff, nx00_, nx0_, grad_phix0_);
+          nx0_ /= gmm::vect_norm2(nx0_);
+        } else gmm::clear(nx0_);
+        nx0_init = true;
+      }
+      return nx0_;
+    }
+
     const base_small_vector &Vs(void) { // relative velocity
       if (!Vs_init) {
         if (alpha != scalar_type(0)) {
+#ifdef CONSIDER_FRAME_INDIFFERENCE
+          gmm::add(y0(), gmm::scaled(x0(), scalar_type(-1)), Vs_);
+          gmm::add(gmm::scaled(nx0(), -g()), Vs_);
+#else
           gmm::add(x(), gmm::scaled(y(), scalar_type(-1)), Vs_);
           gmm::add(gmm::scaled(x0(), scalar_type(-1)), Vs_);
           gmm::add(y0(), Vs_);
+#endif
           gmm::scale(Vs_, alpha);
         } else gmm::clear(Vs_);
         Vs_init = true;
@@ -680,7 +702,7 @@
       lambda_x_(N), lambda_y_(N),
       grad_phix_(N, N), grad_phix_inv_(N, N),
       grad_phiy_(N, N), grad_phiy_inv_(N, N), alpha(alpha_),
-      x0_(N), y0_(N), Vs_(N), grad_phiy0_(N, N), un(N) {}
+      x0_(N), y0_(N), nx0_(N), Vs_(N), grad_phiy0_(N, N), un(N) {}
 
   };
 
@@ -978,7 +1000,7 @@
         // second sub term
         gmm::mult(gpp.I_nxnx(), gmm::scaled(ny, -g), auxN1);
         gmm::mult(gpp.grad_phix_inv(), auxN1, auxN2);
-        gmm::mult_add(graddeltaunx, auxN2, auxUX);
+        gmm::mult_add(graddeltaunx, auxN2, auxUX);    // auxUX -> \delta u(X) 
- g Dn_x
         gmm::rank_one_update(Melem, deltamudgF,  auxUX);
         gmm::scale(Melem, weight*FMULT);
         mat_elem_assembly(M, I_lx, I_ux, Melem, *mf_lx, cvx, *mf_ux, cvx);
@@ -987,7 +1009,7 @@
           // LXUY term
           // third sub term
           gmm::resize(Melem, ndof_lx, ndof_uy); gmm::clear(Melem);
-          gmm::mult(gpp.vbase_uy(), ny, auxUY);
+          gmm::mult(gpp.vbase_uy(), ny, auxUY);       // auxUY -> \delta u(Y)
           gmm::rank_one_update(Melem, deltamudgF, auxUY);
           gmm::scale(Melem, -weight*FMULT);
           mat_elem_assembly(M, I_lx, I_uy, Melem, *mf_lx, cvx, *mf_uy, cvy);
@@ -997,6 +1019,40 @@
           // Term -(1/r) d_Vs F \delta Vs\delta\mu
 
           if (!isrigid) {
+#ifdef CONSIDER_FRAME_INDIFFERENCE
+            base_matrix gphiy0gphiyinv(N, N);
+            gmm::mult(gpp.grad_phiy0(), gpp.grad_phiy_inv(), gphiy0gphiyinv);
+            gmm::mult(gphiy0gphiyinv, gpp.I_nxny(), auxNN1);
+            gmm::rank_one_update(auxNN1, gpp.nx0(),
+                                 
gmm::scaled(gpp.ny(),scalar_type(1)/gpp.nxdotny()));
+            gmm::mult(dVsF, auxNN1, auxNN2);
+            // Caution: auxNN2 re-used in the second sub term
+
+            // LXUX term
+            // first sub term
+            gmm::resize(Melem, ndof_lx, ndof_ux); gmm::clear(Melem);
+            gmm::mult(gpp.vbase_lx(), gmm::transposed(auxNN2), auxLXN1);
+            // Caution: auxLXN1 re-used in the third sub term
+            gmm::mult(auxLXN1, gmm::transposed(gpp.vbase_ux()), Melem);
+
+            // second sub term
+            base_matrix auxLXUX(ndof_lx, ndof_ux);
+            gmm::mult(auxNN2, gpp.I_nxnx(), auxNN1);
+            gmm::mult(auxNN1, gmm::transposed(gpp.grad_phix_inv()), auxNN2);
+            gmm::mult(gpp.vbase_lx(), gmm::transposed(auxNN2), auxLXN2);
+            gmm::mult(auxLXN2, gmm::transposed(graddeltaunx), auxLXUX);
+            gmm::scale(auxLXUX, -g);
+            gmm::add(auxLXUX, Melem);
+            gmm::scale(Melem, -weight*alpha*FMULT/r);
+            mat_elem_assembly(M, I_lx, I_ux, Melem, *mf_lx, cvx, *mf_ux, cvx);
+
+            // LXUY term
+            // third sub term
+            gmm::resize(Melem, ndof_lx, ndof_uy); gmm::clear(Melem);
+            // Caution: auxLXN1 re-used
+            gmm::mult(auxLXN1, gmm::transposed(gpp.vbase_uy()), Melem);
+            gmm::scale(Melem, weight*alpha*FMULT/r);
+#else
             base_matrix I_gphiy0gphiyinv(N, N);
             gmm::mult(gmm::scaled(gpp.grad_phiy0(), scalar_type(-1)),
                       gpp.grad_phiy_inv(), I_gphiy0gphiyinv);
@@ -1027,14 +1083,11 @@
             // LXUY term
             // third sub term
             gmm::resize(Melem, ndof_lx, ndof_uy); gmm::clear(Melem);
-//             gmm::mult(I_gphiy0gphiyinv, gpp.I_nxny(), auxNN1);
-//             for (size_type j = 0; j < N; ++j) auxNN1(j,j) -= scalar_type(1);
-//             gmm::mult(dVsF, auxNN1, auxNN2);
-//             gmm::mult(gpp.vbase_lx(), gmm::transposed(auxNN2), auxLXN2);
             // Caution: auxLXN2 re-used
             gmm::mult(auxLXN2, gmm::transposed(gpp.vbase_uy()), Melem);
             gmm::scale(Melem, -weight*alpha*FMULT/r);
             mat_elem_assembly(M, I_lx, I_uy, Melem, *mf_lx, cvx, *mf_uy, cvy);
+#endif
           } else {
             // LXUX term
             gmm::mult(gpp.vbase_lx(), gmm::transposed(dVsF), auxLXN1);

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=4688&r1=4687&r2=4688&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Thu Jun 19 09:40:14 2014
@@ -335,7 +335,7 @@
       if (test0 && test1 && (test0 == test1 ||
                              test0 >= 3 || test1 >= 3))
         ga_throw_error(expr, pos, "Incompatibility of test functions "
-                        " in product.");
+                        "in product.");
       GMM_ASSERT1(test0 != size_type(-1) && test1 != size_type(-1),
                   "internal error");