Merge branch 'master' into 'change'

Conflicts:
	AUTHORS
	benchmarks/bessel/main.cpp
	benchmarks/kellogg/main.cpp
	benchmarks/layer-2/main.cpp
	benchmarks/layer/main.cpp
	benchmarks/line-singularity/main.cpp
	benchmarks/lshape/main.cpp
	benchmarks/neutronics-2-group-adapt/main.cpp
	benchmarks/neutronics-2-group-adapt/scaled_scalar_view.h
	benchmarks/screen/main.cpp
	benchmarks/smooth-aniso-x/main.cpp
	benchmarks/smooth-aniso-y/main.cpp
	benchmarks/smooth-iso/main.cpp
	examples/CMakeLists.txt
	examples/bracket/main.cpp
	examples/crack/main.cpp
	examples/iron-water/main.cpp
	examples/linear-advection-diffusion/main.cpp
	examples/neutronics-4-group-adapt/main.cpp
	examples/newton-np-timedep-adapt-system/main.cpp
	examples/ns-bearing/forms.cpp
	examples/ns-bearing/main.cpp
	examples/ns-two-phase/main.cpp
	examples/richards-seepage-adapt/main.cpp
	examples/richards-tracy-adapt/forms.cpp
	examples/richards-tracy-adapt/main.cpp
	examples/richards/forms.cpp
	examples/richards/main.cpp
	examples/saphir/main.cpp
	examples/singular-perturbation/main.cpp
	examples/thermoelasticity/main.cpp
	examples/waveguide/main.cpp
	examples/wire/main.cpp
	hermes_common/CMakeLists.txt
	hermes_common/_hermes_common.cpp
	hermes_common/_hermes_common_api_new.h
	hermes_common/common_time_period.h
	hermes_common/matrix.h
	hermes_common/python_api.cpp
	hermes_common/python_solvers.cpp
	python/hermes2d/_hermes2d.cpp
	src/feproblem.cpp
	src/feproblem.h
	src/h1_adapt.cpp
	src/h1_adapt.h
	src/linsystem.cpp
	src/linsystem.h
	src/matrix_old.h
	src/nonlinsystem.cpp
	src/nonlinsystem.h
	src/ref_selectors/selector.h
	src/refsystem.cpp
	src/solution.cpp
	src/space.cpp
	src/space.h
	src/weakform.cpp
	tests/benchmarks/smooth-aniso-x/main.cpp
	tests/benchmarks/smooth-aniso-y/main.cpp
	tests/benchmarks/smooth-iso/main.cpp
	tests/examples/bracket/main.cpp
	tests/examples/iron-water/main.cpp
	tests/tutorial/09-timedep/main.cpp
	tests/tutorial/11-system-adapt/main.cpp
	tutorial/03-poisson/main.cpp
	tutorial/10-adapt/main.cpp
	tutorial/11-system-adapt/main.cpp
	tutorial/12-general-adapt/main.cpp
	tutorial/15-newton-elliptic-adapt/main.cpp
	tutorial/33-exact-adapt/main.cpp
	tutorial/40-trilinos-linear/main.cpp
	tutorial/41-trilinos-nonlinear/main.cpp
	tutorial/44-trilinos-coupled/main.cpp

AUTHORS

@@ -19,4 +19,5 @@ Ma Zhonghua <mazhonghua83@gmail.com>
 Milan Hanus <mhanus@kma.zcu.cz>
 Damien Lebrun-Grandie <lebrungd@neo.tamu.edu>
 Lukas Korous <lukas.korous@gmail.com>
+Martin Lanzendorfer <lanz@cs.cas.cz>
 Valmor de Almeida <val.hpfem@gmail.com>

CMakeLists.txt

@@ -19,7 +19,7 @@ set(H2D_REAL     YES)
 set(H2D_COMPLEX  YES)
 
 # optional functionality
-set(WITH_PYTHON     NO)
+set(WITH_PYTHON     YES)
 set(WITH_EXAMPLES   YES)
 set(WITH_BENCHMARKS YES)
 set(WITH_TUTORIAL   YES)

benchmarks/neutronics-2-group-adapt/CMakeLists.txt

@@ -5,6 +5,6 @@ endif(NOT H2D_REAL)
 project(neutronics-2-group-adapt)
 set(HERMES ${HERMES_REAL_BIN})
 
-add_executable(${PROJECT_NAME} main.cpp)
+add_executable(${PROJECT_NAME} main.cpp scaled_scalar_view.cpp)
 include (../CMake.common)
 

