[Getfem-commits] r4760 - in /trunk/getfem: doc/sphinx/source/userdoc/ interface/src/ interface/src/scilab/demos/ interface/tests/matlab/ src/...

[Yves . Renard]

Author: renard
Date: Fri Aug 29 21:48:50 2014
New Revision: 4760

URL: http://svn.gna.org/viewcvs/getfem?rev=4760&view=rev
Log:
interface for time integration scheme and adaptation of demos

Modified:
    trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst
    trunk/getfem/interface/src/gf_model_get.cc
    trunk/getfem/interface/src/gf_model_set.cc
    trunk/getfem/interface/src/scilab/demos/demo_wave_equation.sce
    trunk/getfem/interface/tests/matlab/demo_wave_equation.m
    trunk/getfem/src/getfem/getfem_models.h
    trunk/getfem/src/getfem_models.cc
    trunk/getfem/tests/heat_equation.cc
    trunk/getfem/tests/heat_equation.param

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=4760&r1=4759&r2=4760&view=diff
==============================================================================
--- trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst        (original)
+++ trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst        Fri Aug 29 
21:48:50 2014
@@ -11,7 +11,7 @@
 Compute arbitrary terms - high-level generic assembly procedures
 ================================================================
 
-This section presents the second version of generic assembly of |gf|. It is a 
high-level generic assembly in the sense that the language used to describe the 
assembly is quite close to the weak formulation of boundary value problems of 
partial differential equations. It mainly has been developed to circumvent the 
difficulties with the low-level generic assembly (see  :ref:`ud-gasm-low`) for 
which nonlinear terms are quite difficult to take into account. Conversely, a 
symbolic differentiation algorithm is used with this version to simplify the 
writing of new nonlinear terms. Moreover, the assembly language is compiled 
into optimized instructions before the evaluation on each integration point in 
order to obtain a rather optimal computational cost.
+This section presents what is now the main generic assembly of |gf|. It is a 
high-level generic assembly in the sense that the language used to describe the 
assembly is quite close to the weak formulation of boundary value problems of 
partial differential equations. It mainly has been developed to circumvent the 
difficulties with the previous low-level generic assembly (see  
:ref:`ud-gasm-low`) for which nonlinear terms were quite difficult to describe. 
Conversely, a symbolic differentiation algorithm is used with this version to 
simplify the writing of new nonlinear terms. Moreover, the assembly language is 
compiled into optimized instructions before the evaluation on each integration 
point in order to obtain a rather optimal computational cost.
 
 The header file to be included to use the high-level generic assembly 
procedures in C++ is :file:`getfem/generic\_assembly.h`.
 

Modified: trunk/getfem/interface/src/gf_model_get.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/gf_model_get.cc?rev=4760&r1=4759&r2=4760&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_model_get.cc  (original)
+++ trunk/getfem/interface/src/gf_model_get.cc  Fri Aug 29 21:48:50 2014
@@ -112,13 +112,27 @@
        out.pop().from_integer(int(md->model().nb_dof()));
        );
 
+    /address@hidden dt = ('get time step')
+      Gives the value of the time step. @*/
+    sub_command
+      ("get time step", 0, 0, 0, 1,
+       out.pop().from_scalar(md->model().get_time_step());
+       );
+
+    /address@hidden t = ('get time')
+      Give the value of the data `t` corresponding to the current time.
+      @*/
+    sub_command
+      ("get time", 0, 0, 0, 1,
+       out.pop().from_scalar(md->model().get_time());
+       );
+
     /address@hidden T = ('tangent_matrix')
       Return the tangent matrix stored in the model address@hidden/
     sub_command
       ("tangent_matrix", 0, 0, 0, 1,
        RETURN_SPARSE(real_tangent_matrix(), complex_tangent_matrix());
        );
-
 
     /address@hidden ('rhs')
       Return the right hand side of the tangent address@hidden/

Modified: trunk/getfem/interface/src/gf_model_set.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/gf_model_set.cc?rev=4760&r1=4759&r2=4760&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_model_set.cc  (original)
+++ trunk/getfem/interface/src/gf_model_set.cc  Fri Aug 29 21:48:50 2014
@@ -1937,6 +1937,89 @@
        );
 
 
