Added tutorial 2 in statics on plate bending, step-by-step solution.

tutorials/statics/tutorial2/tutorial2_step1.m

+%% Example Plate bending
+% example on how to use the integrated mesh generator for
+% a specifically supported L-shaped plate domain under area and point loads
+%
+% step 1 (task 1):
+% 
+% generate the geometrical description and finite element mesh
+%
+addpath([ '..' filesep '..' filesep '..' filesep ]);
+%  -----------------------------------------------------------MESH GENERATION---
+%% Describe the goemetry of the problem for the mesh generator
+% helfpul parameters (problem specific)
+d = 2; n = 1; % d: length measure; n = number of grid nodes per d
+% function that controls grid node distribution along geo lines
+nf1 = mgen.MgNormFunction( mgen.MgNormFunctionType.Pow   , [0.0, 1.0] );
+nf2 = mgen.MgNormFunction( mgen.MgNormFunctionType.PowRev, [0.0, 1.0] );
+% create all geo points %% geo point coordinates specified in brackets []
+gp1 = mgen.MgGeoPoint( [0*d, 0*d] );
+gp2 = mgen.MgGeoPoint( [4*d, 0*d] );
+gp3 = mgen.MgGeoPoint( [0*d, 3*d] );
+gp4 = mgen.MgGeoPoint( [4*d, 3*d] );
+gp5 = mgen.MgGeoPoint( [8*d, 3*d] );
+gp6 = mgen.MgGeoPoint( [0*d, 8*d] );
+gp7 = mgen.MgGeoPoint( [4*d, 8*d] );
+gp8 = mgen.MgGeoPoint( [8*d, 8*d] );
+% create all geo lines %% geo line identified as continuous line between specified geo points, n presenting number of nodes per element
+gl1  = mgen.MgGeoLine2Point( [gp1, gp2], 4*n, nf2 );
+gl2  = mgen.MgGeoLine2Point( [gp3, gp4],   2, nf2 );
+gl3  = mgen.MgGeoLine2Point( [gp6, gp7],   2, nf2 );
+gl4  = mgen.MgGeoLine2Point( [gp4, gp5], 4*n, nf1 );
+gl5  = mgen.MgGeoLine2Point( [gp7, gp8],   2, nf1 );
+gl6  = mgen.MgGeoLine2Point( [gp1, gp3], 3*n, nf2 );
+gl7  = mgen.MgGeoLine2Point( [gp2, gp4],   2, nf2 );
+gl8  = mgen.MgGeoLine2Point( [gp3, gp6], 5*n, nf1 );
+gl9  = mgen.MgGeoLine2Point( [gp4, gp7],   2, nf1 );
+gl10 = mgen.MgGeoLine2Point( [gp8, gp5],   2, nf2 );
+% create all geo areas %% geo area identified as area enclosed by specified geo lines
+ga1 = mgen.MgGeoArea4Line( [gl1, gl7, gl2, gl6] );
+ga2 = mgen.MgGeoArea4Line( [gl2, gl9, gl3, gl8] );
+ga3 = mgen.MgGeoArea4Line( [gl4, gl10, gl5, gl9] );
+% create geo object and generate the mesh %% include all geo points, geo lines and geo areas
+geo = mgen.MgGeo( [gp1, gp2, gp3, gp4, gp5, gp6, gp7, gp8], ...
+                  [gl1, gl2, gl3, gl4, gl5, gl6, gl7, gl8, gl9, gl10], ...
+                  [ga1, ga2, ga3] );
+% print geo and mesh data info
+% fprintf(geo.print());
+% plot mesh
+cla; clf;
+fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on; %% this is used to properly position the figure
+geo.plot();

tutorials/statics/tutorial2/tutorial2_step2.m

