%--------------------------------------------------------------------------
% PURPOSE
% Test multifiber beams - incremental solving - elastoplastic and damage
% laws - bending loading
%--------------------------------------------------------------------------
% REFERENCES
% Chaymaa LEJOUAD
% 15-02-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
%% Global variables
global ME TP
%% Definition of the section
FILE = '94_1.mail';
ME = INPUT.ACQU(FILE,'MAIL');
%% Definition of the model - cross section scale
MO1 = MODEL('SBET' ,'MECHANICS','ELASTICITY','ISOTROPIC','DAMAGE','MAZARS','QUAS');
MO2 = MODEL('SAINF','MECHANICS','ELASTICITY','ISOTROPIC','PLASTICITY','VONMISES','POJS');
MO3 = MODEL('SASUP','MECHANICS','ELASTICITY','ISOTROPIC','PLASTICITY','VONMISES','POJS');
MOT = [MO1 [MO2 MO3]];
%% Definition of the materials - cross section scale
MA1 = CHAMELEM.MATE(MO1,'youn',36000e6,'nu',0.2,'rho',...
2500,'ay',0.83,'az',0.83,'at',0.8,'k0',1.0e-4,'bt',12000,...
'beta',1.06,'ac',1.2,'bc',1500);
MA2 = CHAMELEM.MATE(MO2,'youn',200000e6,'nu',0.3,'rho',...
7100,'sect',pi*0.012^2/4,'sigy',500e6,'K',0.01*210000e6);
MA3 = CHAMELEM.MATE(MO3,'youn',200000e6,'nu',0.3,'rho',...
7100,'sect',pi*0.012^2/4,'sigy',500e6,'K',0.01*210000e6);
MAT = [MA1 [MA2 MA3]];
%% Definition of the beam
FILE = '94_2.mail';
ME = INPUT.ACQU(FILE,'MAIL');
%% Definition of the model - beam scale
MOB1 = MODEL('LPOUTRE','MECHANICS','ELASTICITY','SECTION','PLASTICITY',...
'SECTION','TIMO');
MOBT = MOB1;
%% Topology
TP = TOPOLOGY(MOBT);
%% Definition of the material - beam scale
MAB1 = CHAMELEM.MATE(MOB1,'mods',MOT,'mats',MAT);
MABT = MAB1;
%% Dirichlet boundary conditions
% Flexion 4 points
% CL1 = MATRICE('DIRI','P1G',1,2);
% CL2 = MATRICE('DIRI','P1D',2);
% CL3= MATRICE('DIRI','PVERIN1',2);
% CL4= MATRICE('DIRI','PVERIN2',2);
% CL5= MATRICE('DIRI','LPOUTRE',3,4,5);
%
% CLT = [CL1 [CL2 [CL3 [CL4 CL5]]]];
% Flexion 3 points
CL1 = MATRICE('DIRI','P1G',1,2);
CL2 = MATRICE('DIRI','P1D',2);
CL3= MATRICE('DIRI','P0',2);
CL4= MATRICE('DIRI','LPOUTRE',3,4,5);
CLT = [CL1 [CL2 [CL3 CL4]]];
%% Loadings
%Flexion 4 points
% FO1 = CHPOINT.DEPI(CL3,0.1);
% EV1 = EVOL([0 1],[0 1],'Time','Displacement (m)');
% CHT1 = TIMELOAD(FO1,EV1,'DIRI');
%
% FO2 = CHPOINT.DEPI(CL4,0.1);
% EV2 = EVOL([0 1],[0 1],'Time','Displacement (m)');
% CHT2 = TIMELOAD(FO2,EV2,'DIRI');
%
% CHT = [CHT1 CHT2];
%Flexion 3 points
FO1 = CHPOINT.DEPI(CL3,0.1);
EV1 = EVOL([0 1],[0 1],'Time','Displacement (m)');
CHT = TIMELOAD(FO1,EV1,'DIRI');
%% Static analysis
PB1 = PROBLEM('model',MOBT,'mater',MABT,'diric',CLT,'loadt',CHT,'comp_time',0:0.01:1,...
'solve_type','SECANT_NEWTON','stif_update',10);
SOL = SOLVE(PB1);
%% Reaction curve
% Results from CastLab
EV_OUT = EVOL.REAC(SOL,CL2,EV1,2,'P0');
EV_OUT.absc = EV_OUT.absc;
EV_OUT.ordo = -EV_OUT.ordo;
plot(EV_OUT);
% Reference results (Cast3M)
ref = ...
[
0.0000 0.0000
1.00000E-03 3421.0
2.00000E-03 6842.1
3.00000E-03 8740.1
4.00000E-03 9032.3
5.00000E-03 8430.8
6.00000E-03 7731.4
7.00000E-03 7447.0
8.00000E-03 7285.4
9.00000E-03 7298.0
1.00000E-02 7333.6
1.10000E-02 7407.9
1.20000E-02 7570.4
1.30000E-02 7799.8
1.40000E-02 8073.2
1.50000E-02 8381.5
1.60000E-02 8693.7
1.70000E-02 9032.5
1.80000E-02 9380.8
1.90000E-02 9738.1
2.00000E-02 10124.
2.10000E-02 10477.
2.20000E-02 10809.
2.30000E-02 11129.
2.40000E-02 11444.
2.50000E-02 11761.
2.60000E-02 12085.
2.70000E-02 12417.
2.80000E-02 12757.
2.90000E-02 13102.
3.00000E-02 13450.
3.10000E-02 13797.
3.20000E-02 14142.
3.30000E-02 14424.
3.40000E-02 14386.
3.50000E-02 14359.
3.60000E-02 14365.
3.70000E-02 14390.
3.80000E-02 14426.
3.90000E-02 14468.
4.00000E-02 14508.
4.10000E-02 14544.
4.20000E-02 14574.
4.30000E-02 14599.
4.40000E-02 14623.
4.50000E-02 14648.
4.60000E-02 14675.
4.70000E-02 14705.
4.80000E-02 14737.
4.90000E-02 14769.
5.00000E-02 14802.
5.10000E-02 14835.
5.20000E-02 14867.
5.30000E-02 14898.
5.40000E-02 14928.
5.50000E-02 14958.
5.60000E-02 14989.
5.70000E-02 15019.
5.80000E-02 15050.
5.90000E-02 15081.
6.00000E-02 15112.
6.10000E-02 15163.
6.20000E-02 15191.
6.30000E-02 15213.
6.40000E-02 15158.
6.50000E-02 15099.
6.60000E-02 15130.
6.70000E-02 15161.
6.80000E-02 15186.
6.90000E-02 15201.
7.00000E-02 15205.
7.10000E-02 15198.
7.20000E-02 15182.
7.30000E-02 15163.
7.40000E-02 15145.
7.50000E-02 15133.
7.60000E-02 15129.
7.70000E-02 15141.
7.80000E-02 15160.
7.90000E-02 15184.
8.00000E-02 15210.
8.10000E-02 15226.
8.20000E-02 15241.
8.30000E-02 15256.
8.40000E-02 15271.
8.50000E-02 15285.
8.60000E-02 15298.
8.70000E-02 15310.
8.80000E-02 15323.
8.90000E-02 15335.
9.00000E-02 15347.
9.10000E-02 15359.
9.20000E-02 15371.
9.30000E-02 15385.
9.40000E-02 15398.
9.50000E-02 15413.
9.60000E-02 15428.
9.70000E-02 15443.
9.80000E-02 15458.
9.90000E-02 15473.
0.10000 15488.];
hold on
plot(ref(:,1),ref(:,2),'go')
legend('CastLab','Cast3M','Location','Best')
close all
%% Non regression test
if abs(EV_OUT.ordo(end) - ref(end,2))/ref(end,2) > 5.0e-3
error('TEST IS NOT SUCCESSFUL')
else
disp('---------------------------------')
disp('094_INCR_DAM_PLAS_F3P_MULTIFIBER_TIMO')
disp('TEST IS SUCCESSFUL')
disp('---------------------------------')
end