[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] r5233 - in /trunk/getfem: doc/sphinx/source/userdoc/gas
From: |
Yves . Renard |
Subject: |
[Getfem-commits] r5233 - in /trunk/getfem: doc/sphinx/source/userdoc/gasm_high.rst src/getfem_generic_assembly.cc |
Date: |
Mon, 01 Feb 2016 13:13:52 -0000 |
Author: renard
Date: Mon Feb 1 14:13:50 2016
New Revision: 5233
URL: http://svn.gna.org/viewcvs/getfem?rev=5233&view=rev
Log:
Adding Skew(m) operator
Modified:
trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst
trunk/getfem/src/getfem_generic_assembly.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=5233&r1=5232&r2=5233&view=diff
==============================================================================
--- trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst (original)
+++ trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst Mon Feb 1
14:13:50 2016
@@ -541,18 +541,20 @@
The command ``Reshape(t, i, j, ...)`` reshapes the tensor ``t`` (which could
be an expression). The only constraint is that the number of components should
be compatible. For instance ``Reshape(Grad_u, 1, meshdim)`` is equivalent to
``Grad_u'`` for u a scalar variable. Note that the order of the components
remain unchanged and are classically stored in Fortran order for compatibility
with Blas/Lapack.
-Trace, Deviator and Sym operators
----------------------------------
-
-Trace, Deviator and Sym operators are linear operators acting on square
matrices:
+Trace, Deviator, Sym and Skew operators
+---------------------------------------
+
+Trace, Deviator, Sym and Skew operators are linear operators acting on square
matrices:
- ``Trace(m)`` gives the trace (sum of diagonal components) of a square
matrix ``m``.
- - ``Deviator(m)`` gives the deviator of a square matrix ``m``. It is
equivalent to ``m - Trace(m)*Id(meshdim)/meshdim``.
+ - ``Deviator(m)`` gives the deviator of a square matrix ``m``. It is
equivalent to ``m - Trace(m)*Id(dim)/dim``.
- ``Sym(m)`` gives the symmetric part of a square matrix ``m``, i.e. ``(m +
m')/2``.
-The three operators can be applied on test functions. Which means that for
instance both ``Trace(Grad_u)`` and ``Trace(Grad_Test_u)`` is valid when
``Grad_u`` is a square matrix (i.e. ``u`` a vector field of the same dimension
as the mesh).
+ - ``Skew(m)`` gives the skew-symmetric part of a square matrix ``m``, i.e.
``(m - m')/2``.
+
+The four operators can be applied on test functions. Which means that for
instance both ``Trace(Grad_u)`` and ``Trace(Grad_Test_u)`` is valid when
``Grad_u`` is a square matrix (i.e. ``u`` a vector field of the same dimension
as the mesh).
Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=5233&r1=5232&r2=5233&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Mon Feb 1 14:13:50 2016
@@ -90,6 +90,7 @@
GA_COLON, // ':'
GA_QUOTE, // ''' transpose
GA_SYM, // 'Sym' operator
+ GA_SKEW, // 'Skew' operator
GA_TRACE, // 'Trace' operator
GA_DEVIATOR, // 'Deviator' operator
GA_INTERPOLATE, // 'Interpolate' operation
@@ -145,6 +146,7 @@
ga_operator_priorities[GA_QUOTE] = 3;
ga_operator_priorities[GA_UNARY_MINUS] = 3;
ga_operator_priorities[GA_SYM] = 4;
+ ga_operator_priorities[GA_SKEW] = 4;
ga_operator_priorities[GA_TRACE] = 4;
ga_operator_priorities[GA_DEVIATOR] = 4;
ga_operator_priorities[GA_PRINT] = 4;
@@ -215,6 +217,8 @@
}
if (expr.compare(token_pos, token_length, "Sym") == 0)
return GA_SYM;
+ if (expr.compare(token_pos, token_length, "Skew") == 0)
+ return GA_SKEW;
if (expr.compare(token_pos, token_length, "Trace") == 0)
return GA_TRACE;
if (expr.compare(token_pos, token_length, "Deviator") == 0)
@@ -647,6 +651,7 @@
pga_tree_node new_node = new ga_tree_node(op_type, pos);
if (current_node) {
if (op_type == GA_UNARY_MINUS || op_type == GA_SYM
+ || op_type == GA_SKEW
|| op_type == GA_TRACE || op_type == GA_DEVIATOR
|| op_type == GA_PRINT) {
current_node->children.push_back(new_node);
@@ -1159,7 +1164,10 @@
} else if (pnode->op_type == GA_QUOTE) {
GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
ga_print_node(pnode->children[0], str); str << "'";
- } else if (pnode->op_type == GA_SYM) {
+ } else if (pnode->op_type == GA_SKEW) {
+ GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
+ str << "Skew("; ga_print_node(pnode->children[0], str); str << ")";
+ }else if (pnode->op_type == GA_SYM) {
GMM_ASSERT1(pnode->children.size() == 1, "Invalid tree");
str << "Sym("; ga_print_node(pnode->children[0], str); str << ")";
} else if (pnode->op_type == GA_TRACE) {
@@ -1462,6 +1470,10 @@
tree.add_op(GA_SYM, token_pos);
state = 1; break;
+ case GA_SKEW:
+ tree.add_op(GA_SKEW, token_pos);
+ state = 1; break;
+
case GA_TRACE:
tree.add_op(GA_TRACE, token_pos);
state = 1; break;
@@ -4135,7 +4147,7 @@
base_tensor &t;
const base_tensor &tc1;
virtual int exec(void) {
- GA_DEBUG_INFO("Instruction: transpose");
+ GA_DEBUG_INFO("Instruction: symmetric part");
GA_DEBUG_ASSERT(t.size() == tc1.size(), "Wrong sizes");
size_type order = t.sizes().size();
size_type s1 = t.sizes()[order-2], s2 = t.sizes()[order-1];
@@ -4150,6 +4162,28 @@
return 0;
}
ga_instruction_sym(base_tensor &t_, const base_tensor &tc1_)
+ : t(t_), tc1(tc1_) {}
+ };
+
+ struct ga_instruction_skew : public ga_instruction {
+ base_tensor &t;
+ const base_tensor &tc1;
+ virtual int exec(void) {
+ GA_DEBUG_INFO("Instruction: skew-symmetric part");
+ GA_DEBUG_ASSERT(t.size() == tc1.size(), "Wrong sizes");
+ size_type order = t.sizes().size();
+ size_type s1 = t.sizes()[order-2], s2 = t.sizes()[order-1];
+ size_type s = t.size() / (s1*s2);
+ for (size_type i = 0; i < s1; ++i)
+ for (size_type j = 0; j < s2; ++j) {
+ base_tensor::iterator it = t.begin() + s*(i + s1*j);
+ base_tensor::const_iterator it1 = tc1.begin() + s*(i + s1*j),
+ it1T = tc1.begin() + s*(j + s2*i);
+ for (size_type k = 0; k < s; ++k) *it++ = 0.5*(*it1++ - *it1T++);
+ }
+ return 0;
+ }
+ ga_instruction_skew(base_tensor &t_, const base_tensor &tc1_)
: t(t_), tc1(tc1_) {}
};
@@ -6674,13 +6708,21 @@
}
break;
- case GA_SYM:
+ case GA_SYM: case GA_SKEW:
if (dim0 != 2 || size0.back() != size0[size0.size()-2])
- ga_throw_error(expr, pnode->pos, "Sym operator is for "
+ ga_throw_error(expr, pnode->pos, "Sym and Skew operators are for "
"square matrices only.");
mi = size0;
- if (child0->tensor_proper_size() == 1)
- { tree.replace_node_by_child(pnode, 0); pnode = child0; break; }
+ if (child0->tensor_proper_size() == 1) {
+ if (pnode->op_type == GA_SYM)
+ { tree.replace_node_by_child(pnode, 0); pnode = child0; break; }
+ else {
+ pnode->node_type = GA_NODE_ZERO;
+ gmm::clear(pnode->t.as_vector());
+ tree.clear_children(pnode);
+ break;
+ }
+ }
pnode->t.adjust_sizes(mi);
pnode->test_function_type = child0->test_function_type;
@@ -6695,7 +6737,10 @@
pnode->test_function_type = 0;
for (size_type i = 0; i < mi.back(); ++i)
for (size_type j = 0; j < mi.back(); ++j)
- pnode->t(j, i) = 0.5*(child0->t(j,i) + child0->t(i,j));
+ if (pnode->op_type == GA_SYM)
+ pnode->t(j, i) = 0.5*(child0->t(j,i) + child0->t(i,j));
+ else
+ pnode->t(j, i) = 0.5*(child0->t(j,i) - child0->t(i,j));
tree.clear_children(pnode);
} else if (child0->node_type == GA_NODE_ZERO) {
pnode->node_type = GA_NODE_ZERO;
@@ -6708,6 +6753,9 @@
{
mi = size0;
size_type N = (child0->tensor_proper_size() == 1) ? 1 : mi.back();
+
+ if (child0->tensor_proper_size() == 1)
+ { tree.replace_node_by_child(pnode, 0); pnode = child0; break; }
if ((dim0 != 2 && child0->tensor_proper_size() != 1) ||
(dim0 == 2 && mi[mi.size()-2] != N))
@@ -6751,6 +6799,13 @@
(dim0 == 2 && mi[mi.size()-2] != N))
ga_throw_error(expr, pnode->pos,
"Deviator operator is for square matrices only.");
+
+ if (child0->tensor_proper_size() == 1) {
+ pnode->node_type = GA_NODE_ZERO;
+ gmm::clear(pnode->t.as_vector());
+ tree.clear_children(pnode);
+ break;
+ }
pnode->t.adjust_sizes(mi);
pnode->test_function_type = child0->test_function_type;
@@ -7834,7 +7889,7 @@
if (child1 == pnode_child) minus_sign = !(minus_sign);
// A remaining minus sign is added at the end if necessary.
break;
- case GA_UNARY_MINUS: case GA_QUOTE: case GA_SYM:
+ case GA_UNARY_MINUS: case GA_QUOTE: case GA_SYM: case GA_SKEW:
case GA_TRACE: case GA_DEVIATOR: case GA_PRINT:
// Copy of the term
result_tree.insert_node(result_tree.root, pnode->node_type);
@@ -8008,7 +8063,7 @@
{ is_constant = false; break; }
break;
- case GA_UNARY_MINUS: case GA_QUOTE: case GA_SYM:
+ case GA_UNARY_MINUS: case GA_QUOTE: case GA_SYM: case GA_SKEW:
case GA_TRACE: case GA_DEVIATOR: case GA_PRINT:
is_constant = child_0_is_constant;
break;
@@ -8664,7 +8719,7 @@
}
break;
- case GA_UNARY_MINUS: case GA_QUOTE: case GA_SYM:
+ case GA_UNARY_MINUS: case GA_QUOTE: case GA_SYM: case GA_SKEW:
case GA_TRACE: case GA_DEVIATOR: case GA_PRINT:
ga_node_derivation(tree, workspace, m, child0, varname,
interpolatename, order);
@@ -9051,7 +9106,7 @@
if (mark0) return ga_node_is_affine(child0);
return ga_node_is_affine(child1);
- case GA_UNARY_MINUS: case GA_QUOTE: case GA_SYM:
+ case GA_UNARY_MINUS: case GA_QUOTE: case GA_SYM: case GA_SKEW:
case GA_TRACE: case GA_DEVIATOR:
case GA_PRINT: case GA_NODE_INTERPOLATE_FILTER:
return ga_node_is_affine(child0);
@@ -10164,6 +10219,14 @@
}
break;
+ case GA_SKEW:
+ {
+ pgai = std::make_shared<ga_instruction_skew>
+ (pnode->t, child0->t);
+ rmi.instructions.push_back(std::move(pgai));
+ }
+ break;
+
case GA_TRACE:
{
size_type N = (child0->tensor_proper_size() == 1) ? 1:size0.back();
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] r5233 - in /trunk/getfem: doc/sphinx/source/userdoc/gasm_high.rst src/getfem_generic_assembly.cc,
Yves . Renard <=