[Getfem-commits] [getfem-commits] branch devel-logari81-internal-variables updated: Add unit test for internal variable condensation

[Konstantinos Poulios]

This is an automated email from the git hooks/post-receive script.

logari81 pushed a commit to branch devel-logari81-internal-variables
in repository getfem.

The following commit(s) were added to 
refs/heads/devel-logari81-internal-variables by this push:
     new 9d1c444  Add unit test for internal variable condensation
9d1c444 is described below

commit 9d1c4449796fe7b3c3adabeed57977c5d6ac7507
Author: Konstantinos Poulios <address@hidden>
AuthorDate: Sat Feb 1 11:28:19 2020 +0100

    Add unit test for internal variable condensation
---
 tests/Makefile.am          |   4 +
 tests/test_condensation.cc | 184 +++++++++++++++++++++++++++++++++++++++++++++
 tests/test_condensation.pl | 105 ++++++++++++++++++++++++++
 3 files changed, 293 insertions(+)

diff --git a/tests/Makefile.am b/tests/Makefile.am
index ba8b7a3..296934f 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -45,6 +45,7 @@ check_PROGRAMS =                   \
        test_assembly_assignment   \
        test_interpolated_fem      \
        test_internal_variables    \
+       test_condensation          \
        test_range_basis           \
        laplacian                  \
        laplacian_with_bricks      \
@@ -100,6 +101,7 @@ 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_condensation_SOURCES = test_condensation.cc
 test_tree_sorted_SOURCES = test_tree_sorted.cc
 test_mat_elem_SOURCES = test_mat_elem.cc
 test_slice_SOURCES = test_slice.cc
@@ -141,6 +143,7 @@ TESTS = \
        test_assembly_assignment.pl   \
        test_interpolated_fem.pl      \
        test_internal_variables.pl    \
+       test_condensation.pl          \
        test_range_basis.pl           \
        laplacian.pl                  \
        laplacian_with_bricks.pl      \
@@ -182,6 +185,7 @@ EXTRA_DIST =                                                
\
        test_int_set.pl                                         \
        test_interpolated_fem.pl                                \
        test_internal_variables.pl                              \
+       test_condensation.pl                                    \
        test_slice.pl                                           \
        test_mesh_im_level_set.pl                               \
        thermo_elasticity_electrical_coupling.pl                \
