Examples updated.

examples/dynamics/fluids/example_ins_cylinder.m

@@ -5,7 +5,7 @@
 addpath([ '..' filesep '..' filesep '..' filesep ]);
 %% Definition of sections
 s1 = mafe.Section2D();
-s1.rho   = 1000.0;    % [kg/m^3]
+s1.rho   = 1.0;    % [kg/m^3]
 s1.mu    = 1.0;    % [Ns/m^2]
 s1.t     = 1.0;    % [m]
 %  -----------------------------------------------------------MESH GENERATION---

examples/statics/plates/example_static_plate1b.m

@@ -53,7 +53,7 @@ case 'bend'
               mafe.ConstraintNode( n2, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 ),
               mafe.ConstraintNode( n7, mafe.DofType.Rota1, mafe.ConstraintType.Neumann,   1.0 ),
               mafe.ConstraintNode( n8, mafe.DofType.Rota1, mafe.ConstraintType.Neumann,   1.0 ),
-              ];
+              ]';
 case 'shear'
     % impose constant shear
     %
@@ -77,7 +77,7 @@ case 'shear'
               mafe.ConstraintNode( n8, mafe.DofType.Rota2, mafe.ConstraintType.Dirichlet, 0.0 ),
               mafe.ConstraintNode( n8, mafe.DofType.Disp3, mafe.ConstraintType.Neumann,   1.0 ),
               mafe.ConstraintNode( n7, mafe.DofType.Disp3, mafe.ConstraintType.Neumann,   1.0 ),
-              ];
+              ]';
 case 'twist'
     % impose constant twist
     %
@@ -85,7 +85,7 @@ case 'twist'
               mafe.ConstraintNode( n2, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 ),
               mafe.ConstraintNode( n8, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 ),
               mafe.ConstraintNode( n7, mafe.DofType.Disp3, mafe.ConstraintType.Neumann,   1.0 ),
-              ];
+              ]';
 otherwise
     disp('Unknown test case! Constraints not set');
     const = [];

examples/statics/plates/example_static_plate2b.m

@@ -13,7 +13,7 @@ p1 = mafe.EleLoad( mafe.DofType.Disp3, -1.e-3*s1.E*s1.t^3 ); % vertical load in
 %  -----------------------------------------------------------MESH GENERATION---
 %% Describe the goemetry of the problem for the mesh generator
 % helfpul parameters (problem specific)
-lx = 8.0; ly = 4.0; nx = 8*4; ny = nx*ly/lx; % lx, ly: total lengths ; nx, ny: number of elements per direction
+lx = 8.0; ly = 4.0; nx = 4*4; ny = nx*ly/lx; % lx, ly: total lengths ; nx, ny: number of elements per direction
 % function that controls grid node distribution along geo lines
 nf1 = mgen.MgNormFunction( mgen.MgNormFunctionType.Pow   , [0.0, 0.1] );
 nf2 = mgen.MgNormFunction( mgen.MgNormFunctionType.PowRev, [0.0, 0.1] );
@@ -49,9 +49,11 @@ const = mafe.Constraint.empty(0,0);
 % left edge: vertical support (Navier-type support)
 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.Rota2, mafe.ConstraintType.Dirichlet, 0.0 ); % hard support (Kirchhoff)
 % lower edge: vertical support (Navier-type support)
 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 );
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 ); % hard support (Kirchhoff)
 % right edge: bending angle Rota1!=0 (because of symmetry about y-axis, 1/4 plate)
 nn = nodes( [gl2.getGridNodes.id] ); % get nodes on GeoLine 2 (right)
 const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 );
@@ -86,7 +88,7 @@ G = ( s1.E / (2*(1.0 + s1.nue)) ); % shear modulus
 x = lx/2; y = ly/2; % point of evaluation: at centre
 
 for m = 1:2:tx %
