example_static_plate1c.m 4.83 KB
Newer Older
Andreas Zilian's avatar
Andreas Zilian committed
1
2
%% Example Plate Bending 2a
% example on how to use the integrated mesh generator for
3
% A two-dimensional quadratic plate structure composed of distorted elements
Andreas Zilian's avatar
Andreas Zilian committed
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
% -> extended patch test (using many distorted quadrilateral elements)
addpath([ '..' filesep '..' filesep '..' filesep ]);
%% Definition of sections
s1 = mafe.Section2D();
s1.E     = 2.1e6;  % [N/m^2]
s1.nue   = 0.3;    % [-]
s1.t     = 1.0e-3; % [m]
%  -----------------------------------------------------------MESH GENERATION---
%% Describe the goemetry of the problem for the mesh generator
% function that controls grid node distribution along geo lines
nf1 = mgen.MgNormFunction( mgen.MgNormFunctionType.Pow   , [0.0, 0.85] );
% create all geo points
gp1 = mgen.MgGeoPoint( [ 0.0,  0.0] );
gp2 = mgen.MgGeoPoint( [ 1.0,  0.0] );
gp3 = mgen.MgGeoPoint( [ 1.0,  1.0] );
gp4 = mgen.MgGeoPoint( [ 0.0,  1.0] );
% create all geo lines
n = 24;
gl1 = mgen.MgGeoLine2Point( [gp2, gp1],  n, nf1 ); % ask for 4 grid nodes along this line
gl2 = mgen.MgGeoLine2Point( [gp3, gp2],  n, nf1 ); % ask for 6 grid nodes along this line
gl3 = mgen.MgGeoLine2Point( [gp4, gp3],  n, nf1 ); % this will adapt automatically
gl4 = mgen.MgGeoLine2Point( [gp1, gp4],  n, nf1 ); % this will adapt automatically
% create all geo areas
ga1 = mgen.MgGeoArea4Line( [gl1, gl2, gl3, gl4] );
% create geo object and generate the mesh
geo = mgen.MgGeo( [gp1, gp2, gp3, gp4], [gl1, gl2, gl3, gl4], [ga1] );
% print geo and mesh data info
%fprintf(geo.print());
%  -----------------------------------------------------------DEFINE FE ITEMS---
%% Definition of nodes in the system
nodes = mafe.Node.empty(0,0);
for gn = geo.getAllGridNodes()
    nodes(end+1) = mafe.Node2D( gn.coord );
end
%% Definition of elements in the system
elems = mafe.Element.empty(0,0);
for ge = geo.getAllGridElems()
    elems(end+1) = mafe.Plate2D( nodes( [ ge.nodes.id ] ), s1 );
end
%% Definition of constraints in the system
%
testcase = 'bend';
%
% vector of constraints
const = mafe.Constraint.empty(0,0);
%
switch(testcase)
case 'bend'
    % impose constant m11, support in Disp3
    %
    % left edge: supported
    nn = nodes( [gl4.getGridNodes.id] ); % get nodes on GeoLine 4 (left)
    const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
    const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 );
    % right edge: distributed moment = 1.0 [Nm/m]
    nn = nodes( [gl2.getGridNodes.id] ); % get nodes on GeoLine 2 (right)
    const(end+1) = mafe.ConstraintEdge( nn, [mafe.DofType.Rota1, mafe.DofType.Rota2], mafe.ConstraintType.Neumann,  [0.0, 1.0] );
    %
case 'shear'
    % impose constant shear
    %
    % left edge: supported
    nn = nodes( [gl4.getGridNodes.id] ); % get nodes on GeoLine 4 (left)
    const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
    % all nodes: fix rotations
    nn = nodes( [ga1.getGridNodes.id] ); % get nodes on GeoArea 1
    const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 );
    const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota2, mafe.ConstraintType.Dirichlet, 0.0 );
    % right edge: vertical distributed force = 1.0 [N/m]
    nn = nodes( [gl2.getGridNodes.id] ); % get nodes on GeoLine 2 (right)
    const(end+1) = mafe.ConstraintEdge( nn, [mafe.DofType.Disp3, mafe.DofType.Disp3], mafe.ConstraintType.Neumann,  [0.0, 1.0] );
    %
case 'twist'
    % impose constant twist
    %
    % left edge: supported
    nn = nodes( [gl4.getGridNodes.id] ); % get nodes on GeoLine 4 (left)
    const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
    % lower edge: supported
    nn = nodes( [gl1.getGridNodes.id] ); % get nodes on GeoLine 1 (bottom)
    const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
    % top right corner: single force = 1.0 [N]
    nn = nodes( [gp3.getGridNodes.id] ); % get node at GeoPoint 3 (top right corner)
    const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Neumann,   1.0 );
    %
otherwise
    disp('Unknown test case! Constraints not set');
    const = [];
end
%  --------------------------------------------DEFINE FE PROBLEM AND ANALYSIS---
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
% %% Select an analysis type: static
ana = mafe.FeAnalysisStatic(fep);
% %% Calculate response
ana.analyse(); disp( sprintf('---> number dof =  %5d', fep.ndofs) );
%  -----------------------------------------------------------INSPECT RESULTS---
% %% Visualise system response
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
fep.plotSystem('reference');
fep.plotSystem('deformed', 1e-4);
%fep.plotSystem('shear13'); colorbar();
%fep.plotSystem('tensor', 1.e-2);
fep.printSystem();