diff --git a/tests/test_condensation.cc b/tests/test_condensation.cc
new file mode 100644
index 0000000..2a3ec22
--- /dev/null
+++ b/tests/test_condensation.cc
@@ -0,0 +1,184 @@
+/*===========================================================================
+
+ Copyright (C) 2019-2020 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_export.h"
+#include "getfem/getfem_model_solvers.h"
+
+using bgeot::dim_type;
+using bgeot::size_type;
+using bgeot::scalar_type;
+using bgeot::base_node;
+
+static bool debug=false;
+
+int main(int argc, char *argv[]) {
+
+  gmm::set_traces_level(1);
+
+  bgeot::md_param PARAM;
+  PARAM.add_int_param("NX", 1);
+  PARAM.add_int_param("NY", 1);
+  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");
+
+  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(98) = getfem::select_faces_of_normal(m, outer_faces, base_node(-1, 
0), 0.001);
+  m.region(99) = getfem::select_faces_of_normal(m, outer_faces, base_node(1, 
0), 0.001);
+
+  const dim_type N(m.dim());
+  getfem::mesh_fem mf(m, N);
+  mf.set_classical_finite_element(FEM_ORDER);
+  if (DIFFICULTY % 1000 <= 99) {
+    dal::bit_vector kept;
+    kept.add(0, mf.nb_basic_dof());
+    kept.setminus(mf.basic_dof_on_region(98));
+    mf.reduce_to_basic_dof(kept);
+  }
+
+  getfem::mesh_im mim(m);
+  mim.set_integration_method(dim_type(2*IM_ORDER+1));
+
+  getfem::im_data mimd1(mim), mimd3(mim);
+  {
+    bgeot::multi_index vecsizes(1);
+    vecsizes[0] = 3;
+    mimd3.set_tensor_size(vecsizes);
+  }
+
+  getfem::model md1, md2;
+  md1.add_fem_variable("u", mf);
+  md2.add_fem_variable("u", mf);
+  if (DIFFICULTY % 100 > 19) {
+    md1.add_im_variable("eps", mimd3);
+    md2.add_internal_im_variable("eps", mimd3);
+  } else if (DIFFICULTY % 100 > 9) {
+    md1.add_im_variable("eps", mimd3);
+    md2.add_im_variable("eps", mimd3);
+  }
+  md1.add_im_variable("p", mimd1);
+  if (DIFFICULTY % 100 > 19 && DIFFICULTY % 100 <= 29)
+    md2.add_im_variable("p", mimd1);
+  else
+    md2.add_internal_im_variable("p", mimd1);
+
+
+  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, exprC;
+  if (DIFFICULTY % 100 > 9) {
+    exprA = 
"(-p*Id(2)+2*G*([[2/3,0;0,-1/3],[-1/3,0;0,2/3],[0,1;1,0]].eps)):Grad_Test_u";
+    exprB = "(eps-Grad_u:[[1,0;0,0],[0,0;0,1],[0,0.5;0.5,0]]).Test_eps";
+  } else
+    exprA = "(-p*Id(2)+2*G*(Sym(Grad_u)-Div_u*Id(2)/3)):Grad_Test_u";
+  if (DIFFICULTY % 10 == 0)
+    exprC = "(p+K*Div_u)*Test_p";
+  else
+    exprC = "((p+1e3*pow(p,3))+K*Div_u)*Test_p";
+  getfem::add_nonlinear_generic_assembly_brick(md1, mim, exprA);
+  getfem::add_nonlinear_generic_assembly_brick(md2, mim, exprA);
+  if (DIFFICULTY % 100 > 9) {
+    getfem::add_nonlinear_generic_assembly_brick(md1, mim, exprB);
+    getfem::add_nonlinear_generic_assembly_brick(md2, mim, exprB);
+  }
+  getfem::add_nonlinear_generic_assembly_brick(md1, mim, exprC);
+  getfem::add_nonlinear_generic_assembly_brick(md2, mim, exprC);
+
+  if (DIFFICULTY % 1000 > 99) {
+    md1.add_filtered_fem_variable("dirmult", mf, 98);
+    md2.add_filtered_fem_variable("dirmult", mf, 98);
+    getfem::add_linear_generic_assembly_brick(md1, mim, "u.dirmult", 98);
+    getfem::add_linear_generic_assembly_brick(md2, mim, "u.dirmult", 98);
+  }
+
+  // load
+  getfem::add_linear_generic_assembly_brick(md1, mim, "1e-3*Test_u(2)", 99);
+  getfem::add_linear_generic_assembly_brick(md2, mim, "1e-3*Test_u(2)", 99);
+
+  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;
+
+  getfem::default_newton_line_search ls1A, ls2A;
+  getfem::newton_search_with_step_control ls1B, ls2B;
+
+  std::cout<<"SOLVING MODEL 1 (without internal variables)"<<std::endl;
+  gmm::iteration iter(1E-9, 1, 30);
+  if (DIFFICULTY % 10000 > 999)
+    getfem::standard_solve(md1, iter, getfem::rselect_linear_solver(md1, 
"mumps"), ls1B);
+  else
+    getfem::standard_solve(md1, iter, getfem::rselect_linear_solver(md1, 
"mumps"), ls1A);
+
+  std::cout<<"SOLVING MODEL 2 (with internal variables)"<<std::endl;
+  iter.init();
+  if (DIFFICULTY % 10000 > 999)
+    getfem::standard_solve(md2, iter, getfem::rselect_linear_solver(md2, 
"mumps"), ls2B);
+  else
+    getfem::standard_solve(md2, iter, getfem::rselect_linear_solver(md2, 
"mumps"), ls2A);
+
+  if (debug) {
+    std::cout<<std::endl<<"u1:"<<std::endl;
+    for (const scalar_type &val : md1.real_variable("u")) 
std::cout<<val<<std::endl;
+    std::cout<<std::endl<<"u2:"<<std::endl;
+    for (const scalar_type &val : md2.real_variable("u")) 
std::cout<<val<<std::endl;
+
+    if (DIFFICULTY % 100 > 9) {
+      std::cout<<std::endl<<"eps1:"<<std::endl;
+      for (const scalar_type &val : md1.real_variable("eps")) 
std::cout<<val<<std::endl;
+      std::cout<<std::endl<<"eps2:"<<std::endl;
+      for (const scalar_type &val : md2.real_variable("eps")) 
std::cout<<val<<std::endl;
+    }
+
+    std::cout<<std::endl<<"p1:"<<std::endl;
+    for (const scalar_type &val : md1.real_variable("p")) 
std::cout<<val<<std::endl;
+    std::cout<<std::endl<<"p2:"<<std::endl;
+    for (const scalar_type &val : md2.real_variable("p"))
+    std::cout<<val<<std::endl;
+
+    getfem::vtk_export exp("test_internal_variable_condensation.vtk", true);
+    exp.exporting(mf);
+    exp.write_point_data(mf, md1.real_variable("u"), "displ1");
+    exp.write_point_data(mf, md2.real_variable("u"), "displ2");
+    exp.write_point_data(mf, gmm::sub_vector(md2.real_rhs(true), 
md2.interval_of_variable("u")), "u residual");
+  }
+
+  int ret=0;
+  ret += gmm::vect_dist1(md1.real_variable("u"), md2.real_variable("u")) < 
1e-9 ? 0 : 1;
+  if (DIFFICULTY % 100 > 9)
+    ret += gmm::vect_dist1(md1.real_variable("eps"), md2.real_variable("eps")) 
< 1e-9 ? 0 : 2;
+  ret += gmm::vect_dist1(md1.real_variable("p"), md2.real_variable("p")) < 
1e-9 ? 0 : 4;
+
+  std::cout<<"Test with difficulty "<<DIFFICULTY<<" returned "<<ret<<std::endl;
+  return ret;
+}
diff --git a/tests/test_condensation.pl b/tests/test_condensation.pl
new file mode 100644
index 0000000..b3da896
--- /dev/null
+++ b/tests/test_condensation.pl
@@ -0,0 +1,105 @@
+# Copyright (C) 2019-2020 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_condensation $def 2>&1 |" or die;
+  while (<F>) {
+    if ($_ =~ /FAILED/) {
+      $er = 1;
+      print "============================================\n";
+      print $_, <F>;
+    }
+    print $_;
+  }
+  close(F); if ($?) { exit(1); }
+}
+
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=0");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=10");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=11");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=20");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=21");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=30");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=31");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=100");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=101");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=110");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=111");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=120");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=121");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=130");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=131");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1000");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1001");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1010");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1011");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1020");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1021");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1030");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1031");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1100");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1101");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1110");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1111");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1120");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1121");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1130");
+print ".";
+start_program("-d NX=2 -d NX=3 -d DIFFICULTY=1131");
+print ".\n";
+
+if ($er == 1) { exit(1); }
+
+