benchmarks/neutronics-2-group-adapt/scaled_scalar_view.cpp

@@ -0,0 +1,28 @@
+#include "scaled_scalar_view.h"
+
+void ScaledScalarView::scale(double sc)
+{
+    /* FIXME: If I uncomment the lines below, I am getting:
+/home/certik1/repos/hermes2d/benchmarks/neutronics-2-group-adapt/scaled_scalar_view.cpp: In member function ‘void ScaledScalarView::scale(double)’:
+/home/certik1/repos/hermes2d/benchmarks/neutronics-2-group-adapt/scaled_scalar_view.cpp:4: error: ‘class ScaledScalarView’ has no member named ‘lin’
+/home/certik1/repos/hermes2d/benchmarks/neutronics-2-group-adapt/scaled_scalar_view.cpp:5: error: ‘class ScaledScalarView’ has no member named ‘mode3d’
+/home/certik1/repos/hermes2d/benchmarks/neutronics-2-group-adapt/scaled_scalar_view.cpp:6: error: ‘class ScaledScalarView’ has no member named ‘yscale’
+/home/certik1/repos/hermes2d/benchmarks/neutronics-2-group-adapt/scaled_scalar_view.cpp:7: error: ‘class ScaledScalarView’ has no member named ‘contours’
+/home/certik1/repos/hermes2d/benchmarks/neutronics-2-group-adapt/scaled_scalar_view.cpp:8: error: ‘class ScaledScalarView’ has no member named ‘cont_step’
+/home/certik1/repos/hermes2d/benchmarks/neutronics-2-group-adapt/scaled_scalar_view.cpp:9: error: ‘class ScaledScalarView’ has no member named ‘lin’
+/home/certik1/repos/hermes2d/benchmarks/neutronics-2-group-adapt/scaled_scalar_view.cpp:10: error: ‘class ScaledScalarView’ has no member named ‘refresh’
+
+
+However, the "lin" is declared in src/views/scalar_view.h as a protected member of the class ScalarView, so it should work.
+*/
+
+    /*
+    this->lin.lock_data();
+    if (this->mode3d)
+        this->yscale *= sc;
+    else if (this->contours)
+        this->cont_step *= sc;
+    this->lin.unlock_data();
+    this->refresh();
+    */
+}

doc/src/examples.rst

@@ -1625,3 +1625,30 @@ and hp-FEM are shown below.
 
 
 
+Richards Equation I
+-------------------
+
+**Git reference:** Example `richards 
+<http://git.hpfem.org/hermes2d.git/tree/HEAD:/examples/richards>`_.
+
+Description coming soon.
+
+
+Richards Equation II
+--------------------
+
+**Git reference:** Example `richards-tracy-adapt 
+<http://git.hpfem.org/hermes2d.git/tree/HEAD:/examples/richards-tracy-adapt>`_.
+
+Description coming soon.
+
+
+
+Richards Equation III
+---------------------
+
+**Git reference:** Example `richards-seepage-adapt 
+<http://git.hpfem.org/hermes2d.git/tree/HEAD:/examples/richards-seepage-adapt>`_.
+
+Description coming soon.
+

examples/CMake.common

@@ -25,7 +25,7 @@ include_directories(${hermes2d_SOURCE_DIR}/hermes_common/)
 include_directories(${PYTHON_INCLUDE_PATH} ${NUMPY_INCLUDE_PATH})
 include_directories(${TRILINOS_INCLUDE_DIR})
 target_link_libraries(${PROJECT_NAME} hermes_common ${HERMES}
-    ${PYTHON_LIBRARIES}
-    ${UMFPACK_LIBRARY} ${AMD_LIBRARY} ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES}
-    ${TRILINOS_LIBRARIES})
+        ${PYTHON_LIBRARIES}
+        ${UMFPACK_LIBRARY}
+        ${AMD_LIBRARY} ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${TRILINOS_LIBRARIES})
 add_dependencies(${PROJECT_NAME} ${HERMES})

examples/CMakeLists.txt

@@ -6,6 +6,7 @@ add_subdirectory(bracket)
 add_subdirectory(crack)
 add_subdirectory(crack-long)
 add_subdirectory(newton-np-timedep-adapt-system)
+add_subdirectory(newton-np-timedep-adapt-system/nonadaptive)
 add_subdirectory(iron-water)
 add_subdirectory(saphir)
 add_subdirectory(thermoelasticity)

examples/newton-np-timedep-adapt-system/CMakeLists.txt

