[Getfem-commits] (no subject)

[Yves Renard]

branch: devel-yves-generic-assembly-modifs
commit c62857285b7da6ce66a3ceccd70fd35b1d8c777a
Author: Yves Renard <address@hidden>
Date:   Wed May 9 17:15:17 2018 +0200

    Grad(expression) working
---
 .../source/project/libdesc_high_gen_assemb.rst     | 12 +++-
 doc/sphinx/source/userdoc/gasm_high.rst            | 12 ++--
 interface/tests/python/check_asm.py                | 33 +++++-----
 src/getfem_generic_assembly_semantic.cc            | 76 ++++++++++++++--------
 4 files changed, 79 insertions(+), 54 deletions(-)

diff --git a/doc/sphinx/source/project/libdesc_high_gen_assemb.rst 
b/doc/sphinx/source/project/libdesc_high_gen_assemb.rst
index 594d53f..57d75ab 100644
--- a/doc/sphinx/source/project/libdesc_high_gen_assemb.rst
+++ b/doc/sphinx/source/project/libdesc_high_gen_assemb.rst
@@ -24,11 +24,15 @@ Files
    :header: "File(s)", "Description"
    :widths: 8, 15
 
-   :file:`getfem_generic_assembly.h` and :file:`getfem_generic_assembly.cc`, 
"In order not to export all the internal of the generic assembly, all is 
implemented in a single file"
+   :file:`getfem_generic_assembly.h`, "Main header for exported definitions. 
Only this header has to be included to use the generic assembly. Other headers 
of the module are for internal use only."
+   :file:`getfem_generic_assembly_tree.h` and 
:file:`getfem_generic_assembly_tree.cc`, "Definition of the tree structure and 
basic operations on it, including reading an assembly string and transform it 
in a syntax tree and make the invert transformation of a tree into a string."
+   :file:`getfem_generic_assembly_fonction_and_operators.h` and 
:file:`getfem_generic_assembly_fonction_and_operators.cc`, "Definition of 
redefined function and nonlinear operator of the assembly language."
+   :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"
 
 
-For the moment, having a single (giant) file simplify the maintenance.
-
 A few implementation details
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -42,6 +46,8 @@ The assembly string is transformed in an assembly tree by a 
set of function in :
 
  - Symbolic (automatic) differentiation of an assembly tree by 
``ga_derivative(...)``
 
+ - Symbolic (automatic) gradient computation of an assembly tree by 
``ga_gradient(...)``
+
  - Compilation in a sequence of instructions with optimizations by 
``ga_compile(...)``.
 
  - Execution of the sequence of instructions and assembly by ``ga_exec(...)``.
diff --git a/doc/sphinx/source/userdoc/gasm_high.rst 
b/doc/sphinx/source/userdoc/gasm_high.rst
index 2148d63..ae221b4 100644
--- a/doc/sphinx/source/userdoc/gasm_high.rst
+++ b/doc/sphinx/source/userdoc/gasm_high.rst
@@ -502,23 +502,19 @@ Parentheses can be used in a standard way to change the 
operation order. If no p
 Explicit vectors
 ----------------
 
-The assembly language allows to define explicit vectors (i.e. order 1 tensors) 
with the notation ``[a;b;c;d;e]``, i.e. an arbitrary number of components 
separated by a semicolon, the whole vector beginning with a right bracket and 
ended by a left bracket. The components can be some numeric constants, some 
valid expressions and may also contain test functions. In the latter case, the 
vector has to be homogeneous with respect to the test functions. This means 
that a construction of the typ [...]
+The assembly language allows to define explicit vectors (i.e. order 1 tensors) 
with the notation ``[a,b,c,d,e]``, i.e. an arbitrary number of components 
separated by a comma (note the separation with a semicolon ``[a;b;c;d;e]`` is 
also permitted), the whole vector beginning with a right bracket and ended by a 
left bracket. The components can be some numeric constants, some valid 
expressions and may also contain test functions. In the latter case, the vector 
has to be homogeneous with res [...]
 
 
 Explicit matrices
 -----------------
 
-Similarly to explicit vectors, it is possible to define explicit matrices 
(i.e. order 2 tensors) with the notation ``[a,c;b,d]``,  i.e. an arbitrary 
number of lines separated by a semicolon, each line having the same number of 
components separated by a comma. Alternatively the nested format 
``[[a,b],[c,d]]`` provides an equivalent result. For instance 
``[11,12,13;21,22,23]`` and ``[[11,21],[12,22],[13,23]]`` both represent the 
same 2x3 matrix. The components can be some numeric constants [...]
+Similarly to explicit vectors, it is possible to define explicit matrices 
(i.e. order 2 tensors) with the notation ``[[a,b],[c,d]]``, i.e. an arbitrary 
number of columns vectors separated by a comma (the syntax ``[a,c;b,d]`` of 
lines separated by a semicolon is also permitted). For instance 
``[[11,21],[12,22],[13,23]]`` and ``[11,12,13;21,22,23]`` both represent the 
same 2x3 matrix. The components can be some numeric constants, some valid 
expressions and may also contain test functions.
 
 
 Explicit tensors
----------------------------
-
-Explicit order four tensors are also allowed. To this aim, the two 
supplementary dimensions compared to matrices are separated by  ``,,`` and 
``;;``. For instance ``[1,1;1,1,,1,1;1,1;;2,2;2,2,,2,2;2,2;;3,3;3,3,,3,3;3,3]`` 
is a valid 3x2x2x2 tensor. Note that constant fourth order tensors can also be 
obtained by the tensor product of two constant matrices or by the Reshape 
instruction. 
-
-The nested format can also be used for defining order four tensors. The 
previous example is equivalent to writting 
``[[[[1,2,3],[1,2,3]],[[1,2,3],[1,2,3]]],[[[1,2,3],[1,2,3]],[[1,2,3],[1,2,3]]]]``.
 The nested format also allows the definition of order three tensors.
+----------------
 
-Explicit order five or six tensors are not directly supported by the assembly 
language. However, they can be easily obtained via the Reshape instruction.
+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.
 
 
 Access to tensor components
diff --git a/interface/tests/python/check_asm.py 
b/interface/tests/python/check_asm.py
index 4c88cea..0e82c44 100644
--- a/interface/tests/python/check_asm.py
+++ b/interface/tests/python/check_asm.py
@@ -19,7 +19,7 @@
 # Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
 #
 ############################################################################
-"""  test high generic assembly language.
+"""  Test of the high generic assembly language.
 
   This program is used to check that Python-GetFEM interface, and more
   generally GetFEM are working. It focuses on testing some operations
@@ -158,14 +158,11 @@ print 'Assembly string "Grad([[u,2,u],[u,1,u]])" gives:'
 res = gf.asm('expression analysis', "Grad([[u,2,u],[u,1,u]])",  mim, 0, md)
 if (res != 
"([[[Grad_u(1),0,Grad_u(1)],[Grad_u(1),0,Grad_u(1)]],[[Grad_u(2),0,Grad_u(2)],[Grad_u(2),0,Grad_u(2)]]])"):
   print "Bad gradient"; exit(1)
-  
-print 'Assembly string "Grad([u,u])" gives:'
-res = gf.asm('expression analysis', "Grad([u,u])",  mim, 0, md)
 
 print 'Assembly string "Grad([u;u])" gives:'
 res = gf.asm('expression analysis', "Grad([u,u])",  mim, 0, md)
-
-
+if (res != "([[Grad_u(1),Grad_u(1)],[Grad_u(2),Grad_u(2)]])"):
+  print "Bad gradient"; exit(1)
 
 print 'Assembly string "Grad(Reshape(Grad_w, 1, 4))" gives:'
 res = gf.asm('expression analysis', "Grad(Reshape(Grad_w, 1, 4))",  mim, 0, md)