+    /address@hidden ('shift variables for time integration')
+      Function used to shift the variables of a model to the data
+      corresponding of ther value on the previous time step for time
+      integration schemes. For each variable for which a time integration
+      scheme has been declared, the scheme is called to perform the shift.
+      This function has to be called between two time steps. @*/
+    sub_command
+      ("shift variables for time integration", 0, 0, 0, 0,
+       md->model().shift_variables_for_time_integration();
+       );
+
+    /address@hidden ('perform init time derivative', @scalar ddt)
+      By calling this function, indicates that the next solve will compute
+      the solution for a (very) small time step `ddt` in order to initalize
+      the data corresponding to the derivatives needed by time integration
+      schemes (mainly the initial time derivative for order one in time
+      problems  and the second order time derivative for second order in time
+      problems). The next solve will not change the value of the variables. @*/
+    sub_command
+      ("perform init time derivative", 1, 1, 0, 0,
+       double ddt = in.pop().to_scalar();
+       md->model().perform_init_time_derivative(ddt);
+       );
+
+    /address@hidden ('set time step', @scalar dt)
+      Set the value of the time step to `dt`. This value can be change
+      from a step to another for all one-step schemes (i.e for the moment
+      to all proposed time integration schemes). @*/
+    sub_command
+      ("set time step", 1, 1, 0, 0,
+       double dt = in.pop().to_scalar();
+       md->model().set_time_step(dt);
+       );
+
+    /address@hidden ('set time', @scalar t)
+      Set the value of the data `t` corresponding to the current time to `t`.
+      @*/
+    sub_command
+      ("set time", 1, 1, 0, 0,
+       double t = in.pop().to_scalar();
+       md->model().set_time(t);
+       );
+
+
+    /address@hidden ('add theta method for first order', @str varname, @scalar 
theta)
+      Attach a theta method for the time discretization of the variable
+      `varname`. Valid only if there is at most first order time derivative
+      of the variable. @*/
+    sub_command
+      ("add theta method for first order", 2, 2, 0, 0,
+       std::string varname = in.pop().to_string();
+       double theta = in.pop().to_scalar();
+       getfem::add_theta_method_for_first_order(md->model(), varname, theta);
+       );
+
+    /address@hidden ('add theta method for second order', @str varname, 
@scalar theta)
+      Attach a theta method for the time discretization of the variable
+      `varname`. Valid only if there is at most second order time derivative
+      of the variable. @*/
+    sub_command
+      ("add theta method for second order", 2, 2, 0, 0,
+       std::string varname = in.pop().to_string();
+       double theta = in.pop().to_scalar();
+       getfem::add_theta_method_for_second_order(md->model(), varname, theta);
+       );
+
+    /address@hidden ('add Newmark scheme', @str varname, @scalar beta, @scalar 
gamma)
+      Attach a theta method for the time discretization of the variable
+      `varname`. Valid only if there is at most second order time derivative
+      of the variable. @*/
+    sub_command
+      ("add Newmark scheme", 3, 3, 0, 0,
+       std::string varname = in.pop().to_string();
+       double beta = in.pop().to_scalar();
+       double gamma = in.pop().to_scalar();
+       getfem::add_Newmark_scheme(md->model(), varname, beta, gamma);
+       );
+
+
+
+
+    
+
     /address@hidden ind = ('add basic d on dt brick', @tmim mim, @str 
varnameU, @str dataname_dt[, @str dataname_rho[, @int region]])
     Add the standard discretization of a first order time derivative on
     `varnameU`. The parameter `dataname_rho` is the density which could

Modified: trunk/getfem/interface/src/scilab/demos/demo_wave_equation.sce
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/scilab/demos/demo_wave_equation.sce?rev=4760&r1=4759&r2=4760&view=diff
==============================================================================
--- trunk/getfem/interface/src/scilab/demos/demo_wave_equation.sce      
(original)
+++ trunk/getfem/interface/src/scilab/demos/demo_wave_equation.sce      Fri Aug 
29 21:48:50 2014
@@ -1,5 +1,5 @@
-// Simple demo of a wave equation solved with transient bricks
-// trace on;
+// Simple demo of a wave equation solved with the
+// Getfem tool for time integration schemes
 
 lines(0);
 stacksize('max');
@@ -38,6 +38,7 @@
 
 // interpolate the initial data
 U0 = gf_mesh_fem_get_eval(mf, list(list('y.*(y-1).*x.*(x-1).*x.*x')));
+V0 = 0*U0;
 
 md = gf_model('real');
 gf_model_set(md, 'add fem variable', 'u', mf, 2);
@@ -47,18 +48,22 @@
 // transient part.
 T     = 1.0; // For a good animation, choose 15 here (the computation is quite 
long then)
 dt    = 0.025;