@@ -1,4 +1,5 @@
 project(newton-np-timedep-adapt-system)
 
 add_executable(${PROJECT_NAME} main.cpp)
+
 include (../CMake.common)

examples/newton-np-timedep-adapt-system/circular.mesh

@@ -0,0 +1,46 @@
+r = 0.5e-3
+rr = 0.5*r
+vertices=
+{
+	{ 0, 0 },	#0
+	{ rr, 0 },	#1
+	{ r, 0 },	#2
+	{ 0, rr },	#3
+	{ 0, r },	#4
+	{ -rr, 0 },	#5
+	{ -r, 0 },	#6
+	{ 0, -rr },	#7
+	{ 0, -r }	#8
+}
+
+elements =
+{
+	{ 0, 1, 3, 0 },
+	{ 1, 2, 4, 3, 0},
+	{ 0, 3, 5, 0 },
+	{ 3, 4, 6, 5, 0},
+	{ 0, 5, 7, 0 },
+	{ 5, 6, 8, 7, 0},
+	{ 0, 7, 1, 0 },
+	{ 7, 8, 2, 1, 0 }
+}
+
+boundaries =
+{
+	{ 2, 4, 2 },
+	{ 4, 6, 1 },
+	{ 6, 8, 3 },
+	{ 8, 2, 1 }
+}
+
+curves =
+{
+	{ 2, 4, 90 },
+	{ 4, 6, 90 },
+	{ 6, 8, 90 },
+	{ 8, 2, 90 },
+	{ 1, 3, 90 },
+	{ 3, 5, 90 },
+	{ 5, 7, 90 },
+	{ 7, 1, 90 }
+}

examples/newton-np-timedep-adapt-system/nonadaptive/CMakeLists.txt

@@ -0,0 +1,3 @@
+project(newton-np-timedep-system)
+add_executable(${PROJECT_NAME} main.cpp)
+include (../../CMake.common)

examples/newton-np-timedep-adapt-system/nonadaptive/main.cpp

@@ -0,0 +1,225 @@
+#define H2D_REPORT_WARN
+#define H2D_REPORT_INFO
+#define H2D_REPORT_VERBOSE
+#define H2D_REPORT_FILE "application.log"
+
+#include "hermes2d.h"
+
+/** \addtogroup e_newton_np_timedep_adapt_system Newton Time-dependant System with Adaptivity
+ \{
+ \brief This example shows how to combine the automatic adaptivity with the Newton's method for a nonlinear time-dependent PDE system.
+
+ This example shows how to combine the automatic adaptivity with the
+ Newton's method for a nonlinear time-dependent PDE system.
+ The time discretization is done using implicit Euler or
+ Crank Nicholson method (see parameter TIME_DISCR).
+ The following PDE's are solved:
+ Nernst-Planck (describes the diffusion and migration of charged particles):
+ \f[dC/dt - D*div[grad(C)] - K*C*div[grad(\phi)]=0,\f]
+ where D and K are constants and C is the cation concentration variable,
+ phi is the voltage variable in the Poisson equation:
+ \f[ - div[grad(\phi)] = L*(C - C_0),\f]
+ where \f$C_0\f$, and L are constant (anion concentration). \f$C_0\f$ is constant
+ anion concentration in the domain and L is material parameter.
+ So, the equation variables are phi and C and the system describes the
+ migration/diffusion of charged particles due to applied voltage.
+ The simulation domain looks as follows:
+ \verbatim
+      2
+  +----------+
+  |          |
+ 1|          |1
+  |          |
+  +----------+
+      3
+ \endverbatim
+ For the Nernst-Planck equation, all the boundaries are natural i.e. Neumann.
+ Which basically means that the normal derivative is 0:
+ \f[ BC: -D*dC/dn - K*C*d\phi/dn = 0 \f]
+ For Poisson equation, boundary 1 has a natural boundary condition
+ (electric field derivative is 0).
+ The voltage is applied to the boundaries 2 and 3 (Dirichlet boundaries)
+ It is possible to adjust system paramter VOLT_BOUNDARY to apply
+ Neumann boundary condition to 2 (instead of Dirichlet). But by default:
+  - BC 2: \f$\phi = VOLTAGE\f$
+  - BC 3: \f$\phi = 0\f$
+  - BC 1: \f$\frac{d\phi}{dn} = 0\f$
+ */
+
+#define SIDE_MARKER 1
+#define TOP_MARKER 2
+#define BOT_MARKER 3
+#define NONCIRCULAR
+
+/*** Fundamental coefficients ***/
+const double D = 10e-11; 	            // [m^2/s] Diffusion coefficient
+const double R = 8.31; 		            // [J/mol*K] Gas constant
+const double T = 293; 		            // [K] Aboslute temperature
+const double F = 96485.3415;	        // [s * A / mol] Faraday constant
+const double eps = 2.5e-2; 	          // [F/m] Electric permeability
+const double mu = D / (R * T);        // Mobility of ions
+const double z = 1;		                // Charge number
+const double K = z * mu * F;          // Constant for equation
+const double L =  F / eps;	          // Constant for equation
+const double VOLTAGE = 1;	            // [V] Applied voltage
+const scalar C0 = 1200;	              // [mol/m^3] Anion and counterion concentration
+
+/* For Neumann boundary */
+const double height = 180e-6;	              // [m] thickness of the domain
+const double E_FIELD = VOLTAGE / height;    // Boundary condtion for positive voltage electrode
+
+
+/* Simulation parameters */
+const int PROJ_TYPE = 1;              // For the projection of the initial condition 
+                                      // on the initial mesh: 1 = H1 projection, 0 = L2 projection
+const double TAU = 0.05;              // Size of the time step
+const int T_FINAL = 5;                // Final time
+const int P_INIT = 3;       	        // Initial polynomial degree of all mesh elements.
+const int REF_INIT = 5;     	        // Number of initial refinements
+const bool MULTIMESH = false;	        // Multimesh?
+const int TIME_DISCR = 1;             // 1 for implicit Euler, 2 for Crank-Nicolson
+const int VOLT_BOUNDARY = 1;          // 1 for Dirichlet, 2 for Neumann
+
+/* Nonadaptive solution parameters */
+const double NEWTON_TOL = 1e-6;       // Stopping criterion for nonadaptive solution
+const int NEWTON_MAX_ITER = 20;       // Maximum allowed number of Newton iterations
+
+// Weak forms
+#include "../forms.cpp"
+
+/*** Boundary types and conditions ***/
+
+// Poisson takes Dirichlet and Neumann boundaries
+BCType phi_bc_types(int marker) {
+  return (marker == SIDE_MARKER)
+      ? BC_NATURAL : BC_ESSENTIAL;
+}
+
+// Nernst-Planck takes Neumann boundaries
+BCType C_bc_types(int marker) {
+  return BC_NATURAL;
+}
+
+// Diricleht Boundary conditions for Poisson equation.
+scalar C_essential_bc_values(int ess_bdy_marker, double x, double y) {
+  return 0;
+}
+
+// Diricleht Boundary conditions for Poisson equation.
+scalar phi_essential_bc_values(int ess_bdy_marker, double x, double y) {
+  return ess_bdy_marker == TOP_MARKER ? VOLTAGE : 0.0;
+}
+
+scalar voltage_ic(double x, double y, double &dx, double &dy) {
+  // y^2 function for the domain.
+  //return (y+100e-6) * (y+100e-6) / (40000e-12);
+  return 0.0;
+}
+
+scalar concentration_ic(double x, double y, double &dx, double &dy) {
+  return C0;
+}
+
+
+
+int main (int argc, char* argv[]) {
+  // load the mesh file
+  Mesh Cmesh, phimesh, basemesh;
+
+  H2DReader mloader;
+
+#ifdef CIRCULAR
+  mloader.load("../circular.mesh", &basemesh);
+  basemesh.refine_all_elements(0);
+  basemesh.refine_towards_boundary(TOP_MARKER, 2);
+  basemesh.refine_towards_boundary(BOT_MARKER, 4);
+  MeshView mview("Mesh", 0, 600, 400, 400);
+  mview.show(&basemesh);
+#else
+  mloader.load("../small.mesh", &basemesh);
+  basemesh.refine_all_elements(1);
+  //basemesh.refine_all_elements(1); // when only p-adapt is used
+  //basemesh.refine_all_elements(1); // when only p-adapt is used
+  basemesh.refine_towards_boundary(TOP_MARKER, REF_INIT);
+  //basemesh.refine_towards_boundary(BOT_MARKER, (REF_INIT - 1) + 8); // when only p-adapt is used
+  basemesh.refine_towards_boundary(BOT_MARKER, REF_INIT - 1);
+#endif
+
+  Cmesh.copy(&basemesh);
+  phimesh.copy(&basemesh);
+
+  // Spaces for concentration and the voltage
+  H1Space Cspace(&Cmesh, C_bc_types, C_essential_bc_values, P_INIT);
+  H1Space phispace(MULTIMESH ? &phimesh : &Cmesh, phi_bc_types, phi_essential_bc_values, P_INIT);
+
+  // The weak form for 2 equations
+  WeakForm wf(2);
+
+  Solution C_prev_time,    // prveious time step solution, for the time integration
+    C_prev_newton,   // solution convergin during the Newton's iteration
+    phi_prev_time,
+    phi_prev_newton;
+
+  // Add the bilinear and linear forms
+  // generally, the equation system is described:
+  if (TIME_DISCR == 1) {  // Implicit Euler.
+    wf.add_vector_form(0, callback(Fc_euler), H2D_ANY,
+		  Tuple<MeshFunction*>(&C_prev_time, &C_prev_newton, &phi_prev_newton));
+    wf.add_vector_form(1, callback(Fphi_euler), H2D_ANY, Tuple<MeshFunction*>(&C_prev_newton, &phi_prev_newton));
+    wf.add_matrix_form(0, 0, callback(J_euler_DFcDYc), H2D_UNSYM, H2D_ANY, &phi_prev_newton);
+    wf.add_matrix_form(0, 1, callback(J_euler_DFcDYphi), H2D_UNSYM, H2D_ANY, &C_prev_newton);
+    wf.add_matrix_form(1, 0, callback(J_euler_DFphiDYc), H2D_UNSYM);
+    wf.add_matrix_form(1, 1, callback(J_euler_DFphiDYphi), H2D_UNSYM);
+  } else {
+    wf.add_vector_form(0, callback(Fc_cranic), H2D_ANY, 
+		  Tuple<MeshFunction*>(&C_prev_time, &C_prev_newton, &phi_prev_newton, &phi_prev_time));
+    wf.add_vector_form(1, callback(Fphi_cranic), H2D_ANY, Tuple<MeshFunction*>(&C_prev_newton, &phi_prev_newton));
+    wf.add_matrix_form(0, 0, callback(J_cranic_DFcDYc), H2D_UNSYM, H2D_ANY, Tuple<MeshFunction*>(&phi_prev_newton, &phi_prev_time));
+    wf.add_matrix_form(0, 1, callback(J_cranic_DFcDYphi), H2D_UNSYM, H2D_ANY, Tuple<MeshFunction*>(&C_prev_newton, &C_prev_time));
+    wf.add_matrix_form(1, 0, callback(J_cranic_DFphiDYc), H2D_UNSYM);
+    wf.add_matrix_form(1, 1, callback(J_cranic_DFphiDYphi), H2D_UNSYM);
+  }
+
+  // Nonlinear solver
+  NonlinSystem nls(&wf, Tuple<Space*>(&Cspace, &phispace));
+
+  phi_prev_time.set_exact(MULTIMESH ? &phimesh : &Cmesh, voltage_ic);
+  C_prev_time.set_exact(&Cmesh, concentration_ic);
+
+  C_prev_newton.copy(&C_prev_time);
+  phi_prev_newton.copy(&phi_prev_time);
+
+  // Project the function init_cond() on the FE space
+  // to obtain initial coefficient vector for the Newton's method.
+  info("Projecting initial conditions to obtain initial vector for the Newton'w method.");
+  nls.project_global(Tuple<MeshFunction*>(&C_prev_time, &phi_prev_time), 
+      Tuple<Solution*>(&C_prev_newton, &phi_prev_newton));
+
+
+  //VectorView vview("electric field [V/m]", 0, 0, 600, 600);
+  ScalarView Cview("Concentration [mol/m3]", 0, 0, 800, 800);
+  ScalarView phiview("Voltage [V]", 650, 0, 600, 600);
+  phiview.show(&phi_prev_time);
+  Cview.show(&C_prev_time);
+  char title[100];
+
+  int nstep = (int) (T_FINAL / TAU + 0.5);
+  for (int n = 1; n <= nstep; n++) {
+    verbose("\n---- Time step %d ----", n);
+    bool verbose = true; // Default is false.
+    if (!nls.solve_newton(Tuple<Solution*>(&C_prev_newton, &phi_prev_newton),
+        NEWTON_TOL, NEWTON_MAX_ITER, verbose)) 
+          error("Newton's method did not converge.");
+    sprintf(title, "time step = %i", n);
+    phiview.set_title(title);
+    phiview.show(&phi_prev_newton);
+    Cview.set_title(title);
+    Cview.show(&C_prev_newton);
+    phi_prev_time.copy(&phi_prev_newton);
+    C_prev_time.copy(&C_prev_newton);
+  }
+  View::wait();
+
+  return 0;
+}
+/// \}