[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] r4535 - in /trunk/getfem: ./ doc/sphinx/source/userdoc/
|
From: |
Yves . Renard |
|
Subject: |
[Getfem-commits] r4535 - in /trunk/getfem: ./ doc/sphinx/source/userdoc/ src/ src/getfem/ tests/ |
|
Date: |
Fri, 14 Mar 2014 14:46:34 -0000 |
Author: renard
Date: Fri Mar 14 15:46:33 2014
New Revision: 4535
URL: http://svn.gna.org/viewcvs/getfem?rev=4535&view=rev
Log:
Allowing the use of tensor fields in generic assembly and allowing up to order
6 tensors
Modified:
trunk/getfem/configure.ac
trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst
trunk/getfem/src/getfem/getfem_generic_assembly.h
trunk/getfem/src/getfem_generic_assembly.cc
trunk/getfem/tests/test_assembly.cc
Modified: trunk/getfem/configure.ac
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/configure.ac?rev=4535&r1=4534&r2=4535&view=diff
==============================================================================
--- trunk/getfem/configure.ac (original)
+++ trunk/getfem/configure.ac Fri Mar 14 15:46:33 2014
@@ -159,6 +159,8 @@
AC_CHECK_CXX_FLAG([-fmessage-length=0],CXXFLAGS)
AC_CHECK_CXX_FLAG([-ftemplate-depth-40],CXXFLAGS)
AC_CHECK_CXX_FLAG([-pedantic],CXXFLAGS)
+ AC_CHECK_CXX_FLAG([-std=c++0x],CXXFLAGS)
+ AC_CHECK_CXX_FLAG([-std=c++11],CXXFLAGS)
AC_CHECK_CXX_FLAG([-Wshadow],CXXFLAGS)
AC_CHECK_CXX_FLAG([-Wno-unknown-pragmas],CXXFLAGS)
AC_CHECK_CXX_FLAG([-Wpointer-arith],CXXFLAGS)
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=4535&r1=4534&r2=4535&view=diff
==============================================================================
--- trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst (original)
+++ trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst Fri Mar 14
15:46:33 2014
@@ -44,7 +44,7 @@
- A certain number of operators: ``+``, ``-``, ``*``, ``/``, ``:``, ``.``,
``.*``, ``./``, address@hidden, ``'``.
- - Some constants : ``pi``, ``meshdim`` (the dimension of the current mesh),
``qdim(u)`` the dimension of the variable ``u`` (the size for fixed size
variables and the dimension of the vector field for f.e.m. variables),
``Id(n)`` the identity :math:`n\times n` matrix.
+ - 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 f.e.m.
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 correct.
@@ -347,7 +347,7 @@
The tensors
***********
-Basically, what is manipulated in the generic assembly language are tensors.
This can be order 0 tensors in scalar expressions (for instance in
``3+sin(pi/2)``), order 1 tensors in vector expressions (such as ``x.x`` or
``Grad_u`` if u is a scalar variable), order 2 tensors for matrix expressions
and so on. For efficiency reasons, the language manipulates tensors up to order
four. The language could be easily extended to support tensors of order greater
than four but it may lead to inefficient computations. When an expression
contains tests functions (as in ``Trace(Grad_Test_u)`` for a vector field
``u``), the computation is done for each test functions, which means that the
tensor implicitly have a supplementary component. This means that, implicitly,
the maximal order of manipulated tensors are in fact six (in
``Grad_Test_u:Grad_Test2_u`` there are two components implicitly added for
first and second order test functions).
+Basically, what is manipulated in the generic assembly language are tensors.
This can be order 0 tensors in scalar expressions (for instance in
``3+sin(pi/2)``), order 1 tensors in vector expressions (such as ``x.x`` or
``Grad_u`` if u is a scalar variable), order 2 tensors for matrix expressions
and so on. For efficiency reasons, the language manipulates tensors up to order
six. The language could be easily extended to support tensors of order greater
than six but it may lead to inefficient computations. When an expression
contains tests functions (as in ``Trace(Grad_Test_u)`` for a vector field
``u``), the computation is done for each test functions, which means that the
tensor implicitly have a supplementary component. This means that, implicitly,
the maximal order of manipulated tensors are in fact six (in
``Grad_Test_u:Grad_Test2_u`` there are two components implicitly added for
first and second order test functions).
Order four tensors are necessary for instance to express elasticity tensors or
in general to obtain the tangent term for vector valued unknowns.
@@ -487,7 +487,12 @@
Explicit order four 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,2,,1,1;1,2;;1,1;1,2,,1,1;1,2]`` is a 2x2x2x2
valid tensor. Note that constant fourth order tensors can also be obtained by
the tensor product of two constant matrices. Note also that constant third
order tensors are not supported.
+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,2,,1,1;1,2;;1,1;1,2,,1,1;1,2]`` is a 2x2x2x2
valid tensor. Note that constant fourth order tensors can also be obtained by
the tensor product of two constant matrices.
+
+Explicit order five or six 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.
Access to tensor components
@@ -501,7 +506,8 @@
- ``pi``: the constant Pi.
- ``meshdim``: the dimension of the current mesh (i.e. size of geometrical
nodes)
- ``Id(n)``: the identity matrix of size :math:`n\times n`. `n` should be an
integer expression. For instance ``Id(meshdim)`` is allowed.
- - ``qdim(u)``: the dimension of the variable ``u`` (i.e. the size for fixed
size variables and the dimension of the vector field for f.e.m. variables)
+ - ``qdim(u)``: the total dimension of the variable ``u`` (i.e. the size for
fixed size variables and the total dimension of the vector/tensor field for
f.e.m. variables)
+ - ``qdims(u)``: the dimensions of the variable ``u`` (i.e. the size for
fixed size variables and the vector of dimensions of the vector/tensor field
for f.e.m. variables)
Special expressions linked to the current position
**************************************************
Modified: trunk/getfem/src/getfem/getfem_generic_assembly.h
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_generic_assembly.h?rev=4535&r1=4534&r2=4535&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_generic_assembly.h (original)
+++ trunk/getfem/src/getfem/getfem_generic_assembly.h Fri Mar 14 15:46:33 2014
@@ -312,6 +312,21 @@
return mf ? associated_mf(name)->get_qdim() * (n / ndof) : n;
}
+ bgeot::multi_index qdims(const std::string &name) const {
+ const mesh_fem *mf = associated_mf(name);
+ size_type n = gmm::vect_size(value(name));
+ if (mf) {
+ bgeot::multi_index mi = mf->get_qdims();
+ size_type qmult = n / mf->nb_dof();
+ if (qmult > 1) {
+ if (mi.back() == 1) mi.back() *= qmult; else mi.push_back(qmult);
+ }
+ return mi;
+ } else {
+ bgeot::multi_index mi(1); mi[0] = n; return mi;
+ }
+ }
+
const model_real_plain_vector &value(const std::string &name) const {
if (model)
return model->real_variable(name);
Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=4535&r1=4534&r2=4535&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Fri Mar 14 15:46:33 2014
@@ -769,6 +769,7 @@
case 0:
str << (nt ? scalar_type(0) : pnode->t[0]);
break;
+
case 1:
str << "[";
for (size_type i = 0; i < pnode->tensor_proper_size(0); ++i) {
@@ -777,6 +778,7 @@
}
str << "]";
break;
+
case 2:
str << "[";
for (size_type i = 0; i < pnode->tensor_proper_size(0); ++i) {
@@ -788,6 +790,7 @@
}
str << "]";
break;
+
case 3:
str << "[";
for (size_type i = 0; i < pnode->tensor_proper_size(0); ++i) {
@@ -802,6 +805,7 @@
}
str << "]";
break;
+
case 4:
str << "[";
for (size_type i = 0; i < pnode->tensor_proper_size(0); ++i) {
@@ -819,6 +823,21 @@
}
str << "]";
break;
+
+ case 5: case 6:
+ str << "Reshape([";
+ for (size_type i = 0; i < pnode->tensor_proper_size(); ++i) {
+ if (i != 0) str << "; ";
+ str << (nt ? scalar_type(0) : pnode->t[i]);
+ }
+ 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");
}
GMM_ASSERT1(pnode->children.size() == 0, "Invalid tree");
@@ -918,9 +937,8 @@
if (pnode->qdim1 == 1)
str << "*Test_" << pnode->name_test1;
else {
- str << "*([ 0";
- for (size_type i = 1; i < pnode->qdim1; ++i) str << ", 0";
- str << "].Test_" << pnode->name_test1 << ")";
+ str << "*(Reshape(Test_" << pnode->name_test1 << ","
+ << pnode->qdim1<< ")(1))";
}
}
if (pnode->name_test2.size()) {
@@ -928,9 +946,8 @@
if (pnode->qdim2 == 1)
str << "*Test2_" << pnode->name_test2;
else {
- str << "*([ 0";
- for (size_type i = 1; i < pnode->qdim2; ++i) str << ", 0";
- str << "].Test2_" << pnode->name_test2 << ")";
+ str << "*(Reshape(Test2_" << pnode->name_test2 << ","
+ << pnode->qdim2<< ")(1))";
}
}
if (pnode->test_function_type) str << ")";
@@ -1932,6 +1949,7 @@
SPEC_FUNCTIONS.insert("pi");
SPEC_FUNCTIONS.insert("meshdim");
SPEC_FUNCTIONS.insert("qdim");
+ SPEC_FUNCTIONS.insert("qdims");
SPEC_FUNCTIONS.insert("Id");
// Predefined operators
@@ -4048,7 +4066,7 @@
pnode->t = child0->t;
gmm::scale(pnode->t.as_vector(), scalar_type(child1->t[0]));
} else {
- if (dim0+dim1 > 4)
+ if (dim0+dim1 > 6)
ga_throw_error(expr, pnode->pos, "Unauthorized "
"tensor multiplication.");
for (size_type i = 0; i < dim0; ++i)
@@ -4075,7 +4093,7 @@
for (size_type i = 0; i < dim0; ++i)
mi.push_back(child0->tensor_proper_size(i));
} else {
- if (dim0+dim1 > 4)
+ if (dim0+dim1 > 6)
ga_throw_error(expr, pnode->pos, "Unauthorized "
"tensor multiplication.");
for (size_type i = 0; i < dim0; ++i)
@@ -4441,35 +4459,30 @@
}
} else {
size_type q = workspace.qdim(name), n = mf->linked_mesh().dim();
+ bgeot::multi_index mii = workspace.qdims(name);
if (!q) ga_throw_error(expr, pnode->pos,
"Invalid null size of variable");
switch (val_grad_or_hess) {
case 0: // value
pnode->node_type = GA_NODE_VAL;
- if (q == 1)
- pnode->init_scalar_tensor(scalar_type(0));
- else
- pnode->init_vector_tensor(q);
break;
case 1: // grad
pnode->node_type = GA_NODE_GRAD;
- if (q == 1 && n == 1)
- pnode->init_scalar_tensor(scalar_type(0));
- else if (q == 1)
- pnode->init_vector_tensor(n);
- else
- pnode->init_matrix_tensor(q, n);
+ if (n > 1) {
+ if (q == 1 && mii.size() <= 1) mii[0] = n;
+ else mii.push_back(n);
+ }
break;
case 2: // Hessian
pnode->node_type = GA_NODE_HESS;
- if (q == 1 && n == 1)
- pnode->init_scalar_tensor(scalar_type(0));
- else if (q == 1)
- pnode->init_matrix_tensor(n,n);
- else
- pnode->init_third_order_tensor(q, n, n);
+ if (n > 1) {
+ if (q == 1 && mii.size() <= 1) { mii[0] = n;
mii.push_back(n); }
+ else { mii.push_back(n); mii.push_back(n); }
+ }
break;
}
+ pnode->t.adjust_sizes(mii);
+ pnode->test_function_type = 0;
}
} else {
if (workspace.is_constant(name))
@@ -4512,36 +4525,48 @@
}
} else {
size_type q = workspace.qdim(name), n = mf->linked_mesh().dim();
+ bgeot::multi_index mii = workspace.qdims(name);
+ if (mii.size() > 6)
+ ga_throw_error(expr, pnode->pos,
+ "Tensor with too much dimensions. Limited to
6");
if (!q)
ga_throw_error(expr, pnode->pos,
"Invalid null size of variable");
switch (val_grad_or_hess) {
case 0: // value
pnode->node_type = GA_NODE_TEST;
- if (q == 1)
+ if (q == 1 && mii.size() <= 1)
pnode->init_vector_tensor(2);
- else
- pnode->init_matrix_tensor(2,q);
+ else {
+ mii.insert(mii.begin(), 2);
+ pnode->t.adjust_sizes(mii);
+ }
pnode->test_function_type = test;
break;
case 1: // grad
pnode->node_type = GA_NODE_GRAD_TEST;
- if (q == 1 && n == 1)
+ if (q == 1 && mii.size() <= 1 && n == 1)
pnode->init_vector_tensor(2);
- else if (q == 1)
- pnode->init_matrix_tensor(2,n);
- else
- pnode->init_third_order_tensor(2,q,n);
+ else if (q == 1 && mii.size() <= 1)
+ pnode->init_matrix_tensor(2, n);
+ else {
+ mii.insert(mii.begin(), 2);
+ if (n > 1) mii.push_back(n);
+ pnode->t.adjust_sizes(mii);
+ }
pnode->test_function_type = test;
break;
case 2: // hessian
pnode->node_type = GA_NODE_HESS_TEST;
- if (q == 1 && n == 1)
+ if (q == 1 && mii.size() <= 1 && n == 1)
pnode->init_vector_tensor(2);
- else if (q == 1)
+ else if (q == 1 && mii.size() <= 1)
pnode->init_third_order_tensor(2,n,n);
- else
- pnode->init_fourth_order_tensor(2,q,n,n);
+ else {
+ mii.insert(mii.begin(), 2);
+ if (n > 1) { mii.push_back(n); mii.push_back(n); }
+ pnode->t.adjust_sizes(mii);
+ }
pnode->test_function_type = test;
break;
}
@@ -4570,7 +4595,7 @@
if (pnode->children.size() < 3)
ga_throw_error(expr, child1->pos,
"Not enough parameters for Reshape");
- if (pnode->children.size() > 6)
+ if (pnode->children.size() > 8)
ga_throw_error(expr, child1->pos,
"Too many parameters for Reshape");
pnode->t = child1->t;
@@ -4693,6 +4718,25 @@
if (pnode->t[0] <= 0)
ga_throw_error(expr, pnode->pos,
"Invalid null size of variable");
+ } else if (!(child0->name.compare("qdims"))) {
+ if (child1->node_type != GA_NODE_VAL)
+ ga_throw_error(expr, pnode->pos, "The argument of qdim "
+ "function can only be a variable name.");
+ pnode->node_type = GA_NODE_CONSTANT;
+ bgeot::multi_index mii = workspace.qdims(child1->name);
+ if (mii.size() > 6)
+ ga_throw_error(expr, pnode->pos,
+ "Tensor with too much dimensions. Limited to 6");
+ if (mii.size() == 0 || scalar_type(mii[0]) <= 0)
+ ga_throw_error(expr, pnode->pos,
+ "Invalid null size of variable");
+ if (mii.size() == 1)
+ pnode->init_scalar_tensor(scalar_type(mii[0]));
+ if (mii.size() >= 1) {
+ pnode->init_vector_tensor(mii.size());
+ for (size_type i = 0; i < mii.size(); ++i)
+ pnode->t[i] = scalar_type(mii[i]);
+ }
} else if (!(child0->name.compare("Id"))) {
bool valid = (child1->node_type == GA_NODE_CONSTANT);
int n = valid ? int(round(child1->t[0])) : -1;
Modified: trunk/getfem/tests/test_assembly.cc
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/tests/test_assembly.cc?rev=4535&r1=4534&r2=4535&view=diff
==============================================================================
--- trunk/getfem/tests/test_assembly.cc (original)
+++ trunk/getfem/tests/test_assembly.cc Fri Mar 14 15:46:33 2014
@@ -1046,8 +1046,6 @@
VEC_TEST_1("Test for source term", ndofu, "u.Test_u", mim, size_type(-1),
Iu, getfem::asm_source_term(V, mim, mf_u, mf_u, U));
-
-
}
if (all) {
@@ -1065,14 +1063,13 @@
if (N == 2)
{VEC_TEST_1("Test for Neumann term", ndofu,
- "([A(1), A(3); A(2), A(4)]'*Normal).Test_u", mim,
+ "(A'*Normal).Test_u", mim,
NEUMANN_BOUNDARY_NUM,
Iu, getfem::asm_normal_source_term(V, mim, mf_u, mf_u,
A, NEUMANN_BOUNDARY_NUM));}
if (N == 3)
{VEC_TEST_1("Test for Neumann term", ndofu,
- "([A(1), A(4), A(7); A(2), A(5), A(8); A(3), A(6), A(9)]'"
- "*Normal).Test_u", mim, NEUMANN_BOUNDARY_NUM,
+ "(A'*Normal).Test_u", mim, NEUMANN_BOUNDARY_NUM,
Iu, getfem::asm_normal_source_term(V, mim, mf_u, mf_u,
A, NEUMANN_BOUNDARY_NUM));}
}
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] r4535 - in /trunk/getfem: ./ doc/sphinx/source/userdoc/ src/ src/getfem/ tests/,
Yves . Renard <=