[Getfem-commits] r5291 - in /trunk/getfem: doc/sphinx/source/userdoc/gasm_high.rst src/getfem_plasticity.cc

[Yves . Renard]

Author: renard
Date: Fri Apr  8 21:19:00 2016
New Revision: 5291

URL: http://svn.gna.org/viewcvs/getfem?rev=5291&view=rev
Log:
Adding Normalized_reg(v, eps) operator

Modified:
    trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst
    trunk/getfem/src/getfem_plasticity.cc

Modified: trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst?rev=5291&r1=5290&r2=5291&view=diff
==============================================================================
--- trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst        (original)
+++ trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst        Fri Apr  8 
21:19:00 2016
@@ -567,7 +567,9 @@
 
   - ``Norm_sqr(v)`` for ``v`` a vector or a matrix gives the square of the 
euclidean norm of a vector or of the |Frobenius| norm of a matrix. For a vector 
this is equivalent to ``v.v`` and for a matrix to ``m:m``.
 
-  - ``Normalized(v)`` for ``v`` a vector or a matrix gives ``v`` divided by 
its euclidean (for vectors) or |Frobenius| (for matrices) norm.
+  - ``Normalized(v)`` for ``v`` a vector or a matrix gives ``v`` divided by 
its euclidean (for vectors) or |Frobenius| (for matrices) norm. In order to 
avoid problems when ``v`` is small, it is implemented as ``Normalized_reg(v, 
1E-25)``.
+
+  - ``Normalized_reg(v, eps)`` for ``v`` a vector or a matrix gives a 
regularized version of ``Normalized(v)`` : ``v/sqrt(|v|*|v|+eps*eps)``.
 
   - ``Det(m)`` gives the determinant of a square matrix ``m``.
 

Modified: trunk/getfem/src/getfem_plasticity.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_plasticity.cc?rev=5291&r1=5290&r2=5291&view=diff
==============================================================================
--- trunk/getfem/src/getfem_plasticity.cc       (original)
+++ trunk/getfem/src/getfem_plasticity.cc       Fri Apr  8 21:19:00 2016
@@ -932,31 +932,108 @@
     }
     
     // Value : u/|u|