-theta = 0.5;
-gf_model_set(md, 'add fem data', 'v', mf, 1, 2);
-gf_model_set(md, 'add initialized data', 'dt', [dt]);
-gf_model_set(md, 'add initialized data', 'theta', [theta]);
-gf_model_set(md, 'add basic d2 on dt2 brick', mim, 'u', 'v', 'dt', 'theta');
-gf_model_set(md, 'add theta method dispatcher', transient_bricks, 'theta');
+beta = 0.25;
+gamma = 0.5;
+
+gf_model_set(md, 'add Newmark scheme', 'u', beta, gamma);
+gf_model_set(md, 'add mass brick', mim, 'Dot2_u');
+gf_model_set(md, 'set time step', dt);
+
 
 // Initial data.
-gf_model_set(md, 'variable', 'u',  U0, 2);
-gf_model_set(md, 'first iter');
+gf_model_set(md, 'variable', 'Previous_u',  U0);
+gf_model_set(md, 'variable', 'Previous_Dot_u',  V0);
 
-gf_model_get(md, 'listbricks');
+
+// Initialisation of the acceleration 'Previous_Dot2_u'
+gf_model_set(md, 'perform init time derivative', dt/20.);
+gf_model_get(md, 'solve');
 
 // Iterations
 h = scf();
@@ -67,10 +72,10 @@
 Index = 0;
 
 for t=0:dt:T
+
   gf_model_get(md, 'solve');
-  gf_model_set(md, 'velocity update for order two theta method', 'u', 'v', 
'dt', 'theta');
   U = gf_model_get(md, 'variable', 'u');
-  V = gf_model_get(md, 'variable', 'v');
+  V = gf_model_get(md, 'variable', 'Dot_u');
 
   drawlater;
   clf();
@@ -89,7 +94,8 @@
 
   xs2png(h.figure_id, path + sprintf('/waveeq%03d.png',Index));
   Index = Index + 1;
-  gf_model_set(md, 'next iter');
+
+  gf_model_set(md, 'shift variables for time integration');
 end
 
 printf('demo wave_equation terminated\n');

Modified: trunk/getfem/interface/tests/matlab/demo_wave_equation.m
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_wave_equation.m?rev=4760&r1=4759&r2=4760&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_wave_equation.m    (original)
+++ trunk/getfem/interface/tests/matlab/demo_wave_equation.m    Fri Aug 29 
21:48:50 2014
@@ -16,8 +16,9 @@
 % Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
 
 
-% Simple demo of a wave equation solved with transient bricks
-% trace on;
+% Simple demo of a wave equation solved with a Newmark scheme with the
+% Getfem tool for time integration schemes
+
 gf_workspace('clear all');
 m = gf_mesh('cartesian',[0:.2:1],[0:.2:1]);
 % 
m=gf_mesh('import','structured','GT="GT_QK(2,1)";SIZES=[1,1];NOISED=1;NSUBDIV=[1,1];')
@@ -37,34 +38,37 @@
 
 % interpolate the initial data
 U0 = gf_mesh_fem_get(mf, 'eval', { 'y.*(y-1).*x.*(x-1).*x.*x' });
+V0 = 0*U0;
 
 md=gf_model('real');
 gf_model_set(md, 'add fem variable', 'u', mf, 2);
-transient_bricks = [gf_model_set(md, 'add Laplacian brick', mim, 'u')];
+gf_model_set(md, 'add Laplacian brick', mim, 'u');
 gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'u', mf, 1);
 
 % transient part.
 T = 15.0;
 dt = 0.025;
-theta = 0.5;
-gf_model_set(md, 'add fem data', 'v', mf, 1, 2);
-gf_model_set(md, 'add initialized data', 'dt', [dt]);
-gf_model_set(md, 'add initialized data', 'theta', [theta]);
-gf_model_set(md, 'add basic d2 on dt2 brick', mim, 'u', 'v', 'dt', 'theta');
-gf_model_set(md, 'add theta method dispatcher', transient_bricks, 'theta');
+beta = 0.25;
+gamma = 0.5;
+
+gf_model_set(md, 'add Newmark scheme', 'u', beta, gamma);
+gf_model_set(md, 'add mass brick', mim, 'Dot2_u');
+gf_model_set(md, 'set time step', dt);
 
 % Initial data.
-gf_model_set(md, 'variable', 'u',  U0, 2);
-gf_model_set(md, 'first iter');
+gf_model_set(md, 'variable', 'Previous_u',  U0);
+gf_model_set(md, 'variable', 'Previous_Dot_u',  V0);
 
-gf_model_get(md, 'listbricks');
+% Initialisation of the acceleration 'Previous_Dot2_u'
+gf_model_set(md, 'perform init time derivative', dt/20.);
+gf_model_get(md, 'solve');
 
 % Iterations
 for t = 0:dt:T
+
   gf_model_get(md, 'solve');
-  gf_model_set(md, 'velocity update for order two theta method', 'u', 'v', 
'dt', 'theta');
   U = gf_model_get(md, 'variable', 'u');