-    for n = 1:2:ty % 
+    for n = 1:2:ty %
         % coefficients from load function
         p_mn = 16 * p1.p(mafe.DofType.Disp3)/(pi^2*m*n); % constant area load p1
         % series contribution to transversal displacement and derivatives
@@ -130,10 +132,10 @@ fprintf( '%6s%30s\t%30s\t%30s\n', sprintf('u_z'), ...
                                   sprintf('%+4.8e', u_z_c), ...
                                   sprintf('%4.8e' , u_z_e) );
 fprintf( '%6s%30s\t%30s\t%30s\n', sprintf('m_xx'), ...
-                                  sprintf('%+4.8e', m_xx_a), ...
+                                  sprintf('[Kirchhoff] %+4.8e', m_xx_a), ...
                                   sprintf('%+4.8e', m_xx_c), ...
                                   sprintf('%4.8e' , m_xx_e) );
 fprintf( '%6s%30s\t%30s\t%30s\n', sprintf('m_yy'), ...
-                                  sprintf('%+4.8e', m_yy_a), ...
+                                  sprintf('[Kirchhoff] %+4.8e', m_yy_a), ...
                                   sprintf('%+4.8e', m_yy_c), ...
                                   sprintf('%4.8e' , m_yy_e) );

examples/statics/plates/example_static_plate2c.m

@@ -13,7 +13,7 @@ F = -1.e-3*s1.E*s1.t^3;
 %  -----------------------------------------------------------MESH GENERATION---
 %% Describe the goemetry of the problem for the mesh generator
 % helfpul parameters (problem specific)
-lx = 8.0; ly = 4.0; nx = 8*4; ny = nx*ly/lx; % lx, ly: total lengths ; nx, ny: number of elements per direction
+lx = 8.0; ly = 4.0; nx = 4*4; ny = nx*ly/lx; % lx, ly: total lengths ; nx, ny: number of elements per direction
 % function that controls grid node distribution along geo lines
 nf1 = mgen.MgNormFunction( mgen.MgNormFunctionType.Pow   , [0.0, 0.1] );
 nf2 = mgen.MgNormFunction( mgen.MgNormFunctionType.PowRev, [0.0, 0.1] );
@@ -49,9 +49,11 @@ const = mafe.Constraint.empty(0,0);
 % left edge: vertical support (Navier-type support)
 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.Rota2, mafe.ConstraintType.Dirichlet, 0.0 ); % hard support (Kirchhoff)
 % lower edge: vertical support (Navier-type support)
 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 );
+const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 ); % hard support (Kirchhoff)
 % right edge: bending angle Rota1!=0 (because of symmetry about y-axis, 1/4 plate)
 nn = nodes( [gl2.getGridNodes.id] ); % get nodes on GeoLine 2 (right)
 const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 );
@@ -133,10 +135,10 @@ fprintf( '%6s%30s\t%30s\t%30s\n', sprintf('u_z'), ...
                                   sprintf('%+4.8e', u_z_c), ...
                                   sprintf('%4.8e' , u_z_e) );
 fprintf( '%6s%30s\t%30s\t%30s\n', sprintf('m_xx'), ...
-                                  sprintf('%+4.8e', m_xx_a), ...
+                                  sprintf('[Kirchhoff] %+4.8e', m_xx_a), ...
                                   sprintf('%+4.8e', m_xx_c), ...
                                   sprintf('%4.8e' , m_xx_e) );
 fprintf( '%6s%30s\t%30s\t%30s\n', sprintf('m_yy'), ...
-                                  sprintf('%+4.8e', m_yy_a), ...
+                                  sprintf('[Kirchhoff] %+4.8e', m_yy_a), ...
                                   sprintf('%+4.8e', m_yy_c), ...
                                   sprintf('%4.8e' , m_yy_e) );

examples/statics/plates/example_static_plate_tpf1.m

