[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] r4760 - in /trunk/getfem: doc/sphinx/source/userdoc/ in
From: |
Yves . Renard |
Subject: |
[Getfem-commits] r4760 - in /trunk/getfem: doc/sphinx/source/userdoc/ interface/src/ interface/src/scilab/demos/ interface/tests/matlab/ src/... |
Date: |
Fri, 29 Aug 2014 19:48:50 -0000 |
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
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] r4760 - in /trunk/getfem: doc/sphinx/source/userdoc/ interface/src/ interface/src/scilab/demos/ interface/tests/matlab/ src/...,
Yves . Renard <=