[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Sun, 24 Mar 2019 14:24:28 -0400 (EDT) |
branch: master
commit 3553430cc06875be4e242a6a5031fbb22203f9f9
Author: Konstantinos Poulios <address@hidden>
Date: Sun Mar 24 19:24:18 2019 +0100
Add support for variables defined at integration points
---
src/getfem/getfem_generic_assembly.h | 3 +
src/getfem/getfem_generic_assembly_tree.h | 11 +++
src/getfem/getfem_models.h | 6 +-
src/getfem_generic_assembly_compile_and_exec.cc | 59 ++++++++++--
src/getfem_generic_assembly_semantic.cc | 102 +++++++++++++--------
src/getfem_generic_assembly_workspace.cc | 6 ++
src/getfem_models.cc | 10 +++
tests/Makefile.am | 4 +
tests/test_internal_variables.cc | 114 ++++++++++++++++++++++++
tests/test_internal_variables.pl | 45 ++++++++++
10 files changed, 314 insertions(+), 46 deletions(-)
diff --git a/src/getfem/getfem_generic_assembly.h
b/src/getfem/getfem_generic_assembly.h
index ece421f..664bdcb 100644
--- a/src/getfem/getfem_generic_assembly.h
+++ b/src/getfem/getfem_generic_assembly.h
@@ -431,6 +431,9 @@ namespace getfem {
void add_fem_variable(const std::string &name, const mesh_fem &mf,
const gmm::sub_interval &I,
const model_real_plain_vector &VV);
+ void add_im_variable(const std::string &name, const im_data &imd,
+ const gmm::sub_interval &I,
+ const model_real_plain_vector &VV);
void add_fixed_size_variable(const std::string &name,
const gmm::sub_interval &I,
const model_real_plain_vector &VV);
diff --git a/src/getfem/getfem_generic_assembly_tree.h
b/src/getfem/getfem_generic_assembly_tree.h
index f7628d6..3be591a 100644
--- a/src/getfem/getfem_generic_assembly_tree.h
+++ b/src/getfem/getfem_generic_assembly_tree.h
@@ -277,6 +277,14 @@ namespace getfem {
void init_matrix_tensor(size_type n, size_type m)
{ set_to_original(); t.adjust_sizes(n, m); }
+ void init_identity_matrix_tensor(size_type n) {
+ init_matrix_tensor(n, n);
+ auto itw = t.begin();
+ for (size_type i = 0; i < n; ++i)
+ for (size_type j = 0; j < n; ++j)
+ *itw++ = (i == j) ? scalar_type(1) : scalar_type(0);
+ }
+
void init_third_order_tensor(size_type n, size_type m, size_type l)
{ set_to_original(); t.adjust_sizes(n, m, l); }
@@ -369,6 +377,9 @@ namespace getfem {
inline void init_matrix_tensor(size_type n, size_type m)
{ t.init_matrix_tensor(n, m); test_function_type = 0; }
+ inline void init_identity_matrix_tensor(size_type n)
+ { t.init_identity_matrix_tensor(n); test_function_type = 0; }
+
inline void init_third_order_tensor(size_type n, size_type m, size_type l)
{ t.init_third_order_tensor(n, m, l); test_function_type = 0; }
diff --git a/src/getfem/getfem_models.h b/src/getfem/getfem_models.h
index c77fca5..49f9e87 100644
--- a/src/getfem/getfem_models.h
+++ b/src/getfem/getfem_models.h
@@ -750,7 +750,11 @@ namespace getfem {
}
- /** Add data, defined at integration points.*/
+ /** Add variable defined at integration points. */
+ void add_im_variable(const std::string &name, const im_data &im_data,
+ size_type niter = 1);
+
+ /** Add data defined at integration points. */
void add_im_data(const std::string &name, const im_data &im_data,
size_type niter = 1);
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index ca6ecc4..f07c821 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -4125,6 +4125,32 @@ namespace getfem {
interpolate(interpolate_) {}
};
+ struct ga_instruction_imd_vector_assembly : public ga_instruction {
+ const base_tensor &t;
+ base_vector &V;
+ const fem_interpolation_context &ctx;
+ const gmm::sub_interval &I;
+ const im_data *imd;
+ scalar_type &coeff;
+ const size_type &ipt;
+ virtual int exec() {
+ GA_DEBUG_INFO("Instruction: vector term assembly for im_data variable");
+ size_type ifirst = I.first();
+ size_type i = t.size()
+ * imd->filtered_index_of_point(ctx.convex_num(), ipt);
+ GMM_ASSERT1(i+t.size() <= I.size(), "Internal error
"<<i<<"+"<<t.size()<<" <= "<<I.size());
+ for (const auto &val : t.as_vector())
+ V[ifirst+(i++)] += coeff*val;
+ return 0;
+ }
+ ga_instruction_imd_vector_assembly
+ (const base_tensor &t_, base_vector &V_,
+ const fem_interpolation_context &ctx_, const gmm::sub_interval &I_,
+ const im_data *imd_, scalar_type &coeff_, const size_type &ipt_)
+ : t(t_), V(V_), ctx(ctx_), I(I_), imd(imd_), coeff(coeff_), ipt(ipt_)
+ {}
+ };
+
struct ga_instruction_vector_assembly : public ga_instruction {
const base_tensor &t;
base_vector &V;
@@ -4305,6 +4331,7 @@ namespace getfem {
const fem_interpolation_context &ctx1, &ctx2;
const gmm::sub_interval &Ir1, &Ir2, &In1, &In2;
const mesh_fem *mfn1, *mfn2, **mfg1, **mfg2;
+ const im_data *imd1, *imd2;
const scalar_type &coeff, &alpha1, &alpha2;
const size_type &nbpt, &ipt;
base_vector elem;
@@ -4313,7 +4340,7 @@ namespace getfem {
virtual int exec() {
GA_DEBUG_INFO("Instruction: matrix term assembly");
bool empty_weight = (coeff == scalar_type(0));
- if (ipt == 0 || interpolate) {
+ if (ipt == 0 || interpolate || imd1 || imd2) { // initialize
if (empty_weight) elem.resize(0);
elem.resize(t.size());
if (!empty_weight)
@@ -4323,7 +4350,7 @@ namespace getfem {
// Faster than a daxpy blas call on my config
add_scaled_4(t, coeff*alpha1*alpha2, elem);
- if (ipt == nbpt-1 || interpolate) { // finalize
+ if (ipt == nbpt-1 || interpolate || imd1 || imd2) { // finalize
const mesh_fem *pmf1 = mfg1 ? *mfg1 : mfn1;
const mesh_fem *pmf2 = mfg2 ? *mfg2 : mfn2;
bool reduced = (pmf1 && pmf1->is_reduced())
@@ -4337,11 +4364,12 @@ namespace getfem {
if (ninf == scalar_type(0)) return 0;
size_type s1 = t.sizes()[0], s2 = t.sizes()[1];
- size_type cv1 = pmf1 ? ctx1.convex_num() : s1;
- size_type cv2 = pmf2 ? ctx2.convex_num() : s2;
+ size_type cv1 = ctx1.convex_num(), cv2 = ctx2.convex_num();
size_type N = 1;
size_type ifirst1 = I1.first(), ifirst2 = I2.first();
+ if (imd1) ifirst1 += s1*imd1->filtered_index_of_point(cv1, ipt);
+ if (imd2) ifirst2 += s2*imd2->filtered_index_of_point(cv2, ipt);
if (pmf1) {
if (!ctx1.is_convex_num_valid()) return 0;
@@ -4353,7 +4381,7 @@ namespace getfem {
} else
populate_contiguous_dofs_vector(dofs1, s1, ifirst1); // --> dofs1
- if (pmf1 == pmf2 && cv1 == cv2) {
+ if (pmf1 == pmf2 && (pmf1 ? (cv1 == cv2) : (s1 == s2))) {
if (ifirst1 == ifirst2) {
add_elem_matrix(K, dofs1, dofs1, dofs1_sort, elem, ninf*1E-14, N);
} else {
@@ -4382,14 +4410,14 @@ namespace getfem {
const fem_interpolation_context &ctx2_,
const gmm::sub_interval &Ir1_, const gmm::sub_interval &In1_,
const gmm::sub_interval &Ir2_, const gmm::sub_interval &In2_,
- const mesh_fem *mfn1_, const mesh_fem **mfg1_,
- const mesh_fem *mfn2_, const mesh_fem **mfg2_,
+ const mesh_fem *mfn1_, const mesh_fem **mfg1_, const im_data *imd1_,
+ const mesh_fem *mfn2_, const mesh_fem **mfg2_, const im_data *imd2_,
const scalar_type &coeff_, const scalar_type &a1, const scalar_type &a2,
const size_type &nbpt_, const size_type &ipt_, bool interpolate_)
: t(t_), Kr(Kr_), Kn(Kn_), ctx1(ctx1_), ctx2(ctx2_),
Ir1(Ir1_), Ir2(Ir2_), In1(In1_), In2(In2_),
mfn1(mfn1_), mfn2(mfn2_), mfg1(mfg1_), mfg2(mfg2_),
- coeff(coeff_), alpha1(a1), alpha2(a2),
+ imd1(imd1_), imd2(imd2_), coeff(coeff_), alpha1(a1), alpha2(a2),
nbpt(nbpt_), ipt(ipt_), interpolate(interpolate_),
dofs1(0), dofs2(0) {}
};
@@ -6771,6 +6799,8 @@ namespace getfem {
"weak form has to be a scalar quantity");
const mesh_fem *mf
= workspace.associated_mf(root->name_test1);
+ const im_data *imd
+ = workspace.associated_im_data(root->name_test1);
add_interval_to_gis(workspace, root->name_test1, gis);
if (mf) {
@@ -6802,6 +6832,14 @@ namespace getfem {
(root->tensor(), workspace.unreduced_vector(),
workspace.assembled_vector(), ctx, *Ir, *In, mf, mfg,
gis.coeff, gis.nbpt, gis.ipt, interpolate);
+ } else if (imd) {
+ GMM_ASSERT1(root->interpolate_name_test1.size() == 0,
+ "Interpolate transformation on integration "
+ "point variable");
+ pgai = std::make_shared<ga_instruction_imd_vector_assembly>
+ (root->tensor(), workspace.assembled_vector(), gis.ctx,
+ workspace.interval_of_variable(root->name_test1),
+ imd, gis.coeff, gis.ipt);
} else {
pgai = std::make_shared<ga_instruction_vector_assembly>
(root->tensor(), workspace.assembled_vector(),
@@ -6818,6 +6856,9 @@ namespace getfem {
const mesh_fem
*mf1=workspace.associated_mf(root->name_test1),
*mf2=workspace.associated_mf(root->name_test2),
**mfg1 = 0, **mfg2 = 0;
+ const im_data
+ *imd1 = workspace.associated_im_data(root->name_test1),
+ *imd2 = workspace.associated_im_data(root->name_test2);
const std::string &intn1 = root->interpolate_name_test1,
&intn2 = root->interpolate_name_test2;
bool secondary1 = intn1.size() &&
@@ -6907,7 +6948,7 @@ namespace getfem {
pgai = std::make_shared<ga_instruction_matrix_assembly>
(root->tensor(), workspace.unreduced_matrix(),
workspace.assembled_matrix(), ctx1, ctx2,
- *Ir1, *In1, *Ir2, *In2, mf1, mfg1, mf2, mfg2,
+ *Ir1, *In1, *Ir2, *In2, mf1, mfg1, imd1, mf2, mfg2,
imd2,
gis.coeff, *alpha1, *alpha2, gis.nbpt, gis.ipt,
interpolate);
}
diff --git a/src/getfem_generic_assembly_semantic.cc
b/src/getfem_generic_assembly_semantic.cc
index 09cc7e5..cf27450 100644
--- a/src/getfem_generic_assembly_semantic.cc
+++ b/src/getfem_generic_assembly_semantic.cc
@@ -471,13 +471,15 @@ namespace getfem {
case GA_NODE_XFEM_MINUS_HESS_TEST: case GA_NODE_XFEM_MINUS_DIVERG_TEST:
{
const mesh_fem *mf = workspace.associated_mf(pnode->name);
+ const im_data *imd = workspace.associated_im_data(pnode->name);
size_type t_type = pnode->test_function_type;
if (t_type == 1) {
pnode->name_test1 = pnode->name;
pnode->interpolate_name_test1 = pnode->interpolate_name;
pnode->interpolate_name_test2 = pnode->name_test2 = "";
- pnode->qdim1 = (mf ? workspace.qdim(pnode->name)
- : gmm::vect_size(workspace.value(pnode->name)));
+ pnode->qdim1 = (mf || imd)
+ ? workspace.qdim(pnode->name)
+ : gmm::vect_size(workspace.value(pnode->name));
if (option == 1)
workspace.test1.insert
(var_trans_pair(pnode->name_test1,
@@ -489,8 +491,9 @@ namespace getfem {
pnode->interpolate_name_test1 = pnode->name_test1 = "";
pnode->name_test2 = pnode->name;
pnode->interpolate_name_test2 = pnode->interpolate_name;
- pnode->qdim2 = (mf ? workspace.qdim(pnode->name)
- : gmm::vect_size(workspace.value(pnode->name)));
+ pnode->qdim2 = (mf || imd)
+ ? workspace.qdim(pnode->name)
+ : gmm::vect_size(workspace.value(pnode->name));
if (option == 1)
workspace.test2.insert
(var_trans_pair(pnode->name_test2,
@@ -499,23 +502,32 @@ namespace getfem {
ga_throw_error(pnode->expr, pnode->pos,
"Invalid null size of variable");
}
- if (!mf) {
- size_type n = workspace.qdim(pnode->name);
- if (!n)
- ga_throw_error(pnode->expr, pnode->pos,
- "Invalid null size of variable");
- if (n == 1) {
+ size_type q = workspace.qdim(pnode->name);
+ if (!q)
+ ga_throw_error(pnode->expr, pnode->pos,
+ "Invalid null size of variable");
+ if (!mf & !imd) { // global variable
+ if (q == 1) {
+ pnode->init_vector_tensor(1);
+ pnode->tensor()[0] = scalar_type(1);
+ } else
+ pnode->init_identity_matrix_tensor(q);
+ pnode->test_function_type = t_type;
+ } else if (imd) {
+ bgeot::multi_index mii = workspace.qdims(pnode->name);
+ if (q == 1 && mii.size() <= 1) {
pnode->init_vector_tensor(1);
pnode->tensor()[0] = scalar_type(1);
- pnode->test_function_type = t_type;
} else {
- pnode->init_matrix_tensor(n,n);
- pnode->test_function_type = t_type;
- for (size_type i = 0; i < n; ++i)
- for (size_type j = 0; j < n; ++j)
- pnode->tensor()(i,j) = (i==j) ? scalar_type(1)
- : scalar_type(0);
+ mii.insert(mii.begin(), q);
+ pnode->t.set_to_original();
+ pnode->t.adjust_sizes(mii);
+ auto itw = pnode->tensor().begin();
+ for (size_type i = 0; i < q; ++i) // set identity matrix
+ for (size_type j = 0; j < q; ++j)
+ *itw++ = (i == j) ? scalar_type(1) : scalar_type(0);
}
+ pnode->test_function_type = t_type;
}
}
break;
@@ -1731,8 +1743,8 @@ namespace getfem {
workspace.test1.insert
(var_trans_pair(pnode->name_test1,
pnode->interpolate_name_test1));
- pnode->qdim1 = mf ? workspace.qdim(name)
- : gmm::vect_size(workspace.value(name));
+ pnode->qdim1 = (mf || imd) ? workspace.qdim(name)
+ : gmm::vect_size(workspace.value(name));
if (!(pnode->qdim1))
ga_throw_error(pnode->expr, pnode->pos,
"Invalid null size of variable");
@@ -1743,14 +1755,14 @@ namespace getfem {
workspace.test2.insert
(var_trans_pair(pnode->name_test2,
pnode->interpolate_name_test2));
- pnode->qdim2 = mf ? workspace.qdim(name)
- : gmm::vect_size(workspace.value(name));
+ pnode->qdim2 = (mf || imd) ? workspace.qdim(name)
+ : gmm::vect_size(workspace.value(name));
if (!(pnode->qdim2))
ga_throw_error(pnode->expr, pnode->pos,
"Invalid null size of variable");
}
- if (!mf && (test || !imd)) {
+ if (!mf && !imd) { // global variable
if (prefix_id)
ga_throw_error(pnode->expr, pnode->pos, "Gradient, Hessian or "
"Divergence cannot be evaluated for fixed size data.");
@@ -1761,8 +1773,8 @@ namespace getfem {
else
pnode->node_type = GA_NODE_VAL;
- size_type n = gmm::vect_size(workspace.value(name));
- if (n == 1) {
+ size_type q = gmm::vect_size(workspace.value(name));
+ if (q == 1) {
if (test) {
pnode->init_vector_tensor(1);
pnode->tensor()[0] = scalar_type(1);
@@ -1771,23 +1783,43 @@ namespace getfem {
pnode->init_scalar_tensor(workspace.value(name)[0]);
} else {
if (test) {
- pnode->init_matrix_tensor(n,n);
- for (size_type i = 0; i < n; ++i)
- for (size_type j = 0; j < n; ++j)
- pnode->tensor()(i,j) = (i == j) ? scalar_type(1)
- : scalar_type(0);
+ pnode->init_identity_matrix_tensor(q);
} else {
pnode->t.adjust_sizes(workspace.qdims(name));
gmm::copy(workspace.value(name), pnode->tensor().as_vector());
}
}
- } else if (!test && imd) {
+ } else if (imd) { // im_data variable
+ size_type q = workspace.qdim(name);
+ bgeot::multi_index mii = workspace.qdims(name);
+
+ if (!q) ga_throw_error(pnode->expr, pnode->pos,
+ "Invalid null size of variable " << name);
+ if (mii.size() > 6)
+ ga_throw_error(pnode->expr, pnode->pos,
+ "Tensor with too many dimensions. Limited to 6");
if (prefix_id)
ga_throw_error(pnode->expr, pnode->pos, "Gradient, Hessian or "
"Divergence cannot be evaluated for im data.");
- pnode->node_type = GA_NODE_VAL;
- pnode->t.adjust_sizes(workspace.qdims(name));
- } else {
+
+ pnode->node_type = test ? GA_NODE_VAL_TEST : GA_NODE_VAL;
+
+ if (test) {
+ if (q == 1 && mii.size() <= 1) {
+ pnode->init_vector_tensor(1);
+ pnode->tensor()[0] = scalar_type(1);
+ } else {
+ mii.insert(mii.begin(), q);
+ pnode->t.set_to_original();
+ pnode->t.adjust_sizes(mii);
+ auto itw = pnode->tensor().begin();
+ for (size_type i = 0; i < q; ++i) // set identity matrix
+ for (size_type j = 0; j < q; ++j)
+ *itw++ = (i == j) ? scalar_type(1) : scalar_type(0);
+ }
+ } else
+ pnode->t.adjust_sizes(mii);
+ } else { // mesh_fem variable
size_type q = workspace.qdim(name);
size_type n = mf->linked_mesh().dim();
bgeot::multi_index mii = workspace.qdims(name);
@@ -1806,10 +1838,8 @@ namespace getfem {
if (test && q == 1 && mii.size() <= 1) {
mii.resize(1);
mii[0] = 2;
- } else if (test) {
+ } else if (test)
mii.insert(mii.begin(), 2);
- pnode->t.adjust_sizes(mii);
- }
break;
case 1: // grad
pnode->node_type = test ? GA_NODE_GRAD_TEST : GA_NODE_GRAD;
diff --git a/src/getfem_generic_assembly_workspace.cc
b/src/getfem_generic_assembly_workspace.cc
index 2e66c16..6ad3bc0 100644
--- a/src/getfem_generic_assembly_workspace.cc
+++ b/src/getfem_generic_assembly_workspace.cc
@@ -48,6 +48,12 @@ namespace getfem {
variables.emplace(name, var_description(true, &mf, 0, I, &VV, 1));
}
+ void ga_workspace::add_im_variable
+ (const std::string &name, const im_data &imd,
+ const gmm::sub_interval &I, const model_real_plain_vector &VV) {
+ variables.emplace(name, var_description(true, 0, &imd, I, &VV, 1));
+ }
+
void ga_workspace::add_fixed_size_variable
(const std::string &name,
const gmm::sub_interval &I, const model_real_plain_vector &VV) {
diff --git a/src/getfem_models.cc b/src/getfem_models.cc
index 0aaa403..6cd0103 100644
--- a/src/getfem_models.cc
+++ b/src/getfem_models.cc
@@ -753,6 +753,16 @@ namespace getfem {
gmm::copy(t.as_vector(), set_complex_variable(name));
}
+ void model::add_im_variable(const std::string &name, const im_data &imd,
+ size_type niter) {
+ check_name_validity(name);
+ variables.emplace(name,
+ var_description(true, is_complex(), 0, &imd, niter));
+ variables[name].set_size();
+ add_dependency(imd);
+ act_size_to_be_done = true;
+ }
+
void model::add_im_data(const std::string &name, const im_data &imd,
size_type niter) {
check_name_validity(name);
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 385006a..76eb7b2 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -44,6 +44,7 @@ check_PROGRAMS = \
test_assembly \
test_assembly_assignment \
test_interpolated_fem \
+ test_internal_variables \
test_range_basis \
laplacian \
laplacian_with_bricks \
@@ -98,6 +99,7 @@ test_mesh_SOURCES = test_mesh.cc
geo_trans_inv_SOURCES = geo_trans_inv.cc
test_int_set_SOURCES = test_int_set.cc
test_interpolated_fem_SOURCES = test_interpolated_fem.cc
+test_internal_variables_SOURCES = test_internal_variables.cc
test_tree_sorted_SOURCES = test_tree_sorted.cc
test_mat_elem_SOURCES = test_mat_elem.cc
test_slice_SOURCES = test_slice.cc
@@ -138,6 +140,7 @@ TESTS = \
test_assembly.pl \
test_assembly_assignment.pl \
test_interpolated_fem.pl \
+ test_internal_variables.pl \
test_range_basis.pl \
laplacian.pl \
laplacian_with_bricks.pl \
@@ -178,6 +181,7 @@ EXTRA_DIST =
\
geo_trans_inv.pl \
test_int_set.pl \
test_interpolated_fem.pl \
+ test_internal_variables.pl \
test_slice.pl \
test_mesh_im_level_set.pl \
thermo_elasticity_electrical_coupling.pl \
diff --git a/tests/test_internal_variables.cc b/tests/test_internal_variables.cc
new file mode 100644
index 0000000..7d0e409
--- /dev/null
+++ b/tests/test_internal_variables.cc
@@ -0,0 +1,114 @@
+/*===========================================================================
+
+ Copyright (C) 2019-2019 Konstantinos Poulios.
+
+ This file is a part of GetFEM++
+
+ GetFEM++ is free software; you can redistribute it and/or modify it
+ under the terms of the GNU Lesser General Public License as published
+ by the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version along with the GCC Runtime Library
+ Exception either version 3.1 or (at your option) any later version.
+ This program is distributed in the hope that it will be useful, but
+ WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+ License and GCC Runtime Library Exception for more details.
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, write to the Free Software Foundation,
+ Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
+
+===========================================================================*/
+#include "getfem/getfem_regular_meshes.h"
+#include "getfem/getfem_model_solvers.h"
+#include "getfem/getfem_generic_assembly.h"
+
+using bgeot::dim_type;
+using bgeot::size_type;
+using bgeot::scalar_type;
+using bgeot::base_node;
+
+int main(int argc, char *argv[]) {
+
+// gmm::set_traces_level(1);
+
+ bgeot::md_param PARAM;
+ PARAM.add_int_param("NX", 3);
+ PARAM.add_int_param("NY", 2);
+ PARAM.add_int_param("FEM_ORDER", 1);
+ PARAM.add_int_param("IM_ORDER", 1);
+ PARAM.add_int_param("DIFFICULTY", 0);
+ PARAM.read_command_line(argc, argv);
+ size_type NX = PARAM.int_value("NX", "Number of elements in X direction");
+ size_type NY = PARAM.int_value("NY", "Number of elements in Y direction");
+ dim_type FEM_ORDER = dim_type(PARAM.int_value("FEM_ORDER", "Degree of finite
element basis"));
+ dim_type IM_ORDER = dim_type(PARAM.int_value("IM_ORDER", "Degree of
integration method"));
+ size_type DIFFICULTY = PARAM.int_value("DIFFICULTY", "Difficulty of test (0
or 1)");
+
+ getfem::mesh m;
+ getfem::regular_unit_mesh(m, {NX, NY},
bgeot::geometric_trans_descriptor("GT_QK(2, 2)"));
+
+ getfem::mesh_region outer_faces;
+ getfem::outer_faces_of_mesh(m, outer_faces);
+ m.region(100) = getfem::select_faces_of_normal(m, outer_faces, base_node(-1,
0), 0.001);
+ m.region(101) = getfem::select_faces_of_normal(m, outer_faces, base_node(1,
0), 0.001);
+ m.region(102) = getfem::mesh_region::merge(m.region(100), m.region(101));
+
+ dim_type N(2);
+ getfem::mesh_fem mf(m, N), mf_intern(m);
+ mf.set_classical_finite_element(FEM_ORDER);
+ if (DIFFICULTY) mf_intern.set_qdim(3,4);
+ mf_intern.set_classical_discontinuous_finite_element(IM_ORDER);
+
+ getfem::mesh_im mim(m);
+ mim.set_integration_method(dim_type(2*IM_ORDER+1));
+
+ getfem::im_data mimd(mim);
+ if (DIFFICULTY) mimd.set_tensor_size(bgeot::multi_index(3,4));
+
+ getfem::model md1, md2;
+ md1.add_fem_variable("u", mf);
+ md2.add_fem_variable("u", mf);
+ md1.add_im_variable("p", mimd);
+ md2.add_fem_variable("p", mf_intern);
+
+ md1.add_initialized_scalar_data("G", 1);
+ md2.add_initialized_scalar_data("G", 1);
+ md1.add_initialized_scalar_data("K", 1);
+ md2.add_initialized_scalar_data("K", 1);
+
+ std::string exprA, exprB;
+ if (DIFFICULTY) {
+ exprA =
"(-1e3*asin(p(1,1))*Id(2)+2*G*(Sym(Grad_u)-Div_u*Id(2)/3)):Grad_Test_u";
+ exprB =
"(p+sin(0.001*K*Trace(Sym(Grad_u)))*[1,1,1,1;1,1,1,1;1,1,1,1]):Test_p";
+ } else {
+ exprA = "(-p*Id(2)+2*G*(Sym(Grad_u)-Div_u*Id(2)/3)):Grad_Test_u";
+ exprB = "(p+K*Trace(Sym(Grad_u)))*Test_p";
+ }
+ getfem::add_nonlinear_generic_assembly_brick(md1, mim, exprA);
+ getfem::add_nonlinear_generic_assembly_brick(md2, mim, exprA);
+ getfem::add_nonlinear_generic_assembly_brick(md1, mim, exprB);
+ getfem::add_nonlinear_generic_assembly_brick(md2, mim, exprB);
+
+ md1.add_filtered_fem_variable("dirmult", mf, 102);
+ md2.add_filtered_fem_variable("dirmult", mf, 102);
+ getfem::add_linear_generic_assembly_brick(md1, mim,
"(u-0.001*X(1)*[1;0]).dirmult", 102);
+ getfem::add_linear_generic_assembly_brick(md2, mim,
"(u-0.001*X(1)*[1;0]).dirmult", 102);
+
+ gmm::iteration iter(1E-9, 1, 100);
+ getfem::standard_solve(md1, iter);
+ iter.init();
+ getfem::standard_solve(md2, iter);
+
+ for (const scalar_type &val : md1.real_variable("u"))
+ std::cout<<val<<std::endl;
+
+ std::cout<<std::endl;
+ for (const scalar_type &val : md2.real_variable("u"))
+ std::cout<<val<<std::endl;
+
+ std::cout << "Displacement dofs: " << mf.nb_dof() << std::endl;
+ std::cout << "Total dofs of model 1: " << md1.nb_dof() << std::endl;
+ std::cout << "Total dofs of model 2: " << md2.nb_dof() << std::endl;
+
+ return gmm::vect_dist2(md1.real_variable("u"), md2.real_variable("u")) <
1e-9 ? 0 : 1;
+}
diff --git a/tests/test_internal_variables.pl b/tests/test_internal_variables.pl
new file mode 100644
index 0000000..11d9526
--- /dev/null
+++ b/tests/test_internal_variables.pl
@@ -0,0 +1,45 @@
+# Copyright (C) 2019-2019 Yves Renard
+#
+# This file is a part of GetFEM++
+#
+# GetFEM++ is free software; you can redistribute it and/or modify it
+# under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation; either version 3 of the License, or
+# (at your option) any later version along with the GCC Runtime Library
+# Exception either version 3.1 or (at your option) any later version.
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+# License and GCC Runtime Library Exception for more details.
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program; if not, write to the Free Software Foundation,
+# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
+
+$er = 0;
+
+sub start_program
+{
+ my $def = $_[0];
+
+ # print ("def = $def\n");
+
+ open F, "./test_internal_variables $def 2>&1 |" or die;
+ while (<F>) {
+ if ($_ =~ /FAILED/) {
+ $er = 1;
+ print "============================================\n";
+ print $_, <F>;
+ }
+ print $_;
+ }
+ close(F); if ($?) { exit(1); }
+}
+
+start_program("-d DIFFICULTY=0");
+print "Test 1 DONE.\n";
+start_program("-d DIFFICULTY=1");
+print "Test 2 DONE.\n";
+
+if ($er == 1) { exit(1); }
+
+