-%% Example Plate Bending 1a
-% A rectangular plate structure composed of 1 element
-addpath([ '..' filesep '..' filesep '..' filesep ]);
-%% Definition of sections
-s1 = mafe.Section2D();
-s1.E     = 2000.0; % [N/m^2]
-s1.nue   = 0.3; % [-]
-s1.t     = 0.1; % [m]
-%% Definition of element loads
-p1 = mafe.EleLoad( mafe.DofType.Disp3, -0.75 );
-%% Definition of nodes in the system
-% Node 1
-n1 = mafe.Node2D( [ 4.0,  3.0] );
-% Node 2
-n2 = mafe.Node2D( [-5.0,  4.0] );
-% Node 3
-n3 = mafe.Node2D( [-3.0, -4.0] );
-% Node 4
-n4 = mafe.Node2D( [ 2.5, -3.5] );
-% Vector of nodes
-nodes = [ n1, n2, n3, n4 ];
-%% Definition of elements in the system
-% Element 1
-e1 = mafe.Plate2D( [n1, n2, n3, n4], s1, p1 );
-% vector of elements
-elems = [ e1 ];
-%% Definition of constraints in the system
-% Constraint 1 (fixed transversal displacement at left end, lower node)
-c1 = mafe.ConstraintNode( n1, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
-% Constraint 2 (fixed transversal displacement at left end, upper node)
-c2 = mafe.ConstraintNode( n4, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
-% Constraint 3 (fixed transversal displacement at left end, lower node)
-c3 = mafe.ConstraintNode( n1, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 );
-% Constraint 4 (fixed transversal displacement at left end, upper node)
-c4 = mafe.ConstraintNode( n4, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 );
-% Constraint 5 (horizontal force at right end, lower node)
-c5 = mafe.ConstraintNode( n2, mafe.DofType.Rota1, mafe.ConstraintType.Neumann,   0.1 ); % 0.1 [N]
-% Constraint 6 (horizontal force at right end, upper node)
-c6 = mafe.ConstraintNode( n3, mafe.DofType.Rota1, mafe.ConstraintType.Neumann,   0.1 ); % 0.1 [N]
-% vector of constraints
-const = [ c1, c2, c3, c4, c5, c6 ];
-%% 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();
-%% 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', 1.e0);
-%fep.plotSystem('tensor');
-fep.printSystem();

examples/statics/plates/example_static_plate_tpf2.m

-%% Example Plate Bending 2b
-% example on how to use the integrated mesh generator for
-% a rectangular plate with Navier-type support and constant area load, 1/4 plate
-% in addition, the FE solution is compared with the analytical solution at midpoint
-addpath([ '..' filesep '..' filesep '..' filesep ]);
-%% Definition of sections
-s1 = mafe.Section2D();
-s1.E   = 2000; % [N/m^2] | concrete (compression)
-s1.nue = 0.30;   % [-]
-s1.t   = 0.01;   % [m]
-s1.kappa = 5/6*1;
-%% Definition of element loads
-%p1 = mafe.EleLoad( mafe.DofType.Disp3, -1.e-3*s1.E*s1.t^3 ); % vertical load in z [N/m^2]
-p1 = mafe.EleLoad( mafe.DofType.Disp3, -0.75 ); % vertical load in z [N/m^2]
-%  -----------------------------------------------------------MESH GENERATION---
-%% Describe the goemetry of the problem for the mesh generator
-% helfpul parameters (problem specific)
-lx = 3.0; ly = 2.0; nx = 30; ny = nx*ly/lx; % lx, ly: total lengths ; nx, ny: number of elements per direction
-% function that controls grid node distribution along geo lines
-nf1 = mgen.MgNormFunction( mgen.MgNormFunctionType.Pow   , [0.0, 0.1] );
-nf2 = mgen.MgNormFunction( mgen.MgNormFunctionType.PowRev, [0.0, 0.1] );
-% create all geo points
-gp1 = mgen.MgGeoPoint( [ 0.0,  0.0] );
-gp2 = mgen.MgGeoPoint( [lx/2,  0.0] );
-gp3 = mgen.MgGeoPoint( [lx/2, ly/2] );
-gp4 = mgen.MgGeoPoint( [ 0.0, ly/2] );
-% create all geo lines
-gl1 = mgen.MgGeoLine2Point( [gp2, gp1],  nx+1 ); % ask for nx+1 grid nodes along this line
-gl2 = mgen.MgGeoLine2Point( [gp3, gp2],  ny+1 ); % ask for ny+1 grid nodes along this line
-gl3 = mgen.MgGeoLine2Point( [gp4, gp3],     2 ); % this will adapt automatically
-gl4 = mgen.MgGeoLine2Point( [gp1, gp4],     2 ); % 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, p1 );
-end
-%% Definition of constraints in the system
-const = mafe.Constraint.empty(0,0);
-% left edge: vertical support (Navier-type support)
-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.Rota2, mafe.ConstraintType.Dirichlet, 0.0 ); % for Kirchhoff solution in q
-% lower edge: vertical support (Navier-type support)
-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 );
-const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 ); % for Kirchhoff solution in q
-% right edge: bending angle Rota1!=0 (because of symmetry about y-axis, 1/4 plate)
-nn = nodes( [gl2.getGridNodes.id] ); % get nodes on GeoLine 2 (right)
-const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 );
-% upper edge: bending angle Rota2!=0 (because of symmetry about x-axis, 1/4 plate)
-nn = nodes( [gl3.getGridNodes.id] ); % get nodes on GeoLine 3 (top)
-const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota2, mafe.ConstraintType.Dirichlet, 0.0 );
-%  --------------------------------------------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', 1.e-3); view(45.0, 25.0);
-%fep.plotSystem('shear13'); colorbar();
-%fep.plotSystem('shear23'); colorbar();
-%fep.plotSystem('moment11'); colorbar();
-%fep.plotSystem('moment22'); colorbar();
-%fep.plotSystem('moment12'); colorbar();
-%fep.plotSystem('tensor', 1.e+2/(s1.E*s1.t^3));
-%fep.printSystem();
-% ------------------------------------------COMPUTE REFERENCE SOLUTION VALUES---
-
-% reference solution based on double fourier series
-u_z_a  = 0.0; m_xx_a = 0.0; m_yy_a = 0.0; % displacement, bending moments
-tx = 64*4; ty = tx/2; % series terms to be evaluated per direction
-
-B = s1.E * s1.t^3 / (12*(1-s1.nue^2)); % plate bending stiffness
-G = ( s1.E / (2*(1.0 + s1.nue)) ); % shear modulus
-
-x = lx/2; y = ly/2; % point of evaluation: at centre
-
-for m = 1:2:tx %
-    for n = 1:2:ty %
-        % coefficients from load function
-        p_mn = 16 * p1.p(mafe.DofType.Disp3)/(pi^2*m*n); % constant area load p1
-        % series contribution to transversal displacement and derivatives
-        zx = (m/lx)^2; zy = (n/ly)^2; zz = zx + zy;
-        f = p_mn / (B*pi^4*zz^2);
-        g = 1.0 + B*pi^2/(G*s1.kappa*s1.t) * zz;
-        %
-        s = sin(m*pi*x/lx) * sin(n*pi*y/ly);
-        %
-        w      =  f * g * s            ;
-        phix_x =  f * 1 * s * zx*(pi^2);
-        phiy_y =  f * 1 * s * zy*(pi^2);
-        % compute values
-        u_z_a  = u_z_a  + w;
-        m_xx_a = m_xx_a + B * (          phix_x + s1.nue * phiy_y );
-        m_yy_a = m_yy_a + B * ( s1.nue * phix_x +          phiy_y );
-    end
-end
-
-% obtain data from the FE analysis: access nodal solution, element solution
-centre_node = nodes( gp3.getGridNodes.id );
-centre_elem = elems( ny ); % by inspection
-
-u_z_c  = centre_node.lhs( centre_node.dof == mafe.DofType.Disp3 );
-m_xx_c = centre_elem.internal(1,2); % mxx (1) at node (2) (by inspection of mesh)
-m_yy_c = centre_elem.internal(2,2); % myy (2) at node (2) (by inspection of mesh)
-intf = centre_elem.internal
-% compute L2 error
-u_z_e  = sqrt( ( u_z_c  - u_z_a  )^2 / u_z_a^2  );
-m_xx_e = sqrt( ( m_xx_c - m_xx_a )^2 / m_xx_a^2 );
-m_yy_e = sqrt( ( m_yy_c - m_yy_a )^2 / m_yy_a^2 );
-
-% output
-fprintf( 'Results for central point (x1,x2) = (%2.2f,%2.2f):\n', lx/2, ly/2 );
-fprintf( '%6s%30s\t%30s\t%30s\n', '---', ...
-                                  sprintf('analytical (%3dx%3d)', tx, ty), ...
-                                  sprintf('computed (%4dx%4d)', nx, ny), ...
-                                  sprintf('relative L2 error') );
-fprintf( '%6s%30s\t%30s\t%30s\n', sprintf('u_z'), ...
-                                  sprintf('%+4.8e', u_z_a), ...
-                                  sprintf('%+4.8e', u_z_c), ...
-                                  sprintf('%4.8e' , u_z_e) );
-fprintf( '%6s%30s\t%30s\t%30s\n', sprintf('m_xx'), ...
-                                  sprintf('%+4.8e', m_xx_a), ...
-                                  sprintf('%+4.8e', m_xx_c), ...
-                                  sprintf('%4.8e' , m_xx_e) );
-fprintf( '%6s%30s\t%30s\t%30s\n', sprintf('m_yy'), ...
-                                  sprintf('%+4.8e', m_yy_a), ...
-                                  sprintf('%+4.8e', m_yy_c), ...
-                                  sprintf('%4.8e' , m_yy_e) );