+%% Example Plate bending
+% example on how to use the integrated mesh generator for
+% a specifically supported L-shaped plate domain under area and point loads
+%
+% step 2 (task 2.1-2.4):
+% 
+% define section data, nodes and elements, 
+% kinematic and static boundary conditions (constraints)
+%
+addpath([ '..' filesep '..' filesep '..' filesep ]);
+%  -----------------------------------------------------------MESH GENERATION---
+%% Describe the goemetry of the problem for the mesh generator
+% helfpul parameters (problem specific)
+d = 2; n = 1; % d: length measure; n = number of grid nodes per d
+% function that controls grid node distribution along geo lines
+nf1 = mgen.MgNormFunction( mgen.MgNormFunctionType.Pow   , [0.0, 1.0] );
+nf2 = mgen.MgNormFunction( mgen.MgNormFunctionType.PowRev, [0.0, 1.0] );
+% create all geo points %% geo point coordinates specified in brackets []
+gp1 = mgen.MgGeoPoint( [0*d, 0*d] );
+gp2 = mgen.MgGeoPoint( [4*d, 0*d] );
+gp3 = mgen.MgGeoPoint( [0*d, 3*d] );
+gp4 = mgen.MgGeoPoint( [4*d, 3*d] );
+gp5 = mgen.MgGeoPoint( [8*d, 3*d] );
+gp6 = mgen.MgGeoPoint( [0*d, 8*d] );
+gp7 = mgen.MgGeoPoint( [4*d, 8*d] );
+gp8 = mgen.MgGeoPoint( [8*d, 8*d] );
+% create all geo lines %% geo line identified as continuous line between specified geo points, n presenting number of nodes per element
+gl1  = mgen.MgGeoLine2Point( [gp1, gp2], 4*n, nf2 );
+gl2  = mgen.MgGeoLine2Point( [gp3, gp4],   2, nf2 );
+gl3  = mgen.MgGeoLine2Point( [gp6, gp7],   2, nf2 );
+gl4  = mgen.MgGeoLine2Point( [gp4, gp5], 4*n, nf1 );
+gl5  = mgen.MgGeoLine2Point( [gp7, gp8],   2, nf1 );
+gl6  = mgen.MgGeoLine2Point( [gp1, gp3], 3*n, nf2 );
+gl7  = mgen.MgGeoLine2Point( [gp2, gp4],   2, nf2 );
+gl8  = mgen.MgGeoLine2Point( [gp3, gp6], 5*n, nf1 );
+gl9  = mgen.MgGeoLine2Point( [gp4, gp7],   2, nf1 );
+gl10 = mgen.MgGeoLine2Point( [gp8, gp5],   2, nf2 );
+% create all geo areas %% geo area identified as area enclosed by specified geo lines
+ga1 = mgen.MgGeoArea4Line( [gl1, gl7, gl2, gl6] );
+ga2 = mgen.MgGeoArea4Line( [gl2, gl9, gl3, gl8] );
+ga3 = mgen.MgGeoArea4Line( [gl4, gl10, gl5, gl9] );
+% create geo object and generate the mesh %% include all geo points, geo lines and geo areas
+geo = mgen.MgGeo( [gp1, gp2, gp3, gp4, gp5, gp6, gp7, gp8], ...
+                  [gl1, gl2, gl3, gl4, gl5, gl6, gl7, gl8, gl9, gl10], ...
+                  [ga1, ga2, ga3] );
+% print geo and mesh data info
+% fprintf(geo.print());
+% plot mesh
+cla; clf;
+fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on; %% this is used to properly position the figure
+geo.plot();
+%  -----------------------------------------------------------DEFINE FE ITEMS---
+%% Definition of sections
+s1 = mafe.Section2D();
+s1.E     = 3.0e10;    % [N/m^2] Young's modulus for concrete
+s1.nue   = 0.2;       % [-]     Poisson's ratio for concrete
+s1.t     = 0.16;      % [m]     Section thickness
+s1.kappa = 5/6*1000;  % [-]     Shear correction factor
+%% Definition of element loads
+p1 = mafe.EleLoad( mafe.DofType.Disp3, -1.25e2 ); % vertical load in z [N/m^2], 
+%% 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, p1 );
+end
+%% Definition of constraints in the system
+const = mafe.Constraint.empty(0,0);
+% Edge Support at the left side
+nn = nodes( [gl6.getGridNodes.id, gl8.getGridNodes.id] ); % get nodes on GeoLine 6+8
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota2, mafe.ConstraintType.Dirichlet, 0.0 );
+% Edge Support at the top
+nn = nodes( [gl3.getGridNodes.id]); % get nodes on GeoLine 3
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota2, mafe.ConstraintType.Dirichlet, 0.0 );
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 );
+% Column Supports at gp4, gp5 and, gp8 (Displacement in z direction is not
+% allowed so Dirichlet constraint is used.
+nn = nodes( [gp4.getGridNodes.id] ); % get nodes at gp4
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 ); % disp3 representing z direction
+nn = nodes( [gp5.getGridNodes.id] ); % get nodes at gp5
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 ); % disp3 representing z direction
+nn = nodes( [gp8.getGridNodes.id] ); % get nodes at gp8
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 ); % disp3 representing z direction
+% Point Load at gp2 which is in the negative direction of Z axis with a
+% value of 25N
+nn = nodes( [gp2.getGridNodes.id] ); % get nodes at gp2
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Neumann, -25 ); % disp3 representing z direction

