%--------------------------------------------------------------------------
% PURPOSE
% Test restart option - 2D test - Linear elastic - test of the restart
% option with non null initial values - two different solvers used
%--------------------------------------------------------------------------
% REFERENCES
% Benjamin RICHARD
% 16-04-2018
%--------------------------------------------------------------------------
% COMMENTS
%
%
%
%--------------------------------------------------------------------------
% MIT License
%
% Copyright (c) 2018 Benjamin Richard
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the "Software
% "), to deal in the Software without restriction, including without
% limitation the rights to use, copy, modify, merge, publish, distribute,
% sublicense, and/or sell copies of the Software, and to permit persons
% to whom the Software is furnished to do so, subject to the following
% conditions:
%
% The above copyright notice and this permission notice shall be included
% in all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
% OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
% MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
% IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
% CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
% TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
% SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%--------------------------------------------------------------------------
%% Clearing off
fclose all;
clear
close all
%% Declaration de variables global
global options ME TP;
%% Definition des options
options.mode = 'PLANE_STRESS';
%% Loading of the input datafile
FILE = '105.mail';
ME = INPUT.ACQU(FILE,'MAIL');
%% Definition of the model
MO1 = MODEL('ST','MECHANICS','ELASTICITY','ISOTROPIC');
MOT = MO1;
%% Topology
TP = TOPOLOGY(MOT);
%% Definition of the material
MA1 = CHAMELEM.MATE(MO1,'youn',36000e6,'nu',0.2);
MAT = MA1;
%% Static analysis - STEP 1
% Set the BCs
CL1 = MATRICE('DIRI','P1',1,2);
CL2 = MATRICE('DIRI','P2',2);
CL3 = MATRICE('DIRI','P4',1);
CL4 = MATRICE('DIRI','L3',2);
CLT = [CL1 [CL2 [CL3 CL4]]];
% Loading
FO1 = CHPOINT.DEPI(CL4,1);
EV1 = EVOL([0 1],[0 5e-3],'Time','Displacement (m)');
CHT = TIMELOAD(FO1,EV1,'DIRI');
LI1 = 0:0.01:1;
% First problem to solve
PB1 = PROBLEM('model',MOT,'mater',MAT,'diric',CLT,'loadt',CHT,'comp_time',LI1,...
'solve_type','QUASI_NEWTON');
SOL1 = SOLVE(PB1);
EV_1 = EVOL.REAC(SOL1,CL4,EV1,2,'P3',2);
%% Static analysis - STEP 2
% Clear off the curent topology
TOPOLOGY.ERASE_DIRICHLET;
% Set up the new BCs
CL1bis = MATRICE('DIRI','P1',1,2);
CL2bis = MATRICE('DIRI','P2',2);
CL3bis = MATRICE('DIRI','P4',1);
CL4bis = MATRICE('DIRI','L3',2);
CL5bis = MATRICE('DIRI','L2',1);
CLTbis = [CL1bis [CL2bis [CL3bis [CL4bis CL5bis]]]];
% Loading
FO1bis = CHPOINT.DEPI(CL4bis,1);
EV1bis = EVOL([0 1],[5e-3 5e-3],'Time','Displacement (m)');
CH1bis = TIMELOAD(FO1bis,EV1bis,'DIRI');
FO2bis = CHPOINT.DEPI(CL5bis,1);
EV2bis = EVOL([0 1],[-1.0e-3 5e-3],'Time','Displacement (m)');
CH2bis = TIMELOAD(FO2bis,EV2bis,'DIRI');
% Second problem to solve
PB2 = PROBLEM('model',MOT,'mater',MAT,'diric',CLTbis,'loadt',[CH1bis CH2bis],'comp_time',LI1,...
'solve_type','SECANT_NEWTON',...
'restart',true,'cham',SOL1.cham_element(end).cham,'displacement',SOL1.displacement(end));
SOL2 = SOLVE(PB2);
for i=1:length(LI1)-1
eps11(i) = SOL1.cham_element(i).cham.cham_val.val_elem(1).epsf{1}(1);
eps22(i) = SOL1.cham_element(i).cham.cham_val.val_elem(1).epsf{1}(2);
eps12(i) = SOL1.cham_element(i).cham.cham_val.val_elem(1).epsf{1}(4);
end
for i=1:length(LI1)-1
eps11(i+length(LI1)-1) = SOL2.cham_element(i).cham.cham_val.val_elem(1).epsf{1}(1);
eps22(i+length(LI1)-1) = SOL2.cham_element(i).cham.cham_val.val_elem(1).epsf{1}(2);
eps12(i+length(LI1)-1) = SOL2.cham_element(i).cham.cham_val.val_elem(1).epsf{1}(4);
end
figure
box on
grid on
hold on
plot([0 eps11],'r')
plot([0 eps22],'c')
plot([0 eps12],'g')
xlabel('Pseudo time')
ylabel('Strain')
legend('\epsilon 11','\epsilon 22','\epsilon 12','Location','Best')
close all
%% Non regression test
if abs(eps12(end) + 8.9778e-20)/8.9778e-20 > 1.0e-4 || ...
(abs(eps11(end)) - 0.0050)/0.0050 > 1.0e-4 || ...
(abs(eps22(end)) - 0.0050)/0.0050 > 1.0e-4
error('TEST IS NOT SUCCESSFUL')
else
disp('---------------------------------')
disp('105_RESTART_TEST_NON_NULL_INITIA_VALUES')
disp('TEST IS SUCCESSFUL')
disp('---------------------------------')
end