examples/statics/slabs/example_static_membr1b.m

@@ -51,14 +51,14 @@ case 'e11'
               mafe.ConstraintNode( n2, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, -nodedisp ),
               mafe.ConstraintNode( n8, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, +nodedisp ),
               mafe.ConstraintNode( n7, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, +nodedisp ),
-              mafe.ConstraintNode( n2, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet,       0.0 ) ];
+              mafe.ConstraintNode( n2, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet,       0.0 ) ]';
 case 'e22'
     % impose vertical displacement, leading to stretching e22
     const = [ mafe.ConstraintNode( n1, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, +nodedisp ),
               mafe.ConstraintNode( n2, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, -nodedisp ),
               mafe.ConstraintNode( n8, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, -nodedisp ),
               mafe.ConstraintNode( n7, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, +nodedisp ),
-              mafe.ConstraintNode( n2, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet,       0.0 ) ];
+              mafe.ConstraintNode( n2, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet,       0.0 ) ]';
 case 'e12'
     % impose combined horizontal/vertical displacement, leading to shearing e12
     const = [ mafe.ConstraintNode( n1, DofType.Disp1, mafe.ConstraintType.Dirichlet, +nodedisp ),
@@ -68,7 +68,7 @@ case 'e12'
               mafe.ConstraintNode( n1, DofType.Disp2, mafe.ConstraintType.Dirichlet, -nodedisp ),
               mafe.ConstraintNode( n2, DofType.Disp2, mafe.ConstraintType.Dirichlet, -nodedisp ),
               mafe.ConstraintNode( n8, DofType.Disp2, mafe.ConstraintType.Dirichlet, +nodedisp ),
-              mafe.ConstraintNode( n7, DofType.Disp2, mafe.ConstraintType.Dirichlet, +nodedisp ) ];
+              mafe.ConstraintNode( n7, DofType.Disp2, mafe.ConstraintType.Dirichlet, +nodedisp ) ]';
 otherwise
     disp('Unknown test case! Constraints not set');
     const = [];