tutorials/statics/tutorial2/tutorial2_step3.m

+%% Example Plate bending
+% example on how to use the integrated mesh generator for
+% a specifically supported L-shaped plate domain under area and point loads
+%
+% step 3 (task 3.1-3.3):
+% 
+% create finite element problem and finite element analysis (object),
+% perform an analysis (only print the computed results)
+%
+addpath([ '..' filesep '..' filesep '..' filesep ]);
+%  -----------------------------------------------------------MESH GENERATION---
+%% Describe the goemetry of the problem for the mesh generator
+% helfpul parameters (problem specific)
+d = 2; n = 1; % d: length measure; n = number of grid nodes per d
+% function that controls grid node distribution along geo lines
+nf1 = mgen.MgNormFunction( mgen.MgNormFunctionType.Pow   , [0.0, 1.0] );
+nf2 = mgen.MgNormFunction( mgen.MgNormFunctionType.PowRev, [0.0, 1.0] );
+% create all geo points %% geo point coordinates specified in brackets []
+gp1 = mgen.MgGeoPoint( [0*d, 0*d] );
+gp2 = mgen.MgGeoPoint( [4*d, 0*d] );
+gp3 = mgen.MgGeoPoint( [0*d, 3*d] );
+gp4 = mgen.MgGeoPoint( [4*d, 3*d] );
+gp5 = mgen.MgGeoPoint( [8*d, 3*d] );
+gp6 = mgen.MgGeoPoint( [0*d, 8*d] );
+gp7 = mgen.MgGeoPoint( [4*d, 8*d] );
+gp8 = mgen.MgGeoPoint( [8*d, 8*d] );
+% create all geo lines %% geo line identified as continuous line between specified geo points, n presenting number of nodes per element
+gl1  = mgen.MgGeoLine2Point( [gp1, gp2], 4*n, nf2 );
+gl2  = mgen.MgGeoLine2Point( [gp3, gp4],   2, nf2 );
+gl3  = mgen.MgGeoLine2Point( [gp6, gp7],   2, nf2 );
+gl4  = mgen.MgGeoLine2Point( [gp4, gp5], 4*n, nf1 );
+gl5  = mgen.MgGeoLine2Point( [gp7, gp8],   2, nf1 );
+gl6  = mgen.MgGeoLine2Point( [gp1, gp3], 3*n, nf2 );
+gl7  = mgen.MgGeoLine2Point( [gp2, gp4],   2, nf2 );
+gl8  = mgen.MgGeoLine2Point( [gp3, gp6], 5*n, nf1 );
+gl9  = mgen.MgGeoLine2Point( [gp4, gp7],   2, nf1 );
+gl10 = mgen.MgGeoLine2Point( [gp8, gp5],   2, nf2 );
+% create all geo areas %% geo area identified as area enclosed by specified geo lines
+ga1 = mgen.MgGeoArea4Line( [gl1, gl7, gl2, gl6] );
+ga2 = mgen.MgGeoArea4Line( [gl2, gl9, gl3, gl8] );
+ga3 = mgen.MgGeoArea4Line( [gl4, gl10, gl5, gl9] );
+% create geo object and generate the mesh %% include all geo points, geo lines and geo areas
+geo = mgen.MgGeo( [gp1, gp2, gp3, gp4, gp5, gp6, gp7, gp8], ...
+                  [gl1, gl2, gl3, gl4, gl5, gl6, gl7, gl8, gl9, gl10], ...
+                  [ga1, ga2, ga3] );
+% print geo and mesh data info
+% fprintf(geo.print());
+% plot mesh
+cla; clf;
+fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on; %% this is used to properly position the figure
+geo.plot();
+%  -----------------------------------------------------------DEFINE FE ITEMS---
+%% Definition of sections
+s1 = mafe.Section2D();
+s1.E     = 3.0e10;    % [N/m^2] Young's modulus for concrete
+s1.nue   = 0.2;       % [-]     Poisson's ratio for concrete
+s1.t     = 0.16;      % [m]     Section thickness
+s1.kappa = 5/6*1000;  % [-]     Shear correction factor
+%% Definition of element loads
+p1 = mafe.EleLoad( mafe.DofType.Disp3, -1e3 ); % vertical load of p = 1 kN/m^2 in z [N/m^2]
+%% 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, p1 );
+end
+%% Definition of constraints in the system
+const = mafe.Constraint.empty(0,0);
+% Edge Support at the left side
+nn = nodes( [gl6.getGridNodes.id, gl8.getGridNodes.id] ); % get nodes on GeoLine 6+8
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota2, mafe.ConstraintType.Dirichlet, 0.0 );
+% Edge Support at the top
+nn = nodes( [gl3.getGridNodes.id]); % get nodes on GeoLine 3
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota2, mafe.ConstraintType.Dirichlet, 0.0 );
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 );
+% Column Supports at gp4, gp5 and, gp8 (Displacement in z direction is not
+% allowed so Dirichlet constraint is used.
+nn = nodes( [gp4.getGridNodes.id] ); % get nodes at gp4
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 ); % disp3 representing z direction
+nn = nodes( [gp5.getGridNodes.id] ); % get nodes at gp5
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 ); % disp3 representing z direction
+nn = nodes( [gp8.getGridNodes.id] ); % get nodes at gp8
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 ); % disp3 representing z direction
+% Point Load at gp2 which is in the negative direction of Z axis with a
+% value of 10kN
+nn = nodes( [gp2.getGridNodes.id] ); % get nodes at gp2
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Neumann, -10e3 ); % [N], disp3 representing z direction
+%  --------------------------------------------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) );
+%% Print numerical values to screen
+fep.printSystem();

