[Getfem-commits] (no subject)

[Yves Renard]

branch: devel-yves-generic-assembly-modifs
commit 5b026af2a87c026692d4b8b4d8de12ceed546b32
Author: Yves Renard <address@hidden>
Date:   Wed May 9 20:34:52 2018 +0200

    Operator Diff(expression, variable, direction)
---
 .../source/project/libdesc_high_gen_assemb.rst     |   2 +-
 doc/sphinx/source/userdoc/gasm_high.rst            |  21 ++-
 interface/tests/python/check_asm.py                |  23 +++
 src/getfem_generic_assembly_semantic.cc            | 164 ++++++++++++++++++---
 4 files changed, 188 insertions(+), 22 deletions(-)

diff --git a/doc/sphinx/source/project/libdesc_high_gen_assemb.rst 
b/doc/sphinx/source/project/libdesc_high_gen_assemb.rst
index 57d75ab..194ef0c 100644
--- a/doc/sphinx/source/project/libdesc_high_gen_assemb.rst
+++ b/doc/sphinx/source/project/libdesc_high_gen_assemb.rst
@@ -30,7 +30,7 @@ Files
    :file:`getfem_generic_assembly_semantic.h` and 
:file:`getfem_generic_assembly_semantic.cc`, "Semantic analysis and enrichment 
of the syntax tree. Include some operations such as making the derivation of a 
tree with respect to a variable or computing the tree corresponding to the 
gradient of an expression."
    :file:`getfem_generic_assembly_workspace.cc`, "Methodes of the workspace 
object (defined in :file:`getfem_generic_assembly.h`)."
    :file:`getfem_generic_assembly_compile_and_exec.h` and 
:file:`getfem_generic_assembly_compile_and_exec.cc`, "Definition of the 
optimized instructions, compilation into a sequel of optimize instructions and 
execution of the instructions on Gauss point/interpolation points."   
-   :file:`getfem_generic_assembly_interpolation.cc`, "Interpolation operations 
and interpolate transformations"
+   :file:`getfem_generic_assembly_interpolation.cc`, "Interpolation operations 
and interpolate transformations."
 
 
 A few implementation details
diff --git a/doc/sphinx/source/userdoc/gasm_high.rst 
b/doc/sphinx/source/userdoc/gasm_high.rst
index ae221b4..7a74a24 100644
--- a/doc/sphinx/source/userdoc/gasm_high.rst
+++ b/doc/sphinx/source/userdoc/gasm_high.rst
@@ -50,7 +50,7 @@ A specific language has been developed to describe the weak 
formulation of bound
 
   - Explicit matrices: For instance ``[1,3;2,4]`` and ``[[1,2],[3,4]]`` denote 
the same 2x2 matrix. Each component can be an expression.
 
-  - Explicit fourth order tensors: Supplementary dimensions are separated with 
``,,`` and ``;;``. For instance 
``[[[[1,2,3],[1,2,3]],[[1,2,3],[1,2,3]]],[[[1,2,3],[1,2,3]],[[1,2,3],[1,2,3]]]]``
 is a 2x2x2x2 valid tensor.
+  - Explicit fourth order tensors: example of explicit 3x2x2x2 fourth order 
tensor in the nested format: 
``[[[[1,2,3],[1,2,3]],[[1,2,3],[1,2,3]]],[[[1,2,3],[1,2,3]],[[1,2,3],[1,2,3]]]]``.
 
   - ``X`` is the current coordinate on the real element, ``X(i)`` is its i-th 
component.
 
@@ -62,6 +62,8 @@ A specific language has been developed to describe the weak 
formulation of bound
 
   - ``Diff(expression, variable)``: The possibility to explicit differentiate 
an expression with respect to a variable (symbolic differentiation). 
 
+  - ``Diff(expression, variable, direction)``: computes the derivative of 
``expression`` with respect to ``variable`` in the direction ``direction``.
+
   - ``Grad(expression)``: When possible, symbolically derive the gradient of 