@@ -186,15 +183,21 @@ if (res != "(Contract(Hess_w, 1, 2, Grad_w, 1, 
2)+Contract(Grad_w, 1, 2, Hess_w,
   print "Bad gradient"; exit(1)
 
 
-str = "[1;2;3]"; print 'Assembly string "%s" gives:' % str
+str = "Grad(sin(u))"; print 'Assembly string "%s" gives:' % str
 res = gf.asm('expression analysis', str,  mim, 0, md)
+if (res != "((cos(u)@[1,1]).*Grad_u)"): print "Bad gradient"; exit(1)
 
-str = "[1,2,3]"; print 'Assembly string "%s" gives:' % str
+str = "Grad(cos(Grad_u))"; print 'Assembly string "%s" gives:' % str
 res = gf.asm('expression analysis', str,  mim, 0, md)
-
-
-str = "[u;u;u].[u;u;u]"; print 'Assembly string "%s" gives:' % str
-res = gf.asm('expression analysis', str,  mim, 2, md)
-
-
-
+if (res != "((DER_PDFUNC_COS(Grad_u)@[1,1]).*Hess_u)"):
+  print "Bad gradient"; exit(1)
+  
+str = "Grad(min(v, 2*v))"; print 'Assembly string "%s" gives:' % str
+res = gf.asm('expression analysis', str,  mim, 0, md)
+if (res != "(((DER_PDFUNC2_MAX(v, (2*v))@[1,1]).*Grad_v)+((DER_PDFUNC1_MAX(v, 
(2*v))@[1,1]).*(2*Grad_v)))"):
+  print "Bad gradient"; exit(1)
+  
+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)
diff --git a/src/getfem_generic_assembly_semantic.cc 
b/src/getfem_generic_assembly_semantic.cc
index 474a6a8..4d4e212 100644
--- a/src/getfem_generic_assembly_semantic.cc
+++ b/src/getfem_generic_assembly_semantic.cc
@@ -1018,6 +1018,7 @@ namespace getfem {
         {
           size_type s0 = dim0 == 0 ? 1 : size0.back();
           size_type s1 = dim1 == 0 ? 1 : size1[size1.size()-dim1];
+        
           if (s0 != s1) ga_throw_error(pnode->expr, pnode->pos, "Dot product "
                                        "of expressions of different sizes ("
                                        << s0 << " != " << s1 << ").");
@@ -3211,7 +3212,8 @@ namespace getfem {
           } else {
             tree.duplicate_with_addition(pnode);
             if ((pnode->op_type == GA_COLON && child0->tensor_order() == 2) ||
-                (pnode->op_type == GA_DOT && child0->tensor_order() == 1) ||
+                (pnode->op_type == GA_DOT && child0->tensor_order() == 1 &&
+                child1->tensor_order() == 1) ||
                 pnode->op_type == GA_DOTMULT ||
                 (child0->tensor_proper_size()== 1 &&
                  child1->tensor_proper_size()== 1))
@@ -3382,9 +3384,12 @@ namespace getfem {
           }
         } else {
           pga_tree_node child2 = pnode->children[2];
-
-          if (child1->marked && child2->marked)
+         pga_tree_node pg2 = pnode;
+         
+          if (child1->marked && child2->marked) {
             tree.duplicate_with_addition(pnode);
+           pg2 = pnode->parent->children[1];
+         }
 
           if (child1->marked) {
             switch (F.dtype()) {
@@ -3419,8 +3424,7 @@ namespace getfem {
                                varname, interpolatename, order);
           }
           if (child2->marked) {
-            if (child1->marked && child2->marked)
-              pnode = pnode->parent->parent->children[1];
+           pnode = pg2;
             child0 = pnode->children[0]; child1 = pnode->children[1];
             child2 = pnode->children[2];
 
@@ -3986,7 +3990,8 @@ namespace getfem {
          
          if (pg1) {
            if ((pg1->op_type == GA_COLON && child0->tensor_order() == 2) ||
-               (pg1->op_type == GA_DOT && child0->tensor_order() == 1) ||
+               (pg1->op_type == GA_DOT && child0->tensor_order() == 1 &&
+                child1->tensor_order() == 1) ||
                pg1->op_type == GA_DOTMULT ||
                (child0->tensor_proper_size()== 1 ||
                 child1->tensor_proper_size()== 1)) {
@@ -4250,7 +4255,6 @@ namespace getfem {
            ga_node_grad(tree, workspace, m, pg2->children[ch2]);
          }
        }
-#ifdef continue_here
 
       } else if (child0->node_type == GA_NODE_PREDEF_FUNC) {
         std::string name = child0->name;
@@ -4309,18 +4313,27 @@ namespace getfem {
             pga_tree_node pnode_op = pnode->parent;
             if (child1->tensor_order() == 0)
               pnode_op->op_type = GA_MULT;
-            else
-              pnode_op->op_type = GA_DOTMULT;
+            else {
+             pnode_op->op_type = GA_DOTMULT;
+             tree.insert_node(pnode, GA_NODE_OP);
+             pnode->parent->op_type = GA_TMULT;
+             tree.add_child(pnode->parent, GA_NODE_CONSTANT);
+             pnode->parent->children[1]->init_vector_tensor(m.dim());
+             gmm::fill(pnode->parent->children[1]->tensor().as_vector(),
+                       scalar_type(1));
+           }
             pnode_op->children.resize(2, nullptr);
             tree.copy_node(child1, pnode_op, pnode_op->children[1]);
-            ga_node_grad(tree, workspace, m, pnode_op->children[1],
-                               varname, interpolatename, order);
+            ga_node_grad(tree, workspace, m, pnode_op->children[1]);
           }
         } else {
           pga_tree_node child2 = pnode->children[2];
-
-          if (child1->marked && child2->marked)
+         pga_tree_node pg2 = pnode;
+         
+          if (child1->marked && child2->marked) {
             tree.duplicate_with_addition(pnode);
+           pg2 = pnode->parent->children[1];
+         }
 
           if (child1->marked) {
             switch (F.dtype()) {
@@ -4347,16 +4360,21 @@ namespace getfem {
             pga_tree_node pnode_op = pnode->parent;
             if (child1->tensor_order() == 0)
               pnode_op->op_type = GA_MULT;
-            else
-              pnode_op->op_type = GA_DOTMULT;
+            else {
+             pnode_op->op_type = GA_DOTMULT;
+             tree.insert_node(pnode, GA_NODE_OP);
+             pnode->parent->op_type = GA_TMULT;
+             tree.add_child(pnode->parent, GA_NODE_CONSTANT);
+             pnode->parent->children[1]->init_vector_tensor(m.dim());
+             gmm::fill(pnode->parent->children[1]->tensor().as_vector(),
+                       scalar_type(1));
+           }
             pnode_op->children.resize(2, nullptr);
             tree.copy_node(child1, pnode_op, pnode_op->children[1]);
-            ga_node_grad(tree, workspace, m, pnode_op->children[1],
-                               varname, interpolatename, order);
+            ga_node_grad(tree, workspace, m, pnode_op->children[1]);
           }
           if (child2->marked) {
-            if (child1->marked && child2->marked)
-              pnode = pnode->parent->parent->children[1];
+           pnode = pg2;
             child0 = pnode->children[0]; child1 = pnode->children[1];
             child2 = pnode->children[2];
 
@@ -4384,12 +4402,18 @@ namespace getfem {
             pga_tree_node pnode_op = pnode->parent;
             if (child2->tensor_order() == 0)
               pnode_op->op_type = GA_MULT;
-            else
-              pnode_op->op_type = GA_DOTMULT;
+            else {
+             pnode_op->op_type = GA_DOTMULT;
+             tree.insert_node(pnode, GA_NODE_OP);
+             pnode->parent->op_type = GA_TMULT;
+             tree.add_child(pnode->parent, GA_NODE_CONSTANT);
+             pnode->parent->children[1]->init_vector_tensor(m.dim());
+             gmm::fill(pnode->parent->children[1]->tensor().as_vector(),
+                       scalar_type(1));
+           }
             pnode_op->children.resize(2, nullptr);
             tree.copy_node(child2, pnode_op, pnode_op->children[1]);
-            ga_node_grad(tree, workspace, m, pnode_op->children[1],
-                               varname, interpolatename, order);
+            ga_node_grad(tree, workspace, m, pnode_op->children[1]);
           }
         }
       } else if (child0->node_type == GA_NODE_SPEC_FUNC) {
@@ -4437,8 +4461,7 @@ namespace getfem {
             pnode_op->children.resize(2, nullptr);
             tree.copy_node(pnode->children[i], pnode_op,
                            pnode_op->children[1]);
-            ga_node_grad(tree, workspace, m, pnode_op->children[1],
-                               varname, interpolatename, order);
+            ga_node_grad(tree, workspace, m, pnode_op->children[1]);
 
             if (pnode2->children[0]->name.compare("Norm_sqr") == 0
                 && pnode2->children[0]->der1 == 1) {
@@ -4447,11 +4470,8 @@ namespace getfem {
               pnode2->children[0]->node_type = GA_NODE_CONSTANT;
               pnode2->children[0]->init_scalar_tensor(scalar_type(2));
             }
-
-
           }
         }
-#endif
       } else {
         ga_node_grad(tree, workspace, m, child0);
        tree.add_child(pnode, GA_NODE_ALLINDICES);