-    void value(const arg_list &args, base_tensor &result) const
-    {
-      scalar_type no = gmm::vect_norm2(args[0]->as_vector());
-      if (no < 1E-15)
-        gmm::clear(result.as_vector());
-      else
-        gmm::copy(gmm::scaled(args[0]->as_vector(), scalar_type(1)/no),
-                  result.as_vector());
-    }
+    void value(const arg_list &args, base_tensor &result) const {
+      const base_tensor &t = *args[0];
+      scalar_type eps = 1E-25;
+      scalar_type no = 
::sqrt(gmm::vect_norm2_sqr(t.as_vector())+gmm::sqr(eps));
+      gmm::copy(gmm::scaled(t.as_vector(), scalar_type(1)/no),
+               result.as_vector());
+    }
+    // void value(const arg_list &args, base_tensor &result) const {
+    //   scalar_type no = gmm::vect_norm2(args[0]->as_vector());
+    //   if (no < 1E-15)
+    //     gmm::clear(result.as_vector());
+    //   else
+    //     gmm::copy(gmm::scaled(args[0]->as_vector(), scalar_type(1)/no),
+    //               result.as_vector());
+    // }
 
     // Derivative : (|u|^2 Id - u x u)/|u|^3
     void derivative(const arg_list &args, size_type,
                     base_tensor &result) const {
       const base_tensor &t = *args[0];
+      scalar_type eps = 1E-25;
+
       size_type N = t.size();
-      scalar_type no = gmm::vect_norm2(t.as_vector());
-
-      gmm::clear(result.as_vector());
-      if (no >= 1E-15) {
-        scalar_type no3 = no*no*no;
-        for (size_type i = 0; i < N; ++i) {
-          result[i*N+i] += scalar_type(1)/no;
-          for (size_type j = 0; j < N; ++j)
-            result[j*N+i] -= t[i]*t[j] / no3;
-        }
+      scalar_type no = 
::sqrt(gmm::vect_norm2_sqr(t.as_vector())+gmm::sqr(eps));
+      scalar_type no3 = no*no*no;
+
+      for (size_type i = 0; i < N; ++i) {
+       result[i*N+i] += scalar_type(1)/no;
+       for (size_type j = 0; j < N; ++j)
+         result[j*N+i] -= t[i]*t[j] / no3;
+      }
+    }
+    // void derivative(const arg_list &args, size_type,
+    //                 base_tensor &result) const {
+    //   const base_tensor &t = *args[0];
+    //   size_type N = t.size();
+    //   scalar_type no = gmm::vect_norm2(t.as_vector());
+
+    //   gmm::clear(result.as_vector());
+    //   if (no >= 1E-15) {
+    //     scalar_type no3 = no*no*no;
+    //     for (size_type i = 0; i < N; ++i) {
+    //       result[i*N+i] += scalar_type(1)/no;
+    //       for (size_type j = 0; j < N; ++j)
+    //         result[j*N+i] -= t[i]*t[j] / no3;
+    //     }
+    //   }
+    // }
+
+    // Second derivative : not implemented
+    void second_derivative(const arg_list &/*args*/, size_type, size_type,
+                           base_tensor &/*result*/) const {
+      GMM_ASSERT1(false, "Sorry, second derivative not implemented");
+    }
+  };
+
+
+  // Normalized_reg vector/matrix operator : Regularized Vector/matrix divided 
by its Frobenius norm
+  struct normalized_reg_operator : public ga_nonlinear_operator {
+    bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
+      if (args.size() != 2 || args[0]->sizes().size() > 2
+          || args[0]->sizes().size() < 1 || args[0]->size() != 1) return false;
+      if (args[0]->sizes().size() == 1)
+        ga_init_vector(sizes, args[0]->sizes()[0]);
+      else
+        ga_init_matrix(sizes, args[0]->sizes()[0], args[0]->sizes()[1]);
+      return true;
+    }
+    
+    // Value : u/(sqrt([u|^2+\eps^2))
+    void value(const arg_list &args, base_tensor &result) const {
+      const base_tensor &t = *args[0];
+      scalar_type eps = (*args[1])[0];
+      scalar_type no = 
::sqrt(gmm::vect_norm2_sqr(t.as_vector())+gmm::sqr(eps));
+      gmm::copy(gmm::scaled(t.as_vector(), scalar_type(1)/no),
+               result.as_vector());
+    }
+
+    // Derivative / u : ((|u|^2+eps^2) Id - u x u)/(|u|^2+eps^2)^(3/2)
+    // Derivative / eps : -eps*u/(|u|^2+eps^2)^(3/2)
+    void derivative(const arg_list &args, size_type nder,
+                    base_tensor &result) const {
+      const base_tensor &t = *args[0];
+      scalar_type eps = (*args[1])[0];
+
+      size_type N = t.size();
+      scalar_type no = 
::sqrt(gmm::vect_norm2_sqr(t.as_vector())+gmm::sqr(eps));
+      scalar_type no3 = no*no*no;
+
+      switch (nder) {
+      case 1:
+       for (size_type i = 0; i < N; ++i) {
+         result[i*N+i] += scalar_type(1)/no;
+         for (size_type j = 0; j < N; ++j)
+           result[j*N+i] -= t[i]*t[j] / no3;
+       }
+       break;
+       
+      case 2:
+       gmm::copy(gmm::scaled(t.as_vector(), -scalar_type(eps)/no3),
+                 result.as_vector());
+       break;
       }
     }
 
@@ -1068,6 +1145,8 @@
                                 std::make_shared<matrix_logarithm_operator>());
     PREDEF_OPERATORS.add_method("Normalized",
                                 std::make_shared<normalized_operator>());
+    PREDEF_OPERATORS.add_method("Normalized_reg",
+                                std::make_shared<normalized_reg_operator>());
     PREDEF_OPERATORS.add_method("Von_Mises_projection",
                                 
std::make_shared<Von_Mises_projection_operator>());
     return true;