[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: |
Mon, 2 Dec 2019 13:42:49 -0500 (EST) |
branch: fixmisspell
commit c5cce1428dd323f7aa108247a698bbfd2292e32f
Author: Yves Renard <address@hidden>
Date: Mon Dec 2 19:38:23 2019 +0100
adding cross product of three-dimensional vectors
---
doc/sphinx/source/userdoc/gasm_high.rst | 8 +-
src/getfem/getfem_generic_assembly_tree.h | 14 ++--
src/getfem_generic_assembly_compile_and_exec.cc | 73 +++++++++++++++++-
src/getfem_generic_assembly_semantic.cc | 99 ++++++++++++++++++++++++-
src/getfem_generic_assembly_tree.cc | 5 ++
5 files changed, 188 insertions(+), 11 deletions(-)
diff --git a/doc/sphinx/source/userdoc/gasm_high.rst
b/doc/sphinx/source/userdoc/gasm_high.rst
index be8f598..3305c9f 100644
--- a/doc/sphinx/source/userdoc/gasm_high.rst
+++ b/doc/sphinx/source/userdoc/gasm_high.rst
@@ -38,8 +38,12 @@ A specific weak form language has been developed to describe
the weak formulatio
- A certain number of predefined scalar functions (``sin(t)``, ``cos(t)``,
``pow(t,u)``, ``sqrt(t)``, ``sqr(t)``, ``Heaviside(t)``, ...). A scalar
function can be applied to scalar or vector/matrix/tensor expressions. It
applies componentwise. For functions having two arguments (``pow(t,u)``,
``min(t,u)`` ...) if two non-scalar arguments are passed, the dimension have to
be the same. For instance "max([1;2],[0;3])" will return "[1;3]".
- - A certain number of operators: ``+``, ``-``, ``*``, ``/``, ``:``, ``.``,
``.*``, ``./``, ``@``, ``'``.
+ - A certain number of operations: ``+``, ``-``, ``*``, ``/``, ``:``, ``.``,
``.*``, ``./``, ``@``, ``'``, ``Cross_product(v1,v2)``.
+ - A certain number of linear operator: ``Trace(M)``, ``Sym(M)``,
``Skew(M)``, ...
+
+ - A certain number of nonlinear operator: ``Norm(V)``, ``Det(M)``,
``Sym(M)``, ``Skew(M)``, ...
+
- Some constants: ``pi``, ``meshdim`` (the dimension of the current mesh),
``qdim(u)`` and ``qdims(u)`` the dimensions of the variable ``u`` (the size for
fixed size variables and the dimension of the vector field for FEM variables),
``Id(n)`` the identity :math:`n\times n` matrix.
- Parentheses can be used to change the operations order in a standard way.
For instance ``(1+2)*4`` or ``(u+v)*Test_u`` are valid expressions.
@@ -477,6 +481,8 @@ A certain number of binary operations between tensors are
available:
- ``@`` stands for the tensor product.
+ - ``Cross_product(V, W)`` stands for the cross product (vector product) of
``V`` and ``W``. Defined only for three-dimensional vectors.
+
- ``Contract(A, i, B, j)`` stands for the contraction of tensors A and B
with respect to the ith index of A and jth index of B. The first index is
numbered 1. For instance ``Contract(V,1,W,1)`` is equivalent to ``V.W`` for two
vectors ``V`` and ``W``.
- ``Contract(A, i, j, B, k, l)`` stands for the double contraction of
tensors A and B with respect to indices i,j of A and indices k,l of B. The
first index is numbered 1. For instance ``Contract(A,1,2,B,1,2)`` is equivalent
to ``A:B`` for two matrices ``A`` and ``B``.
diff --git a/src/getfem/getfem_generic_assembly_tree.h
b/src/getfem/getfem_generic_assembly_tree.h
index 7715299..83f903d 100644
--- a/src/getfem/getfem_generic_assembly_tree.h
+++ b/src/getfem/getfem_generic_assembly_tree.h
@@ -84,17 +84,17 @@ namespace getfem {
GA_QUOTE, // ''' transpose
GA_COLON_EQ, // ':=' macro def
GA_DEF, // 'Def' macro def
- GA_SYM, // 'Sym' operator
- GA_SKEW, // 'Skew' operator
- GA_TRACE, // 'Trace' operator
+ GA_SYM, // 'Sym(M)' operator
+ GA_SKEW, // 'Skew(M)' operator
+ GA_TRACE, // 'Trace(M)' operator
GA_DEVIATOR, // 'Deviator' operator
GA_INTERPOLATE, // 'Interpolate' operation
GA_INTERPOLATE_FILTER, // 'Interpolate_filter' operation
GA_INTERPOLATE_DERIVATIVE, // 'Interpolate_derivative' operation
GA_ELEMENTARY, // 'Elementary' operation (operation at the element level)
GA_SECONDARY_DOMAIN, // For the integration on a product of two domains
- GA_XFEM_PLUS, // �valuation on the + side of a level-set for
fem_level_set
- GA_XFEM_MINUS, // �valuation on the - side of a level-set for
fem_level_set
+ GA_XFEM_PLUS, // Evaluation on the + side of a level-set for
fem_level_set
+ GA_XFEM_MINUS, // Evaluation on the - side of a level-set for
fem_level_set
GA_PRINT, // 'Print' Print the tensor
GA_DOT, // '.'
GA_DOTMULT, // '.*' componentwise multiplication
@@ -128,6 +128,7 @@ namespace getfem {
GA_NODE_MACRO_PARAM,
GA_NODE_PARAMS,
GA_NODE_RESHAPE,
+ GA_NODE_CROSS_PRODUCT,
GA_NODE_SWAP_IND,
GA_NODE_IND_MOVE_LAST,
GA_NODE_CONTRACT,
@@ -371,6 +372,9 @@ namespace getfem {
return true;
}
+ inline bool is_constant()
+ { return (node_type == GA_NODE_CONSTANT || node_type == GA_NODE_ZERO); }
+
inline void init_scalar_tensor(scalar_type v)
{ t.init_scalar_tensor(v); test_function_type = 0; }
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index e53f9c8..dcf0921 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -2112,6 +2112,64 @@ namespace getfem {
: t(t_), tc1(tc1_), c(c_) {}
};
+ // Performs Cross product in the presence of test functions
+ struct ga_instruction_cross_product_tf : public ga_instruction {
+ base_tensor &t, &tc1, &tc2;
+ bool inv;
+ virtual int exec() {
+ GA_DEBUG_INFO("Instruction: Cross product with test functions");
+
+ size_type n1 = tc1.size() / 3, n2 = tc2.size() / 3, nn=n1*n2;
+ GA_DEBUG_ASSERT(t.size() == nn*3, "Bad tensor size for cross product");
+ size_type mm=2*nn, n1_2 = 2*n1, n2_2 = 2*n2;
+ base_tensor::iterator it = t.begin(), it2 = tc2.begin();
+
+ if (inv) {
+ for (size_type i = 0; i < n2; ++i, ++it2) {
+ base_tensor::iterator it1 = tc1.begin();
+ for (size_type j = 0; j < n1; ++j, ++it, ++it1) {
+ *it = - it1[n1] *it2[n2_2] + it1[n1_2]*it2[n2];
+ it[nn] = - it1[n1_2]*it2[0] + it1[0] *it2[n2_2];
+ it[mm] = - it1[0] *it2[n2] + it1[n1] *it2[0];
+ }
+ }
+ } else {
+ for (size_type i = 0; i < n2; ++i, ++it2) {
+ base_tensor::iterator it1 = tc1.begin();
+ for (size_type j = 0; j < n1; ++j, ++it, ++it1) {
+ *it = it1[n1] *it2[n2_2] - it1[n1_2]*it2[n2];
+ it[nn] = it1[n1_2]*it2[0] - it1[0] *it2[n2_2];
+ it[mm] = it1[0] *it2[n2] - it1[n1] *it2[0];
+ }
+ }
+ }
+ return 0;
+ }
+ ga_instruction_cross_product_tf(base_tensor &t_, base_tensor &tc1_,
+ base_tensor &tc2_, bool inv_)
+ : t(t_), tc1(tc1_), tc2(tc2_), inv(inv_) {}
+ };
+
+ // Performs Cross product in the absence of test functions
+ struct ga_instruction_cross_product : public ga_instruction {
+ base_tensor &t, &tc1, &tc2;
+ virtual int exec() {
+ GA_DEBUG_INFO("Instruction: Cross product with test functions");
+ GA_DEBUG_ASSERT(t.size() == 3 && tc1.size() == 3 && tc2.size() == 3,
+ "Bad tensor size for cross product");
+ t[0] = tc1[1]*tc2[2] - tc1[2]*tc2[1];
+ t[1] = tc1[2]*tc2[0] - tc1[0]*tc2[2];
+ t[2] = tc1[0]*tc2[1] - tc1[1]*tc2[0];
+ return 0;
+ }
+ ga_instruction_cross_product(base_tensor &t_, base_tensor &tc1_,
+ base_tensor &tc2_)
+ : t(t_), tc1(tc1_), tc2(tc2_) {}
+ };
+
+
+
+
struct ga_instruction_dotmult : public ga_instruction {
base_tensor &t, &tc1, &tc2;
virtual int exec() {
@@ -4944,7 +5002,8 @@ namespace getfem {
case GA_NODE_PREDEF_FUNC: case GA_NODE_OPERATOR: case GA_NODE_SPEC_FUNC:
case GA_NODE_CONSTANT: case GA_NODE_ALLINDICES: case GA_NODE_ZERO:
- case GA_NODE_RESHAPE: case GA_NODE_SWAP_IND: case GA_NODE_IND_MOVE_LAST:
+ case GA_NODE_RESHAPE: case GA_NODE_CROSS_PRODUCT:
+ case GA_NODE_SWAP_IND: case GA_NODE_IND_MOVE_LAST:
case GA_NODE_CONTRACT: case GA_NODE_INTERPOLATE_FILTER:
break;
@@ -6369,6 +6428,18 @@ namespace getfem {
pgai = std::make_shared<ga_instruction_copy_tensor>(pnode->tensor(),
child1->tensor());
rmi.instructions.push_back(std::move(pgai));
+ } else if (child0->node_type == GA_NODE_CROSS_PRODUCT) {
+ pga_tree_node child2 = pnode->children[2];
+ if (child1->test_function_type==2 && child2->test_function_type==1)
+ pgai = std::make_shared<ga_instruction_cross_product_tf>
+ (pnode->tensor(), child2->tensor(), child1->tensor(), true);
+ else if (child1->test_function_type || child2->test_function_type)
+ pgai = std::make_shared<ga_instruction_cross_product_tf>
+ (pnode->tensor(), child1->tensor(), child2->tensor(), false);
+ else
+ pgai = std::make_shared<ga_instruction_cross_product>
+ (pnode->tensor(), child1->tensor(), child2->tensor());
+ rmi.instructions.push_back(std::move(pgai));
} else if (child0->node_type == GA_NODE_IND_MOVE_LAST) {
size_type ind;
ind = size_type(round(pnode->children[2]->tensor()[0])-1);
diff --git a/src/getfem_generic_assembly_semantic.cc
b/src/getfem_generic_assembly_semantic.cc
index 595afe9..5423326 100644
--- a/src/getfem_generic_assembly_semantic.cc
+++ b/src/getfem_generic_assembly_semantic.cc
@@ -398,10 +398,12 @@ namespace getfem {
for (size_type i = 0; i < pnode->children.size(); ++i) {
ga_node_analysis(tree, workspace, pnode->children[i], me,
ref_elt_dim, eval_fixed_size, ignore_X, option);
- all_cte = all_cte && (pnode->children[i]->node_type == GA_NODE_CONSTANT);
+ all_cte = all_cte && (pnode->children[i]->is_constant());
all_sc = all_sc && (pnode->children[i]->tensor_proper_size() == 1);
- GMM_ASSERT1(pnode->children[i]->test_function_type != size_type(-1),
- "internal error on child " << i);
+ if (pnode->children[i]->test_function_type == size_type(-1)) {
+ cerr << "node : "; ga_print_node(pnode, cerr); cerr << endl;
+ GMM_ASSERT1(false, "internal error on child " << i);
+ }
if (pnode->node_type != GA_NODE_PARAMS)
ga_valid_operand(pnode->children[i]);
}
@@ -426,7 +428,8 @@ namespace getfem {
case GA_NODE_PREDEF_FUNC: case GA_NODE_OPERATOR: case GA_NODE_SPEC_FUNC :
case GA_NODE_CONSTANT: case GA_NODE_X: case GA_NODE_ELT_SIZE:
case GA_NODE_ELT_K: case GA_NODE_ELT_B: case GA_NODE_NORMAL:
- case GA_NODE_RESHAPE: case GA_NODE_IND_MOVE_LAST: case GA_NODE_SWAP_IND:
+ case GA_NODE_RESHAPE: case GA_NODE_CROSS_PRODUCT:
+ case GA_NODE_IND_MOVE_LAST: case GA_NODE_SWAP_IND:
case GA_NODE_CONTRACT: case GA_NODE_INTERPOLATE_X:
case GA_NODE_INTERPOLATE_NORMAL: case GA_NODE_SECONDARY_DOMAIN_X:
case GA_NODE_SECONDARY_DOMAIN_NORMAL:
@@ -1626,6 +1629,11 @@ namespace getfem {
pnode->init_scalar_tensor(0);
break;
}
+ if (!(name.compare("Cross_product"))) {
+ pnode->node_type = GA_NODE_CROSS_PRODUCT;
+ pnode->test_function_type = 0;
+ break;
+ }
if (!(name.compare("element_K"))) {
pnode->node_type = GA_NODE_ELT_K;
if (ref_elt_dim == 1)
@@ -2040,6 +2048,37 @@ namespace getfem {
}
}
+ // Cross product of two vectors
+ else if (child0->node_type == GA_NODE_CROSS_PRODUCT) {
+ if (pnode->children.size() != 3)
+ ga_throw_error(pnode->expr, child1->pos,
+ "Wrong number of parameters for Cross_product");
+ pga_tree_node child2 = pnode->children[2];
+
+ if (false && child1->is_constant() && child2->is_constant()) {
+ pnode->node_type = GA_NODE_CONSTANT;
+ pnode->test_function_type = 0;
+ if (child1->tensor_proper_size() != 3 ||
+ child2->tensor_proper_size() != 3)
+ ga_throw_error(pnode->expr, child1->pos, "Cross_product is only "
+ "defined on three-dimensional vectors");
+ pnode->t = child1->t;
+ base_tensor &t0 = pnode->tensor();
+ base_tensor &t1 = child1->tensor(), &t2 = child2->tensor();
+ t0[0] = t1[1]*t2[2] - t1[2]*t2[1];
+ t0[1] = t1[2]*t2[0] - t1[0]*t2[2];
+ t0[2] = t1[0]*t2[1] - t1[1]*t2[0];
+ if (pnode->tensor_is_zero())
+ pnode->node_type = GA_NODE_ZERO;
+ tree.clear_children(pnode);
+ } else {
+ pnode->mult_test(child1, child2);
+ mi = pnode->t.sizes();
+ mi.push_back(3);
+ pnode->t.adjust_sizes(mi);
+ }
+ }
+
// Swap_indices operation
else if (child0->node_type == GA_NODE_SWAP_IND) {
if (pnode->children.size() != 4)
@@ -2720,6 +2759,27 @@ namespace getfem {
if (i != 1)
result_tree.copy_node(pnode->children[i], result_tree.root,
result_tree.root->children[i]);
+ } else if (child0->node_type == GA_NODE_CROSS_PRODUCT) {
+ pga_tree_node child2 = pnode->children[2];
+ result_tree.insert_node(result_tree.root, pnode->node_type);
+ result_tree.root->pos = pnode->pos;
+ result_tree.root->expr = pnode->expr;
+ result_tree.root->children.resize(3, nullptr);
+ if (child1 == pnode_child) {
+ std::swap(result_tree.root->children[1],
+ result_tree.root->children[0]);
+ result_tree.copy_node(pnode->children[0], result_tree.root,
+ result_tree.root->children[0]);
+ result_tree.copy_node(pnode->children[2], result_tree.root,
+ result_tree.root->children[2]);
+ } else if (child2 == pnode_child) {
+ std::swap(result_tree.root->children[2],
+ result_tree.root->children[0]);
+ result_tree.copy_node(pnode->children[0], result_tree.root,
+ result_tree.root->children[0]);
+ result_tree.copy_node(pnode->children[1], result_tree.root,
+ result_tree.root->children[1]);
+ } else GMM_ASSERT1(false, "Corrupted tree");
} else if (child0->node_type == GA_NODE_SWAP_IND) {
GMM_ASSERT1(child1 == pnode_child, "Cannot extract a factor of a "
"Swap_indices size parameter");
@@ -2835,6 +2895,7 @@ namespace getfem {
case GA_NODE_XFEM_MINUS_HESS_TEST: case GA_NODE_XFEM_MINUS_DIVERG_TEST:
case GA_NODE_VAL_TEST: case GA_NODE_GRAD_TEST: case GA_NODE_PREDEF_FUNC:
case GA_NODE_HESS_TEST: case GA_NODE_DIVERG_TEST: case GA_NODE_RESHAPE:
+ case GA_NODE_CROSS_PRODUCT:
case GA_NODE_SWAP_IND: case GA_NODE_IND_MOVE_LAST: case GA_NODE_CONTRACT:
case GA_NODE_ELT_SIZE: case GA_NODE_ELT_K: case GA_NODE_ELT_B:
case GA_NODE_CONSTANT: case GA_NODE_X: case GA_NODE_NORMAL:
@@ -2930,6 +2991,10 @@ namespace getfem {
child0->node_type == GA_NODE_SWAP_IND ||
child0->node_type == GA_NODE_IND_MOVE_LAST) {
is_constant = child_1_is_constant;
+ } else if (child0->node_type == GA_NODE_CROSS_PRODUCT) {
+ bool child_2_is_constant
+ = ga_node_extract_constant_term(tree,pnode->children[2],workspace,m);
+ is_constant = (child_1_is_constant && child_2_is_constant);
} else if (child0->node_type == GA_NODE_CONTRACT) {
for (size_type i = 1; i < pnode->children.size(); ++i) {
if (!(ga_node_extract_constant_term(tree, pnode->children[i],
@@ -3554,6 +3619,22 @@ namespace getfem {
child0->node_type == GA_NODE_IND_MOVE_LAST) {
ga_node_derivation(tree, workspace, m, pnode->children[1],
varname, interpolatename, order);
+ } else if (child0->node_type == GA_NODE_CROSS_PRODUCT) {
+ pga_tree_node child2 = pnode->children[2];
+ bool mark2 = child2->marked;
+ if (mark1 && mark2) {
+ tree.duplicate_with_addition(pnode);
+ ga_node_derivation(tree, workspace, m, child1, varname,
+ interpolatename, order);
+ ga_node_derivation(tree, workspace, m,
+ pnode->parent->children[1]->children[2],
+ varname, interpolatename, order);
+ } else if (mark1) {
+ ga_node_derivation(tree, workspace, m, child1, varname,
+ interpolatename, order);
+ } else
+ ga_node_derivation(tree, workspace, m, child2, varname,
+ interpolatename, order);
} else if (child0->node_type == GA_NODE_CONTRACT) {
if (pnode->children.size() == 4) {
@@ -4523,6 +4604,9 @@ namespace getfem {
ga_node_grad(tree, workspace, m, pnode->children[1]);
tree.add_child(pnode, GA_NODE_CONSTANT);
pnode->children.back()->init_scalar_tensor(scalar_type(m.dim()));
+ } else if (child0->node_type == GA_NODE_CROSS_PRODUCT) {
+ GMM_ASSERT1(false, "Sorry, the gradient of a cross product"
+ "has not been implemented. To be done.");
} else if (child0->node_type == GA_NODE_IND_MOVE_LAST) {
size_type order = pnode->tensor_order();
ga_node_grad(tree, workspace, m, pnode->children[1]);
@@ -4899,6 +4983,13 @@ namespace getfem {
child0->node_type == GA_NODE_SWAP_IND ||
child0->node_type == GA_NODE_IND_MOVE_LAST)
return ga_node_is_affine(child1);
+ if (child0->node_type == GA_NODE_CROSS_PRODUCT) {
+ pga_tree_node child2 = pnode->children[2];
+ bool mark2 = child2->marked;
+ if (mark1 && mark2) return false;
+ if (mark1) return ga_node_is_affine(child1);
+ return ga_node_is_affine(child2);
+ }
if (child0->node_type == GA_NODE_CONTRACT) {
if (pnode->children.size() == 4) {
return ga_node_is_affine(child1);
diff --git a/src/getfem_generic_assembly_tree.cc
b/src/getfem_generic_assembly_tree.cc
index 6978d90..9bf1998 100644
--- a/src/getfem_generic_assembly_tree.cc
+++ b/src/getfem_generic_assembly_tree.cc
@@ -1128,6 +1128,11 @@ namespace getfem {
GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
break;
+ case GA_NODE_CROSS_PRODUCT:
+ str << "Cross_product";
+ GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
+ break;
+
case GA_NODE_SWAP_IND:
str << "Swap_indices";
GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");