the given expression.
 
   - Possiblility of macro definition (in the model, the ga_workspace object or 
directly in the assembly string). The macros should be some valid expressions 
that are expanded inline at the lexical analysis phase (if they are used 
several times, the computation is automatically factorized at the compilation 
stage).
@@ -514,7 +516,7 @@ Similarly to explicit vectors, it is possible to define 
explicit matrices (i.e.
 Explicit tensors
 ----------------
 
-Explicit tensors of any order are permitted with the nested format. A tensor 
of order ``n`` is written as a succession of tensor of order ``n-1`` of equal 
dimensions and separated by a comma. For instance 
``[[[[1,2,3],[1,2,3]],[[1,2,3],[1,2,3]]],[[[1,2,3],[1,2,3]],[[1,2,3],[1,2,3]]]]``
 is a fourth order tensor. Another possibility is to use the syntax 
``Reshape([1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3], 2, 2, 2, 2)`` 
where the components have to be given in Fortran order.
+Explicit tensors of any order are permitted with the nested format. A tensor 
of order ``n`` is written as a succession of tensor of order ``n-1`` of equal 
dimensions and separated by a comma. For instance 
``[[[[1,2,3],[1,2,3]],[[1,2,3],[1,2,3]]],[[[1,2,3],[1,2,3]],[[1,2,3],[1,2,3]]]]``
 is a fourth order tensor. Another possibility is to use the syntax 
``Reshape([1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3], 3, 2, 2, 2)`` 
where the components have to be given in Fortran order.
 
 
 Access to tensor components
@@ -667,7 +669,20 @@ So that::
 
   Grad_u:Grad_test_u + Diff(u.u, u)
 
-is a valid expression.  
+is a valid expression. A third argument can be added to the ``Diff`` command 
to specify the direction::
+
+  Diff(expression, variable, direction)
+
+in that case, it replaces the ``Test_variable`` by the expression 
``direction`` which has to be of the same dimension as ``variable``. It 
computes the derivative of ``expression`` with respect to ``variable`` in the 
direction ``direction``.
+For instance::
+
+  Diff(u.u, u, v)
+
+will result in::
+
+  2*(u.v)
+
+if ``v`` is any valid expression of the same dimension than ``u``.
 
 Explicit Gradient
 -----------------
diff --git a/interface/tests/python/check_asm.py 
b/interface/tests/python/check_asm.py
index 0e82c44..da7d70f 100644
--- a/interface/tests/python/check_asm.py
+++ b/interface/tests/python/check_asm.py
@@ -201,3 +201,26 @@ str = "Grad(Norm(v))"; print 'Assembly string "%s" gives:' 
% str
 res = gf.asm('expression analysis', str,  mim, 0, md)
 if (res != "(Derivative_1_Norm(v).Grad_v)"):
   print "Bad gradient"; exit(1)
+
+str = "Diff((v*u).v, v)"; print 'Assembly string "%s" gives:' % str
+res = gf.asm('expression analysis', str,  mim, 1, md)
+
+str = "Diff((v*u).v, v, [0, 1, 3])"; print 'Assembly string "%s" gives:' % str
+res = gf.asm('expression analysis', str,  mim, 0, md)
+if (res != "((v.([0,1,3]*u))+((v*u).[0,1,3]))"):
+  print "Bad gradient"; exit(1)
+
+str = "Diff((w*u).Grad_u, u, 3)"; print 'Assembly string "%s" gives:' % str
+res = gf.asm('expression analysis', str,  mim, 0, md)
+if (res != "(Grad_u.(w*3))"):
+  print "Bad gradient"; exit(1)
+
+str = "Diff((w*u).Grad_u, u, X.w)"; print 'Assembly string "%s" gives:' % str
+res = gf.asm('expression analysis', str,  mim, 0, md)
+if (res != "((Grad_u.(w*(X.w)))+((w*u).((w.[[1,0],[0,1]])+(X.Grad_w))))"):
+  print "Bad gradient"; exit(1)
+
+str = "Diff((w*u).Grad_u, u, X(1))"; print 'Assembly string "%s" gives:' % str
+res = gf.asm('expression analysis', str,  mim, 0, md)
+if (res != "((Grad_u.(w*X(1)))+((w*u).[1,0]))"):
+  print "Bad gradient"; exit(1)
diff --git a/src/getfem_generic_assembly_semantic.cc 
b/src/getfem_generic_assembly_semantic.cc
index 4d4e212..1abff77 100644
--- a/src/getfem_generic_assembly_semantic.cc
+++ b/src/getfem_generic_assembly_semantic.cc
@@ -32,6 +32,21 @@ namespace getfem {
   extern bool predef_operators_plasticity_initialized;
   extern bool predef_operators_contact_initialized;
 
+  static void ga_node_derivation
+  (ga_tree &tree, const ga_workspace &workspace, const mesh &m,
+   pga_tree_node pnode, const std::string &varname,
+   const std::string &interpolatename, size_type order);
+
+  static void ga_node_grad(ga_tree &tree, const ga_workspace &workspace,
+                          const mesh &m, pga_tree_node pnode);
+  static bool ga_node_mark_tree_for_grad(pga_tree_node pnode);
+  static void ga_node_analysis(ga_tree &tree,
+                               const ga_workspace &workspace,
+                               pga_tree_node pnode, const mesh &me,
+                               size_type ref_elt_dim, bool eval_fixed_size,
+                               bool ignore_X, int option);
+
+
   bool ga_extract_variables(const pga_tree_node pnode,
                            const ga_workspace &workspace,
                            const mesh &m,
@@ -149,15 +164,107 @@ namespace getfem {
     return marked;
   }
 
-  static void ga_node_derivation
-  (ga_tree &tree, const ga_workspace &workspace, const mesh &m,
-   pga_tree_node pnode, const std::string &varname,
-   const std::string &interpolatename, size_type order);
+  static void ga_node_expand_expression_in_place_of_test
+  (ga_tree &tree, const ga_workspace &workspace,
+   pga_tree_node &pnode, const std::string &varname,
+   pga_tree_node pexpr, ga_tree &grad_expr, ga_tree &hess_expr,
+   size_type order, const mesh &me, size_type ref_elt_dim, bool 
eval_fixed_size,
+   bool ignore_X, int option) {
+    pga_tree_node parent = pnode->parent; 
+    for (size_type i = 0; i < pnode->children.size(); ++i)
+      ga_node_expand_expression_in_place_of_test
+       (tree, workspace, pnode->children[i], varname, pexpr, grad_expr,
+        hess_expr, order, me, ref_elt_dim, eval_fixed_size, ignore_X, option);
+    const std::string &name = pnode->name;
+    size_type loc_order = pnode->test_function_type;
+
+    if (loc_order == order && !(name.compare(varname))) {
+      bool need_grad = (pnode->node_type == GA_NODE_GRAD_TEST ||
+                       pnode->node_type == GA_NODE_DIVERG_TEST ||
+                       pnode->node_type == GA_NODE_HESS_TEST);
+      bool need_hess = (pnode->node_type == GA_NODE_HESS_TEST);
+
+      if (need_grad && grad_expr.root == nullptr) {
+       tree.copy_node(pexpr, nullptr, grad_expr.root);
+       if (ga_node_mark_tree_for_grad(grad_expr.root)) {
+         ga_node_grad(grad_expr, workspace, me, grad_expr.root);
+         ga_node_analysis(grad_expr, workspace, grad_expr.root, me,
+                          ref_elt_dim, eval_fixed_size, ignore_X, option);
+       } else {
+         bgeot::multi_index mi = grad_expr.root->t.sizes();
+         mi.push_back(me.dim());
+         grad_expr.root->t.adjust_sizes(mi);
+         grad_expr.root->node_type = GA_NODE_ZERO;
+         gmm::clear(grad_expr.root->tensor().as_vector());
+         grad_expr.clear_children(grad_expr.root);
+       }
+      }
+
+      if (need_hess && hess_expr.root == nullptr) {
+       tree.copy_node(grad_expr.root, nullptr, hess_expr.root);
+       if (ga_node_mark_tree_for_grad(hess_expr.root)) {
+         ga_node_grad(hess_expr, workspace, me, hess_expr.root);
+         ga_node_analysis(hess_expr, workspace, hess_expr.root, me,
+                          ref_elt_dim, eval_fixed_size, ignore_X, option);
+       } else {
+         bgeot::multi_index mi = hess_expr.root->t.sizes();
+         mi.push_back(me.dim());
+         hess_expr.root->t.adjust_sizes(mi);
+         hess_expr.root->node_type = GA_NODE_ZERO;
+         gmm::clear(hess_expr.root->tensor().as_vector());
+         hess_expr.clear_children(grad_expr.root);
+       }
+      }
+      switch(pnode->node_type) {
+      case GA_NODE_VAL_TEST:
+       delete pnode; pnode = nullptr;
+       tree.copy_node(pexpr, parent, pnode);
+       break;
+      case GA_NODE_GRAD_TEST:
+       delete pnode; pnode = nullptr;
+       tree.copy_node(grad_expr.root, parent, pnode);
+       break;
+      case GA_NODE_HESS_TEST:
+       delete pnode; pnode = nullptr;
+       tree.copy_node(hess_expr.root, parent, pnode);
+       break;
+      case GA_NODE_DIVERG_TEST:
+       {
+         delete pnode; pnode = nullptr;
+         tree.copy_node(grad_expr.root, parent, pnode);
+         tree.insert_node(pnode, GA_NODE_OP);
+         pnode->parent->op_type = GA_COLON;
+         tree.add_child(pnode->parent, GA_NODE_PARAMS);
+         pga_tree_node pid = pnode->parent->children[1];
+         tree.add_child(pid); tree.add_child(pid);
+         pid->children[0]->node_type = GA_NODE_NAME;
+         pid->children[0]->name = "Id";
+         pid->children[1]->node_type = GA_NODE_CONSTANT;
+         pid->children[1]->init_scalar_tensor(me.dim());
+       }
+       break;
+      case GA_NODE_INTERPOLATE_VAL_TEST: case GA_NODE_INTERPOLATE_GRAD_TEST:
+      case GA_NODE_INTERPOLATE_HESS_TEST: case GA_NODE_INTERPOLATE_DIVERG_TEST:
+       GMM_ASSERT1(false,
+                   "Sorry, directional derivative do not work for the moment "
+                   "with interpolate transformations. Future work.");
+      case GA_NODE_ELEMENTARY_VAL_TEST: case GA_NODE_ELEMENTARY_GRAD_TEST:
+      case GA_NODE_ELEMENTARY_HESS_TEST: case GA_NODE_ELEMENTARY_DIVERG_TEST:
+       GMM_ASSERT1(false,
+                   "Sorry, directional derivative do not work for the moment "
+                   "with elementary transformations. Future work.");
+      case GA_NODE_XFEM_PLUS_VAL_TEST: case GA_NODE_XFEM_PLUS_GRAD_TEST:
+      case GA_NODE_XFEM_PLUS_HESS_TEST: case GA_NODE_XFEM_PLUS_DIVERG_TEST:
+      case GA_NODE_XFEM_MINUS_VAL_TEST: case GA_NODE_XFEM_MINUS_GRAD_TEST:
+      case GA_NODE_XFEM_MINUS_HESS_TEST: case GA_NODE_XFEM_MINUS_DIVERG_TEST:
+       GMM_ASSERT1(false,
+                   "Sorry, directional derivative do not work for the moment "
+                   "with Xfem_plus and Xfem_minus operations. Future work.");
+      default: break;
+      }
+    }
+  }
 
-  static void ga_node_grad(ga_tree &tree, const ga_workspace &workspace,
-                          const mesh &m, pga_tree_node pnode);
-  static bool ga_node_mark_tree_for_grad(pga_tree_node pnode);
-  
   //=========================================================================
   // Some hash code functions for node identification
   //=========================================================================
@@ -1696,7 +1803,7 @@ namespace getfem {
                             ref_elt_dim, eval_fixed_size, ignore_X, option);
            child1 = pnode->children[1];
          } else {
-           mi = child1->t.sizes(); mi.push_back(2);
+           mi = child1->t.sizes(); mi.push_back(me.dim());
            child1->t.adjust_sizes(mi);
            child1->node_type = GA_NODE_ZERO;
            gmm::clear(child1->tensor().as_vector());
@@ -1706,7 +1813,7 @@ namespace getfem {
          pnode = child1;
        } else if (child0->name.compare("Diff") == 0) {
          
-         if (pnode->children.size() != 3)
+         if (pnode->children.size() != 3 && pnode->children.size() != 4)
            ga_throw_error(pnode->expr, child0->pos,
                           "Bad number of parameters for Diff operator");
          pga_tree_node child2 = pnode->children[2];
@@ -1733,6 +1840,15 @@ namespace getfem {
            gmm::clear(child1->tensor().as_vector());
            tree.clear_children(child1);
          }
+         if (pnode->children.size() == 4) {
+           ga_tree grad_expr, hess_expr;
+           ga_node_expand_expression_in_place_of_test
+             (tree, workspace, pnode->children[1], vardiff, pnode->children[3],
+              grad_expr, hess_expr, order+1, me, ref_elt_dim, eval_fixed_size,
+              ignore_X, option);
+            ga_node_analysis(tree, workspace, pnode->children[1], me,
+                            ref_elt_dim, eval_fixed_size, ignore_X, option);
+         }
          tree.replace_node_by_child(pnode, 1);
          pnode = child1;
        } else
@@ -3724,10 +3840,16 @@ namespace getfem {
            child1->node_type = GA_NODE_CONSTANT;
          } else if (pnode->node_type == GA_NODE_INTERPOLATE_X) {
            pnode->node_type = GA_NODE_CONSTANT;
-           pnode->init_matrix_tensor(trans_dim, trans_dim);
-           gmm::clear(pnode->tensor().as_vector());
-           for (size_type i = 0; i < trans_dim; ++i)
-             pnode->tensor()(i,i) = scalar_type(1);
+           if (pnode->nbc1) {
+             pnode->init_vector_tensor(trans_dim);
+             gmm::clear(pnode->tensor().as_vector());
+             pnode->tensor()[pnode->nbc1-1] = scalar_type(1);
+           } else {
+             pnode->init_matrix_tensor(trans_dim, trans_dim);
+             gmm::clear(pnode->tensor().as_vector());
+             for (size_type i = 0; i < trans_dim; ++i)
+               pnode->tensor()(i,i) = scalar_type(1);
+           }
          }
          
          tree.clear_node_rec(pnode->parent->children[1]);
@@ -3744,10 +3866,16 @@ namespace getfem {
       
     case GA_NODE_X:
       pnode->node_type = GA_NODE_CONSTANT;
-      pnode->init_matrix_tensor(meshdim, meshdim);
-      gmm::clear(pnode->tensor().as_vector());
-      for (size_type i = 0; i < meshdim; ++i)
-       pnode->tensor()(i,i) = scalar_type(1);
+      if (pnode->nbc1) {
+       pnode->init_vector_tensor(meshdim);
+       gmm::clear(pnode->tensor().as_vector());
+        pnode->tensor()[pnode->nbc1-1] = scalar_type(1);
+      } else {
+       pnode->init_matrix_tensor(meshdim, meshdim);
+       gmm::clear(pnode->tensor().as_vector());
+       for (size_type i = 0; i < meshdim; ++i)
+         pnode->tensor()(i,i) = scalar_type(1);
+      }
       break;
 
     case GA_NODE_NORMAL: