[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Yves Renard |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Thu, 10 May 2018 03:20:43 -0400 (EDT) |
branch: devel-yves-generic-assembly-modifs
commit aa440849c42811b4dc3977d89e49a82ed31087f1
Author: Yves Renard <address@hidden>
Date: Wed May 9 08:48:48 2018 +0200
Change in the storage of explicit tensor : natural order. Accepts now
higher order tensors with the nested format
---
doc/sphinx/source/userdoc/gasm_high.rst | 2 +-
interface/tests/python/check_asm.py | 30 ++-
src/getfem/getfem_generic_assembly_tree.h | 5 +-
src/getfem_generic_assembly_compile_and_exec.cc | 39 +---
src/getfem_generic_assembly_semantic.cc | 123 +++++-----
src/getfem_generic_assembly_tree.cc | 285 +++++++++++++-----------
6 files changed, 249 insertions(+), 235 deletions(-)
diff --git a/doc/sphinx/source/userdoc/gasm_high.rst
b/doc/sphinx/source/userdoc/gasm_high.rst
index 5710cce..2148d63 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,1;1,2,,1,1;1,2;;1,1;1,2,,1,1;1,2]`` is a
2x2x2x2 valid tensor.
+ - 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.
- ``X`` is the current coordinate on the real element, ``X(i)`` is its i-th
component.
diff --git a/interface/tests/python/check_asm.py
b/interface/tests/python/check_asm.py
index 85cd615..d3f7c1f 100644
--- a/interface/tests/python/check_asm.py
+++ b/interface/tests/python/check_asm.py
@@ -56,6 +56,7 @@ md.set_variable('v', V)
md.add_fem_variable('w', mfw)
md.set_variable('w', W)
+
# Simple test on the integral of u
result = gf.asm('generic', mim, 0, "u", -1, md)
if (abs(result-1) > 1e-8) : print "Bad value"; exit(1)
@@ -140,7 +141,7 @@ if (res !=
print 'Assembly string "Grad(Grad_w./Grad_w)" gives:'
res = gf.asm('expression analysis', "Grad(Grad_w./Grad_w)", mim, 0, md)
if (res !=
- "((Hess_w./(address@hidden; 1]))-(((Grad_w./sqr(Grad_w))@[1;
1]).*Hess_w))"):
+ "((Hess_w./(address@hidden;1]))-(((Grad_w./sqr(Grad_w))@[1;1]).*Hess_w))"):
print "Bad gradient"; exit(1)
print 'Assembly string "Grad(Grad_w/u)" gives:'
@@ -148,14 +149,21 @@ res = gf.asm('expression analysis', "Grad(Grad_w/u)",
mim, 0, md)
if (res != "((Hess_w/u)-((Grad_w/sqr(u))@Grad_u))"):
print "Bad gradient"; exit(1)
-# To be controlled after fixing C_MATRIX
print 'Assembly string "Grad([u,u; 2,1; u,u])" gives:'
res = gf.asm('expression analysis', "Grad([u,u; 2,1; u,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)
-
-# To be controlled after fixing C_MATRIX
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)
@@ -176,3 +184,17 @@ print 'Assembly string "Grad(Contract(Grad_w, 1, 2,
Grad_w, 1, 2))" gives:'
res = gf.asm('expression analysis', "Grad(Contract(Grad_w, 1, 2, Grad_w, 1,
2))", mim, 0, md)
if (res != "(Contract(Hess_w, 1, 2, Grad_w, 1, 2)+Contract(Grad_w, 1, 2,
Hess_w, 1, 2))"):
print "Bad gradient"; exit(1)
+
+
+str = "[1;2;3]"; print 'Assembly string "%s" gives:' % str
+res = gf.asm('expression analysis', str, mim, 0, md)
+
+str = "[1,2,3]"; 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)
+
+
+
diff --git a/src/getfem/getfem_generic_assembly_tree.h
b/src/getfem/getfem_generic_assembly_tree.h
index bd0dff1..967cd0d 100644
--- a/src/getfem/getfem_generic_assembly_tree.h
+++ b/src/getfem/getfem_generic_assembly_tree.h
@@ -300,7 +300,9 @@ namespace getfem {
std::string interpolate_name_test1, interpolate_name_test2; // name
// of interpolation transformation if any
size_type qdim1, qdim2; // Qdims when test_function_type > 0.
- size_type nbc1, nbc2, nbc3; // For explicit matrices, X and macros.
+ size_type nbc1, nbc2, nbc3; // For X (nbc1=coordinate number),
+ // macros (nbc1=param number, nbc2,nbc3
type))
+ // and C_MATRIX (nbc1=order).
size_type pos; // Position of the first character in string
pstring expr; // Original string, for error messages.
std::string name; // variable/constant/function/operator name
@@ -415,7 +417,6 @@ namespace getfem {
void add_sub_tree(ga_tree &sub_tree);
void add_params(size_type pos, pstring expr);
void add_matrix(size_type pos, pstring expr);
- void zip_matrix(const pga_tree_node source_node);
void add_op(GA_TOKEN_TYPE op_type, size_type pos, pstring expr);
void clear_node_rec(pga_tree_node pnode);
void clear_node(pga_tree_node pnode);
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index c8f825c..ae858f8 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -5952,47 +5952,16 @@ namespace getfem {
case GA_NODE_C_MATRIX:
{
- size_type nbc1 = pnode->nbc1, nbc2 = pnode->nbc2, nbc3 = pnode->nbc3;
- size_type nbl = pnode->children.size() / (nbc1*nbc2*nbc3);
if (pnode->test_function_type) {
std::vector<const base_tensor *> components(pnode->children.size());
-
- if (nbc1 == 1 && nbc2 == 1 && nbc3 == 1) {
- for (size_type i = 0; i < pnode->children.size(); ++i)
- components[i] = &(pnode->children[i]->tensor());
- } else if (nbc2 == 1 && nbc3 == 1) {
- for (size_type i = 0; i < nbl; ++i)
- for (size_type j = 0; j < nbc1; ++j)
- components[i+j*nbl] = &(pnode->children[i*nbc1+j]->tensor());
- } else {
- size_type n = 0;
- for (size_type i = 0; i < nbl; ++i)
- for (size_type j = 0; j < nbc3; ++j)
- for (size_type k = 0; k < nbc2; ++k)
- for (size_type l = 0; l < nbc1; ++l)
- components[i+j*nbl+k*nbl*nbc3+l*nbc2*nbc3*nbl]
- = &(pnode->children[n++]->tensor());
- }
+ for (size_type i = 0; i < pnode->children.size(); ++i)
+ components[i] = &(pnode->children[i]->tensor());
pgai = std::make_shared<ga_instruction_c_matrix_with_tests>
(pnode->tensor(), components);
} else {
std::vector<scalar_type *> components(pnode->children.size());
- if (nbc1 == 1 && nbc2 == 1 && nbc3 == 1) {
- for (size_type i = 0; i < pnode->children.size(); ++i)
- components[i] = &(pnode->children[i]->tensor()[0]);
- } else if (nbc2 == 1 && nbc3 == 1) {
- for (size_type i = 0; i < nbl; ++i)
- for (size_type j = 0; j < nbc1; ++j)
- components[i+j*nbl] =
&(pnode->children[i*nbc1+j]->tensor()[0]);
- } else {
- size_type n = 0;
- for (size_type i = 0; i < nbl; ++i)
- for (size_type j = 0; j < nbc3; ++j)
- for (size_type k = 0; k < nbc2; ++k)
- for (size_type l = 0; l < nbc1; ++l)
- components[i+j*nbl+k*nbl*nbc3+l*nbc2*nbc3*nbl]
- = &(pnode->children[n++]->tensor()[0]);
- }
+ for (size_type i = 0; i < pnode->children.size(); ++i)
+ components[i] = &(pnode->children[i]->tensor()[0]);
pgai = std::make_shared<ga_instruction_simple_c_matrix>
(pnode->tensor(), components);
}
diff --git a/src/getfem_generic_assembly_semantic.cc
b/src/getfem_generic_assembly_semantic.cc
index 983a70f..9c6f154 100644
--- a/src/getfem_generic_assembly_semantic.cc
+++ b/src/getfem_generic_assembly_semantic.cc
@@ -1348,9 +1348,9 @@ namespace getfem {
"components should be scalar valued.");
}
- size_type nbc1 = pnode->nbc1, nbc2 = pnode->nbc2, nbc3 = pnode->nbc3;
- size_type nbl = pnode->children.size() / (nbc1*nbc2*nbc3);
- if (all_cte) pnode->node_type = GA_NODE_CONSTANT;
+ GMM_ASSERT1(pnode->children.size() == pnode->tensor_proper_size(),
+ "Internal error");
+
pnode->test_function_type = 0;
for (size_type i = 0; i < pnode->children.size(); ++i) {
if (pnode->children[i]->test_function_type) {
@@ -1378,42 +1378,39 @@ namespace getfem {
}
}
}
- mi.resize(0);
- if (pnode->test_function_type) mi.push_back(2);
- if (pnode->test_function_type >= 3) mi.push_back(2);
- if (nbc1 == 1 && nbc2 == 1 && nbc3 == 1 && nbl == 1) {
- pnode->t.adjust_sizes(mi);
- if (all_cte) pnode->tensor()[0] = child0->tensor()[0];
- } else {
- mi.push_back(nbl);
- if (nbc3 != 1) mi.push_back(nbc3);
- if (nbc2 != 1) mi.push_back(nbc2);
- if (nbc1 != 1) mi.push_back(nbc1);
- pnode->t.adjust_sizes(mi);
- if (all_cte) {
- size_type n = 0;
- if (nbc1 == 1 && nbc2 == 1 && nbc3 == 1)
- for (size_type i = 0; i < nbl; ++i)
- pnode->tensor()[i] = pnode->children[i]->tensor()[0];
- else if (nbc2 == 1 && nbc3 == 1)
- for (size_type i = 0; i < nbl; ++i)
- for (size_type j = 0; j < nbc1; ++j)
- pnode->tensor()(i,j) = pnode->children[n++]->tensor()[0];
- else if (nbc3 == 1)
- for (size_type i = 0; i < nbl; ++i)
- for (size_type j = 0; j < nbc2; ++j)
- for (size_type k = 0; k < nbc1; ++k)
- pnode->tensor()(i,j,k) = pnode->children[n++]->tensor()[0];
- else
- for (size_type i = 0; i < nbl; ++i)
- for (size_type j = 0; j < nbc3; ++j)
- for (size_type k = 0; k < nbc2; ++k)
- for (size_type l = 0; l < nbc1; ++l)
- pnode->tensor()(i,j,k,l)
- = pnode->children[n++]->tensor()[0];
- }
- }
- if (all_cte) tree.clear_children(pnode);
+ int to_add = int(pnode->nb_test_functions() + pnode->nbc1)
+ - int(pnode->tensor().sizes().size());
+ GMM_ASSERT1(to_add >= 0 && to_add <=2, "Internal error");
+ if (to_add) {
+ mi = pnode->tensor().sizes();
+ mi.resize(pnode->nbc1+pnode->nb_test_functions());
+ for (int i = int(mi.size()-1); i >= to_add; --i)
+ mi[i] = mi[i-to_add];
+ for (int i = 0; i < to_add; ++i) mi[i] = 2;
+
+ }
+
+ if (pnode->nb_test_functions()) {
+
+
+
+ pnode->tensor().adjust_sizes(mi);
+ // test pour savoir si c'est bon et sinon, ajout de 1 ou deux ...
+ // comment savoir, en fait ??
+
+ }
+ if (all_cte) {
+ bool all_zero = true;
+ for (size_type i = 0; i < pnode->children.size(); ++i) {
+ pnode->tensor()[i] = pnode->children[i]->tensor()[0];
+ if (pnode->tensor()[i] != scalar_type(0)) all_zero = false;
+ }
+ if (all_zero)
+ pnode->node_type = GA_NODE_ZERO;
+ else
+ pnode->node_type = GA_NODE_CONSTANT;
+ tree.clear_children(pnode);
+ }
}
break;
@@ -2864,8 +2861,6 @@ namespace getfem {
for (size_type i = 0; i < N; ++i) {
factor.clear_children(new_pnode);
new_pnode->node_type = GA_NODE_C_MATRIX;
- new_pnode->nbc1 = meshdim;
- new_pnode->nbc2 = new_pnode->nbc3 = 1;
new_pnode->t.adjust_sizes(mi);
new_pnode->children.resize(N*meshdim);
for (size_type j = 0; j < N; ++j) {
@@ -2873,7 +2868,7 @@ namespace getfem {
if (j == i) {
pga_tree_node param_node = new_pnode->children[k*N+j]
= new ga_tree_node(GA_NODE_PARAMS, pnode->pos, pnode->expr);
- new_pnode->children[k*N+j]->parent = new_pnode;
+ new_pnode->children[k+j*meshdim]->parent = new_pnode;
param_node->children.resize(2);
param_node->children[0]
= new ga_tree_node(GA_NODE_NORMAL, pnode->pos, pnode->expr);
@@ -2884,10 +2879,11 @@ namespace getfem {
param_node->children[1]->init_scalar_tensor(scalar_type(k));
} else {
- new_pnode->children[k*N+j]
+ new_pnode->children[k+j*meshdim]
= new ga_tree_node(GA_NODE_ZERO, pnode->pos, pnode->expr);
- new_pnode->children[k*N+j]->init_scalar_tensor(scalar_type(0));
- new_pnode->children[k*N+j]->parent = new_pnode;
+ new_pnode->children[k+j*meshdim]
+ ->init_scalar_tensor(scalar_type(0));
+ new_pnode->children[k+j*meshdim]->parent = new_pnode;
}
}
}
@@ -3313,6 +3309,10 @@ namespace getfem {
tree.clear_children(pnode->children[i]);
}
}
+ // cout << "After : "; ga_print_node(pnode, cout); cout << endl;
+ // cout << "After : "; ga_print_node(pnode->parent, cout); cout << endl;
+ // ga_node_analysis(tree, workspace, pnode, m, 1, true, false, 1);
+ // cout << "After : "; ga_print_node(pnode->parent, cout); cout << endl;
break;
case GA_NODE_PARAMS:
@@ -4209,25 +4209,16 @@ namespace getfem {
}
}
if (m.dim() > 1) {
- size_type nbl = pnode->children.size() /
- (pnode->nbc1*pnode->nbc2*pnode->nbc3);
- if (pnode->nbc1==1 && pnode->nbc2==1 && pnode->nbc3==1)
- pnode->nbc1 = m.dim();
- else if (pnode->nbc2==1 && pnode->nbc3==1)
- { pnode->nbc2 = pnode->nbc1; pnode->nbc1 = m.dim(); }
- else if (pnode->nbc3==1)
- { pnode->nbc3 = pnode->nbc2; pnode->nbc2 = pnode->nbc1;
- pnode->nbc1 = m.dim(); }
- else GMM_ASSERT1(false, "Sorry this exceed the current limit of "
- "constant tensors (limited to order four)");
- pnode->children.resize(pnode->nbc1*pnode->nbc2*pnode->nbc3*nbl,
- nullptr);
- for (size_type i = pnode->children.size()-1; i > 0; --i) {
- if (i % m.dim())
- tree.copy_node(pnode->children[i/m.dim()], pnode,
- pnode->children[i]);
- else
- std::swap(pnode->children[i/m.dim()], pnode->children[i]);
+ size_type orgsize = pnode->children.size();
+ mi = pnode->tensor().sizes();
+ mi.push_back(m.dim());
+ pnode->nbc1 += 1;
+ pnode->tensor().adjust_sizes(mi);
+
+ pnode->children.resize(orgsize*m.dim(), nullptr);
+ for (size_type i = orgsize; i < pnode->children.size(); ++i) {
+ tree.copy_node(pnode->children[i % orgsize], pnode,
+ pnode->children[i]);
}
for (size_type i = 0; i < pnode->children.size(); ++i) {
pga_tree_node child = pnode->children[i];
@@ -4235,7 +4226,7 @@ namespace getfem {
tree.insert_node(child, GA_NODE_PARAMS);
tree.add_child(child->parent, GA_NODE_CONSTANT);
child->parent->children[1]
- ->init_scalar_tensor(scalar_type(1+i%m.dim()));
+ ->init_scalar_tensor(scalar_type(1+i/orgsize));
}
}
}
@@ -4290,8 +4281,6 @@ namespace getfem {
if (mark1) {
ga_node_grad(tree, workspace, m, pg2->children[ch2]);
}
- ga_print_node(pg1, cout); cout << endl;
- ga_print_node(pg2, cout); cout << endl;
}
#ifdef continue_here
diff --git a/src/getfem_generic_assembly_tree.cc
b/src/getfem_generic_assembly_tree.cc
index 7e0320d..71a7eac 100644
--- a/src/getfem_generic_assembly_tree.cc
+++ b/src/getfem_generic_assembly_tree.cc
@@ -354,28 +354,6 @@ namespace getfem {
current_node->nbc1 = current_node->nbc2 = current_node->nbc3 = 0;
}
- void ga_tree::zip_matrix(const pga_tree_node source_node) {
- GMM_ASSERT1(current_node->node_type == GA_NODE_C_MATRIX &&
- source_node->node_type == GA_NODE_C_MATRIX,
- "Internal error");
- size_type target_size = current_node->children.size();
- size_type source_size = source_node->children.size();
- size_type last_dim_size = target_size/source_size;
- GMM_ASSERT1(target_size == source_size*last_dim_size,
- "Internal error, " << target_size << " != " <<
- source_size << "*" << last_dim_size);
- std::vector<pga_tree_node> new_children;
- for (size_type i = 0; i < source_size; ++i) {
- for (size_type j = 0; j < last_dim_size; ++j)
- new_children.push_back(current_node->children[i*last_dim_size+j]);
- new_children.push_back(source_node->children[i]);
- source_node->children[i]->parent = current_node;
- }
- source_node->children.resize(0); // so that the destructor of source_node
- // will not destruct the children
- current_node->children = new_children;
- }
-
void ga_tree::add_op(GA_TOKEN_TYPE op_type, size_type pos,
pstring expr) {
while (current_node && current_node->parent &&
@@ -546,9 +524,10 @@ namespace getfem {
return false;
break;
case GA_NODE_C_MATRIX:
- if (pnode1->nbc1 != pnode2->nbc1 || pnode1->nbc2 != pnode2->nbc2 ||
- pnode1->nbc3 != pnode2->nbc3)
- return false;
+ if (pnode1->t.sizes().size() != pnode2->t.sizes().size()) return false;
+ for (size_type i = 0; i < pnode1->t.sizes().size(); ++i)
+ if (pnode1->t.sizes()[i] != pnode2->t.sizes()[i]) return false;
+ if (pnode1->nbc1 != pnode2->nbc1) return false;
break;
case GA_NODE_INTERPOLATE_FILTER:
if (pnode1->interpolate_name.compare(pnode2->interpolate_name) ||
@@ -666,7 +645,7 @@ namespace getfem {
case 1:
str << "[";
for (size_type i = 0; i < pnode->tensor_proper_size(0); ++i) {
- if (i != 0) str << "; ";
+ if (i != 0) str << ";";
str << (nt ? scalar_type(0) : pnode->tensor()[i]);
}
str << "]";
@@ -705,10 +684,10 @@ namespace getfem {
}
break;
- case 5: case 6:
+ case 5: case 6: case 7: case 8:
str << "Reshape([";
for (size_type i = 0; i < pnode->tensor_proper_size(); ++i) {
- if (i != 0) str << "; ";
+ if (i != 0) str << ";";
str << (nt ? scalar_type(0) : pnode->tensor()[i]);
}
str << "]";
@@ -1060,33 +1039,70 @@ namespace getfem {
break;
case GA_NODE_C_MATRIX:
- {
- GMM_ASSERT1(pnode->children.size(), "Invalid tree");
- size_type nbc1 = pnode->nbc1, nbc2 = pnode->nbc2;
- size_type nbc3 = pnode->nbc3;
- size_type nbcl = pnode->children.size()/(nbc1*nbc2*nbc3);
- if (nbc1 > 1) str << "[";
- for (size_type i1 = 0; i1 < nbc1; ++i1) {
- if (i1 != 0) str << ",";
- if (nbc2 > 1) str << "[";
- for (size_type i2 = 0; i2 < nbc2; ++i2) {
- if (i2 != 0) str << ",";
- if (nbc3 > 1) str << "[";
- for (size_type i3 = 0; i3 < nbc3; ++i3) {
- if (i3 != 0) str << ",";
- if (nbcl > 1) str << "[";
- for (size_type i4 = 0; i4 < nbcl; ++i4) {
- if (i4 != 0) str << ",";
- size_type ii = ((i4*nbc3 + i3)*nbc2 + i2)*nbc1 + i1;
- ga_print_node(pnode->children[ii], str);
- }
- if (nbcl > 1) str << "]";
- }
- if (nbc3 > 1) str << "]";
- }
- if (nbc2 > 1) str << "]";
- }
- if (nbc1 > 1) str << "]";
+ GMM_ASSERT1(pnode->children.size(), "Invalid tree");
+ GMM_ASSERT1(pnode->nbc1 == pnode->tensor_order(), "Invalid C_MATRIX");
+ switch (pnode->tensor_order()) {
+ case 0:
+ ga_print_node(pnode->children[0], str);
+ break;
+
+ case 1:
+ str << "[";
+ for (size_type i = 0; i < pnode->tensor_proper_size(0); ++i) {
+ if (i != 0) str << ";";
+ ga_print_node(pnode->children[i], str);
+ }
+ str << "]";
+ break;
+
+ case 2: case 3: case 4:
+ {
+ size_type ii(0);
+ size_type n0 = pnode->tensor_proper_size(0);
+ size_type n1 = pnode->tensor_proper_size(1);
+ size_type n2 = ((pnode->tensor_order() > 2) ?
+ pnode->tensor_proper_size(2) : 1);
+ size_type n3 = ((pnode->tensor_order() > 3) ?
+ pnode->tensor_proper_size(3) : 1);
+ if (n3 > 1) str << "[";
+ for (size_type l = 0; l < n3; ++l) {
+ if (l != 0) str << ",";
+ if (n2 > 1) str << "[";
+ for (size_type k = 0; k < n2; ++k) {
+ if (k != 0) str << ",";
+ if (n1 > 1) str << "[";
+ for (size_type j = 0; j < n1; ++j) {
+ if (j != 0) str << ",";
+ if (n0 > 1) str << "[";
+ for (size_type i = 0; i < n0; ++i) {
+ if (i != 0) str << ",";
+ ga_print_node(pnode->children[ii++], str);
+ }
+ if (n0 > 1) str << "]";
+ }
+ if (n1 > 1) str << "]";
+ }
+ if (n2 > 1) str << "]";
+ }
+ if (n3 > 1) str << "]";
+ }
+ break;
+
+ case 5: case 6: case 7: case 8:
+ str << "Reshape([";
+ for (size_type i = 0; i < pnode->tensor_proper_size(); ++i) {
+ if (i != 0) str << ";";
+ ga_print_node(pnode->children[i], str);
+ }
+ str << "]";
+ for (size_type i = 0; i < pnode->tensor_order(); ++i) {
+ if (i != 0) str << ", ";
+ str << pnode->tensor_proper_size(i);
+ }
+ str << ")";
+ break;
+
+ default: GMM_ASSERT1(false, "Invalid tensor dimension");
}
break;
@@ -1664,75 +1680,68 @@ namespace getfem {
GA_TOKEN_TYPE r_type;
size_type nbc1(0), nbc2(0), nbc3(0), n1(0), n2(0), n3(0);
size_type tensor_order(1);
- bool foundcomma(false), foundsemi(false), nested_format(false);
-
- tree.add_matrix(token_pos, expr);
- do {
- r_type = ga_read_term(expr, pos, sub_tree, macro_dict);
+ bool foundcomma(false), foundsemi(false);
- if (sub_tree.root->node_type == GA_NODE_C_MATRIX) {
- // nested format
- if (r_type != GA_COMMA && r_type != GA_RBRACKET)
- // in the nested format only "," and "]" are expected
- ga_throw_error(expr, pos-1, "Bad explicit "
- "vector/matrix/tensor format.")
- else if (sub_tree.root->nbc3 != 1)
- // the sub-tensor to be merged cannot be of fourth order
- ga_throw_error(expr, pos-1, "Definition of explicit "
- "tensors is limitted to the fourth order. "
- "Limit exceeded.")
- else if (foundsemi || // Cannot mix with the non-nested format.
- (sub_tree.root->children.size() > 1 &&
- // The sub-tensor cannot be a column vector [a;b],
- sub_tree.root->nbc1 == 1))
- // the nested format only accepts row vectors [a,b]
- // and converts them to column vectors internally
+ r_type = ga_read_term(expr, pos, sub_tree, macro_dict);
+ size_type nb_comp = 0;
+ tree.add_matrix(token_pos, expr);
+
+ if (sub_tree.root->node_type == GA_NODE_C_MATRIX) { // nested format
+ bgeot::multi_index mii;
+ do {
+ if (nb_comp) {
+ sub_tree.clear();
+ r_type = ga_read_term(expr, pos, sub_tree, macro_dict);
+ }
+ // in the nested format only "," and "]" are expected
+ if (sub_tree.root->node_type != GA_NODE_C_MATRIX ||
+ (r_type != GA_COMMA && r_type != GA_RBRACKET))
ga_throw_error(expr, pos-1, "Bad explicit "
- "vector/matrix/tensor format.") // (see
below)
-
- // convert a row vector [a,b] to a column vector [a;b]
- if (sub_tree.root->children.size() == sub_tree.root->nbc1)
- sub_tree.root->nbc1 = 1;
-
- nested_format = true;
-
- size_type sub_tensor_order = 3;
- if (sub_tree.root->nbc1 == 1)
- sub_tensor_order = 1;
- else if (sub_tree.root->nbc2 == 1)
- sub_tensor_order = 2;
-
- if (tensor_order == 1) {
- if (nbc1 != 0 || nbc2 != 0 || nbc3 != 0)
- ga_throw_error(expr, pos-1, "Bad explicit "
- "vector/matrix/tensor format.");
- nbc1 = 1;
- nbc2 = sub_tree.root->nbc1;
- nbc3 = sub_tree.root->nbc2;
- tensor_order = sub_tensor_order + 1;
- } else {
- if ((tensor_order != sub_tensor_order + 1) ||
- (tensor_order > 2 && nbc2 != sub_tree.root->nbc1) ||
- (tensor_order > 3 && nbc3 != sub_tree.root->nbc2))
- ga_throw_error(expr, pos-1, "Bad explicit "
+ "vector/matrix/tensor format.");
+
+ // convert a row vector [a,b] to a column vector [a;b]
+ if (sub_tree.root->marked &&
+ sub_tree.root->tensor().sizes()[0] == 1 &&
+ sub_tree.root->tensor().size() != 1) {
+ bgeot::multi_index mi = sub_tree.root->tensor().sizes();
+ for (size_type i = mi.size()-1; i > 0; i--)
+ mi[i-1] = mi[i];
+ mi.pop_back();
+ sub_tree.root->tensor().adjust_sizes(mi);
+ }
+ if (!nb_comp) mii = sub_tree.root->tensor().sizes();
+ else {
+ const bgeot::multi_index &mi=sub_tree.root->tensor().sizes();
+ bool cmp = true;
+ if (mii.size() == mi.size()) {
+ for (size_type i = 0; i < mi.size(); ++i)
+ if (mi[i] != mii[i]) cmp = false;
+ } else cmp = false;
+ if (!cmp)
+ ga_throw_error(expr, pos-1, "Bad explicit "
"vector/matrix/tensor format.");
- nbc1 += 1;
- }
-
- tree.zip_matrix(sub_tree.root);
- sub_tree.clear();
- if (r_type == GA_RBRACKET) {
- tree.current_node->nbc1 = nbc1;
- tree.current_node->nbc2 = nbc2;
- tree.current_node->nbc3 = nbc3;
- }
-
- } else {
-
- // non-nested format
-
- tree.add_sub_tree(sub_tree);
-
+ }
+ for (size_type i = 0; i < sub_tree.root->children.size(); ++i) {
+ sub_tree.root->children[i]->parent = tree.current_node;
+ tree.current_node->children.push_back
+ (sub_tree.root->children[i]);
+ }
+ sub_tree.root->children.resize(0);
+ nb_comp++;
+ } while (r_type != GA_RBRACKET);
+ tree.current_node->marked = false;
+ mii.push_back(nb_comp);
+ tree.current_node->tensor().adjust_sizes(mii);
+ } else { // non nested format
+ do {
+ if (nb_comp) {
+ sub_tree.clear();
+ r_type = ga_read_term(expr, pos, sub_tree, macro_dict);
+ }
+ nb_comp++;
+
+ tree.add_sub_tree(sub_tree);
+
++n1; ++n2; ++n3;
if (tensor_order < 2) ++nbc1;
if (tensor_order < 3) ++nbc2;
@@ -1781,12 +1790,36 @@ namespace getfem {
"separated by ',', ';', ',,' and ';;' and "
"be ended by ']'.");
}
- }
-
- } while (r_type != GA_RBRACKET);
-
- state = 2;
+
+ } while (r_type != GA_RBRACKET);
+ bgeot::multi_index mi;
+ nbc1 = tree.current_node->nbc1; nbc2 = tree.current_node->nbc2;
+ nbc3 = tree.current_node->nbc3;
+
+ size_type nbl = tree.current_node->children.size()
+ / (nbc2 * nbc1 * nbc3);
+ switch(tensor_order) {
+ case 1:
+ mi.push_back(1); mi.push_back(nbc1); break;
+ case 2:
+ mi.push_back(nbl); if (nbc1 > 1) mi.push_back(nbc1); break;
+ case 3: mi.push_back(nbl); mi.push_back(nbc2);
mi.push_back(nbc1); break;
+ case 4: mi.push_back(nbl); mi.push_back(nbc3);
mi.push_back(nbc2); mi.push_back(nbc1); break;
+ default: GMM_ASSERT1(false, "Internal error");
+ }
+ tree.current_node->tensor().adjust_sizes(mi);
+ std::vector<pga_tree_node> children = tree.current_node->children;
+ auto it = tree.current_node->children.begin();
+ for (size_type i = 0; i < nbc1; ++i)
+ for (size_type j = 0; j < nbc2; ++j)
+ for (size_type k = 0; k < nbc3; ++k)
+ for (size_type l = 0; l < nbl; ++l, ++it)
+ *it = children[i+nbc1*(j+nbc2*(k+nbc3*l))];
+ tree.current_node->marked = true;
+ }
}
+ tree.current_node->nbc1 = tree.current_node->tensor().sizes().size();
+ state = 2;
break;
default:
- [Getfem-commits] [getfem-commits] devel-yves-generic-assembly-modifs updated (7df7e67 -> d561544), Yves Renard, 2018/05/10
- [Getfem-commits] (no subject), Yves Renard, 2018/05/10
- [Getfem-commits] (no subject), Yves Renard, 2018/05/10
- [Getfem-commits] (no subject), Yves Renard, 2018/05/10
- [Getfem-commits] (no subject), Yves Renard, 2018/05/10
- [Getfem-commits] (no subject), Yves Renard, 2018/05/10
- [Getfem-commits] (no subject), Yves Renard, 2018/05/10
- [Getfem-commits] (no subject),
Yves Renard <=