-  V = gf_model_get(md, 'variable', 'v');
+  V = gf_model_get(md, 'variable', 'Dot_u');
 
   subplot(2,1,1); gf_plot(mf, U, 'mesh', 'on', 'contour', .01:.01:.1); 
   colorbar; title(sprintf('computed solution u for t=%g', t));
@@ -73,11 +77,7 @@
   colorbar; title(sprintf('computed solution du/dt for t=%g', t));
   pause(0.1);
 
-  gf_model_set(md, 'next iter');
+  gf_model_set(md, 'shift variables for time integration');
 end;
 
 
-
-
-
-

Modified: trunk/getfem/src/getfem/getfem_models.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_models.h?rev=4760&r1=4759&r2=4760&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_models.h     (original)
+++ trunk/getfem/src/getfem/getfem_models.h     Fri Aug 29 21:48:50 2014
@@ -1058,12 +1058,6 @@
         return real_variable(varname)[0];
     }
 
-
-    /** Initialisation of the time integration. */
-    void init_time_integration(scalar_type t = scalar_type(0)) {
-      /* init_step = false; */ time_integration = true; set_time(t);
-    }
-
     void set_time_step(scalar_type dt) { time_step = dt; }
     scalar_type get_time_step(void) const { return time_step; }
     scalar_type get_init_time_step(void) const { return init_time_step; }
@@ -1149,7 +1143,7 @@
       init(); complex_version = comp_version;
       is_linear_ = is_symmetric_ = is_coercive_ = true;
       leading_dim = 0;
-      time_integration = false; init_step = false; time_step = scalar_type(1);
+      time_integration = 0; init_step = false; time_step = scalar_type(1);
     }
 
     /** check consistency of RHS and Stiffness matrix for brick with

Modified: trunk/getfem/src/getfem_models.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_models.cc?rev=4760&r1=4759&r2=4760&view=diff
==============================================================================
--- trunk/getfem/src/getfem_models.cc   (original)
+++ trunk/getfem/src/getfem_models.cc   Fri Aug 29 21:48:50 2014
@@ -763,6 +763,7 @@
 //                    "the precomputation of time derivative will not be "
 //                    "executed. Caution: You have to care by yourself of "
 //                    "the compatbility of the operation");
+    time_integration = 1;
   }
 
   void model::copy_init_time_derivative(void) {

Modified: trunk/getfem/tests/heat_equation.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/tests/heat_equation.cc?rev=4760&r1=4759&r2=4760&view=diff
==============================================================================
--- trunk/getfem/tests/heat_equation.cc (original)
+++ trunk/getfem/tests/heat_equation.cc Fri Aug 29 21:48:50 2014
@@ -247,7 +247,7 @@
   
   gmm::iteration iter(residual, 0, 40000);
   
-  model.init_time_integration(0.);
+  model.set_time(0.);
   model.set_time_step(dt);
   
   // Null initial value for the temperature.

Modified: trunk/getfem/tests/heat_equation.param
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/tests/heat_equation.param?rev=4760&r1=4759&r2=4760&view=diff
==============================================================================
--- trunk/getfem/tests/heat_equation.param      (original)
+++ trunk/getfem/tests/heat_equation.param      Fri Aug 29 21:48:50 2014
@@ -24,17 +24,12 @@
 if (N == 1)
   MESH_TYPE = 'GT_PK(1,1)';          % segments
   FEM_TYPE = 'FEM_PK(1,2)';          % P1 for segments
-%  FEM_TYPE = 'FEM_HERMITE(1)'; DATA_FEM_TYPE = 'FEM_PK(1,3)';
   INTEGRATION = 'IM_GAUSS1D(6)';     % Gauss-Legendre integration on segments
 end
 
 if (N == 2 && ~QUAD)
   MESH_TYPE = 'GT_PK(2,1)';          % triangles
   FEM_TYPE = 'FEM_PK(2,1)';          % P1 for triangles
-  %FEM_TYPE = 'FEM_REDUCED_HCT_TRIANGLE'; DATA_FEM_TYPE = 'FEM_PK(2,3)';
-  %INTEGRATION = 'IM_HCT_COMPOSITE(IM_TRIANGLE(6))';
-  %FEM_TYPE = 'FEM_HERMITE(2)'; DATA_FEM_TYPE = 'FEM_PK(2,2)';
-  %FEM_TYPE = 'FEM_ARGYRIS'; DATA_FEM_TYPE = 'FEM_PK(2,5)';
   INTEGRATION = 'IM_TRIANGLE(6)';    % quadrature rule for polynomials up
                                      % to degree 6 on triangles
 end