tutorials/statics/tutorial2/tutorial2_step4.m

+%% Example Plate bending
+% example on how to use the integrated mesh generator for
+% a specifically supported L-shaped plate domain under area and point loads
+%
+% step 4 (task 4.1):
+% 
+% visualise the analysis results for the deformed state, the moment
+% states,the shear states the tensor field representation and, the von
+% Mises stress state.
+%
+addpath([ '..' filesep '..' filesep '..' filesep ]);
+%  -----------------------------------------------------------MESH GENERATION---
+%% Describe the goemetry of the problem for the mesh generator
+% helfpul parameters (problem specific)
+d = 2; n = 3; % d: length measure; n = number of grid nodes per d
+% function that controls grid node distribution along geo lines
+nf1 = mgen.MgNormFunction( mgen.MgNormFunctionType.Pow   , [0.0, 1.0] );
+nf2 = mgen.MgNormFunction( mgen.MgNormFunctionType.PowRev, [0.0, 1.0] );
+% create all geo points %% geo point coordinates specified in brackets []
+gp1 = mgen.MgGeoPoint( [0*d, 0*d] );
+gp2 = mgen.MgGeoPoint( [4*d, 0*d] );
+gp3 = mgen.MgGeoPoint( [0*d, 3*d] );
+gp4 = mgen.MgGeoPoint( [4*d, 3*d] );
+gp5 = mgen.MgGeoPoint( [8*d, 3*d] );
+gp6 = mgen.MgGeoPoint( [0*d, 8*d] );
+gp7 = mgen.MgGeoPoint( [4*d, 8*d] );
+gp8 = mgen.MgGeoPoint( [8*d, 8*d] );
+% create all geo lines %% geo line identified as continuous line between specified geo points, n presenting number of nodes per element
+gl1  = mgen.MgGeoLine2Point( [gp1, gp2], 4*n, nf2 );
+gl2  = mgen.MgGeoLine2Point( [gp3, gp4],   2, nf2 );
+gl3  = mgen.MgGeoLine2Point( [gp6, gp7],   2, nf2 );
+gl4  = mgen.MgGeoLine2Point( [gp4, gp5], 4*n, nf1 );
+gl5  = mgen.MgGeoLine2Point( [gp7, gp8],   2, nf1 );
+gl6  = mgen.MgGeoLine2Point( [gp1, gp3], 3*n, nf2 );
+gl7  = mgen.MgGeoLine2Point( [gp2, gp4],   2, nf2 );
+gl8  = mgen.MgGeoLine2Point( [gp3, gp6], 5*n, nf1 );
+gl9  = mgen.MgGeoLine2Point( [gp4, gp7],   2, nf1 );
+gl10 = mgen.MgGeoLine2Point( [gp8, gp5],   2, nf2 );
+% create all geo areas %% geo area identified as area enclosed by specified geo lines
+ga1 = mgen.MgGeoArea4Line( [gl1, gl7, gl2, gl6] );
+ga2 = mgen.MgGeoArea4Line( [gl2, gl9, gl3, gl8] );
+ga3 = mgen.MgGeoArea4Line( [gl4, gl10, gl5, gl9] );
+% create geo object and generate the mesh %% include all geo points, geo lines and geo areas
+geo = mgen.MgGeo( [gp1, gp2, gp3, gp4, gp5, gp6, gp7, gp8], ...
+                  [gl1, gl2, gl3, gl4, gl5, gl6, gl7, gl8, gl9, gl10], ...
+                  [ga1, ga2, ga3] );
+% print geo and mesh data info
+% fprintf(geo.print());
+% plot mesh
+% cla; clf;
+% fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on; %% this is used to properly position the figure
+% geo.plot();
+%  -----------------------------------------------------------DEFINE FE ITEMS---
+%% Definition of sections
+s1 = mafe.Section2D();
+s1.E     = 3.0e10;    % [N/m^2] Young's modulus for concrete
+s1.nue   = 0.2;       % [-]     Poisson's ratio for concrete
+s1.t     = 0.16;      % [m]     Section thickness
+s1.kappa = 5/6*1000;  % [-]     Shear correction factor
+%% Definition of element loads
+p1 = mafe.EleLoad( mafe.DofType.Disp3, -1e3 ); % vertical load of p = 1 kN/m^2 in z [N/m^2]
+%% 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, p1 );
+end
+%% Definition of constraints in the system
+const = mafe.Constraint.empty(0,0);
+% Edge Support at the left side
+nn = nodes( [gl6.getGridNodes.id, gl8.getGridNodes.id] ); % get nodes on GeoLine 6+8
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota2, mafe.ConstraintType.Dirichlet, 0.0 );
+% Edge Support at the top
+nn = nodes( [gl3.getGridNodes.id]); % get nodes on GeoLine 3
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota2, mafe.ConstraintType.Dirichlet, 0.0 );
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 );
+% Column Supports at gp4, gp5 and, gp8 (Displacement in z direction is not
+% allowed so Dirichlet constraint is used.
+nn = nodes( [gp4.getGridNodes.id] ); % get nodes at gp4
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 ); % disp3 representing z direction
+nn = nodes( [gp5.getGridNodes.id] ); % get nodes at gp5
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 ); % disp3 representing z direction
+nn = nodes( [gp8.getGridNodes.id] ); % get nodes at gp8
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 ); % disp3 representing z direction
+% Point Load at gp2 which is in the negative direction of Z axis with a
+% value of 10kN
+nn = nodes( [gp2.getGridNodes.id] ); % get nodes at gp2
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Neumann, -10e3 ); % [N], disp3 representing z direction
+%  --------------------------------------------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) );
+%% Print numerical values to screen
+% fep.printSystem(); % This option can be enabled to print the elements
+% and the moment and shear force values in each element
+%% Visualise system response
+cla; clf;
+fig = figure(1); clf; subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
+fep.plotSystem('reference');
+fep.plotSystem('deformed',1e02); % amplification by factor of 100 for visualisation
+title('Undeformed and deformed state')
+fig = figure(2); clf; subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
+fep.plotSystem('moment11'); colorbar; % with colorbar legend
+title('Moment_{11}')
+fig = figure(3); clf; subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
+fep.plotSystem('moment22'); colorbar; % with colorbar legend
+title('Moment_{22}')
+fig = figure(4); clf; subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
+fep.plotSystem('shear13'); colorbar; % with colorbar legend
+title('Shear_{13}')
+fig = figure(5); clf; subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
+fep.plotSystem('shear23',1e10); colorbar; % with colorbar legend
+title('Shear_{23}')
+fig = figure(6); clf; subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
+fep.plotSystem('vonmises'); colorbar; % with colorbar legend
+title('von Mises stress')
+fig = figure(7); clf; subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
+fep.plotSystem('tensor',2e-5); % deamplification by factor for proper presentation
+title('principal stress orientations')