SolidsPy reference

SolidsPy have the following modules:

  • solids_GUI.py: The main program;
  • preprocesor.py: Pre-processing subroutines including Gmsh convertion functions using meshio
  • assemutil.py: Assembly of elemental stiffness matrices ;
  • femutil.py: Shape functions, its derivatives and general finite element method subroutines;
  • uelutil.py: Elemental or local matrix subroutines for different elements; and
  • postprocesor.py: Several results handling subroutines.

solids_GUI: simple interface

Computes the displacement solution for a finite element assembly of 2D solids under point loads using as input easy-to-create text files containing element, nodal, materials and loads data. The input files are created out of a Gmsh (.msh) generated file using the Python module meshio.

Created by Juan Gomez and Nicolas Guarin-Zapata as part of the courses:

  • IC0283 Computational Modeling
  • IC0602 Introduction to the Finite Element Method

Which are part of the Civil Engineering Department at Universidad EAFIT.

solids_GUI.solids_GUI(plot_contours=True, compute_strains=False, folder=None)[source]

Run a complete workflow for a Finite Element Analysis

Parameters:
  • plot_contours (Bool (optional)) – Boolean variable to plot contours of the computed variables. By default it is True.
  • compute_strains (Bool (optional)) – Boolean variable to compute Strains and Stresses at nodes. By default it is False.
  • folder (string (optional)) – String with the path to the input files. If not provided it would ask for it in a pop-up window.
Returns:

  • UC (ndarray (nnodes, 2)) – Displacements at nodes.
  • E_nodes (ndarray (nnodes, 3), optional) – Strains at nodes. It is returned when compute_strains is True.
  • S_nodes (ndarray (nnodes, 3), optional) – Stresses at nodes. It is returned when compute_strains is True.

Assembly routines

Functions to assemble the system of equations for the Finite Element Analysis.

assemutil.DME(nodes, elements)[source]

Counts active equations, creates BCs array IBC[] and the assembly operator DME[]

Parameters:
  • nodes (ndarray.) – Array with the nodal numbers and coordinates.
  • elements (ndarray) – Array with the number for the nodes in each element.
Returns:

  • DME (ndarray (int)) – Assembly operator.
  • IBC (ndarray (int)) – Boundary conditions array.
  • neq (int) – Number of active equations in the system.
assemutil.assembler(elements, mats, nodes, neq, DME, sparse=True, uel=None)[source]

Assembles the global stiffness matrix

Parameters:
  • elements (ndarray (int)) – Array with the number for the nodes in each element.
  • mats (ndarray (float)) – Array with the material profiles.
  • nodes (ndarray (float)) – Array with the nodal numbers and coordinates.
  • DME (ndarray (int)) – Assembly operator.
  • neq (int) – Number of active equations in the system.
  • sparse (boolean (optional)) – Boolean variable to pick sparse assembler. It is True by default.
  • uel (callable function (optional)) – Python function that returns the local stiffness matrix.
Returns:

KG – Array with the global stiffness matrix. It might be dense or sparse, depending on the value of _sparse_
Return type:

ndarray (float)
assemutil.dense_assem(elements, mats, nodes, neq, DME, uel=None)[source]

Assembles the global stiffness matrix _KG_ using a dense storing scheme

Parameters:
  • elements (ndarray (int)) – Array with the number for the nodes in each element.
  • mats (ndarray (float)) – Array with the material profiles.
  • nodes (ndarray (float)) – Array with the nodal numbers and coordinates.
  • DME (ndarray (int)) – Assembly operator.
  • neq (int) – Number of active equations in the system.
  • uel (callable function (optional)) – Python function that returns the local stiffness matrix.
Returns:

KG – Array with the global stiffness matrix in a dense numpy array.
Return type:

ndarray (float)
assemutil.eqcounter(nodes)[source]

Counts active equations and creates BCs array IBC

Parameters:nodes (ndarray) – Array with nodes coordinates and boundary conditions.
Returns:
  • neq (int) – Number of equations in the system after removing the nodes with imposed displacements.
  • IBC (ndarray (int)) – Array that maps the nodes with number of equations.
assemutil.loadasem(loads, IBC, neq)[source]

Assembles the global Right Hand Side Vector RHSG

Parameters:
  • loads (ndarray) – Array with the loads imposed in the system.
  • IBC (ndarray (int)) – Array that maps the nodes with number of equations.
  • neq (int) – Number of equations in the system after removing the nodes with imposed displacements.
Returns:

RHSG – Array with the right hand side vector.
Return type:

ndarray
assemutil.retriever(elements, mats, nodes, i, uel=None)[source]

Computes the elemental stiffness matrix of element i

Parameters:
  • elements (ndarray) – Array with the number for the nodes in each element.
  • mats (ndarray.) – Array with the material profiles.
  • nodes (ndarray.) – Array with the nodal numbers and coordinates.
  • i (int.) – Identifier of the element to be assembled.
Returns:

  • kloc (ndarray (float)) – Array with the local stiffness matrix.
  • ndof (int.) – Number of degrees of fredom of the current element.
assemutil.sparse_assem(elements, mats, nodes, neq, DME, uel=None)[source]

Assembles the global stiffness matrix _KG_ using a sparse storing scheme

The scheme used to assemble is COOrdinate list (COO), and it converted to Compressed Sparse Row (CSR) afterward for the solution phase [1]_.

Parameters:
  • elements (ndarray (int)) – Array with the number for the nodes in each element.
  • mats (ndarray (float)) – Array with the material profiles.
  • nodes (ndarray (float)) – Array with the nodal numbers and coordinates.
  • DME (ndarray (int)) – Assembly operator.
  • neq (int) – Number of active equations in the system.
  • uel (callable function (optional)) – Python function that returns the local stiffness matrix.
Returns:

KG – Array with the global stiffness matrix in a sparse Compressed Sparse Row (CSR) format.
Return type:

ndarray (float)

References

[1]Sparse matrix. (2017, March 8). In Wikipedia, The Free Encyclopedia. https://en.wikipedia.org/wiki/Sparse_matrix

FEM routines

Functions to compute kinematics variables for the Finite Element Analysis.

The elements included are:
  1. 4 node bilinear quadrilateral.
  2. 6 node quadratic triangle.
  3. 3 node linear triangle.

The notation used is similar to the one used by Bathe [1]_.

References

[1]Bathe, Klaus-Jürgen. Finite element procedures. Prentice Hall, Pearson Education, 2006.
femutil.eletype(iet)[source]

Assigns number to degrees of freedom

According to iet assigns number of degrees of freedom, number of nodes and minimum required number of integration points.

Parameters:iet (int) –
Type of element. These are:
  1. 4 node bilinear quadrilateral.
  2. 6 node quadratic triangle.

3. 3 node linear triangle. 5. 2 node spring. 6. 2 node truss element. 7. 2 node beam (3 DOF per node).
Returns:
  • ndof (int) – Number of degrees of freedom for the selected element.
  • nnodes (int) – Number of nodes for the selected element.
  • ngpts (int) – Number of Gauss points for the selected element.
femutil.jacoper(dhdx, coord)[source]

Compute the Jacobian of the transformation evaluated at the Gauss point

Parameters:
  • dhdx (ndarray) – Derivatives of the interpolation function with respect to the natural coordinates.
  • coord (ndarray) – Coordinates of the nodes of the element (nn, 2).
Returns:

jaco_inv – Jacobian of the transformation evaluated at (r, s).
Return type:

ndarray (2, 2)
femutil.sha3(x, y)[source]

Shape functions for a 3-noded triangular element

Parameters:
  • x (float) – x coordinate for a point within the element.
  • y (float) – y coordinate for a point within the element.
Returns:

N – Array of interpolation functions.
Return type:

Numpy array

Examples

We can check evaluating at two different points, namely (0, 0) and (0, 0.5). Thus

>>> N = sha3(0, 0)
>>> N_ex = np.array([
...    [1, 0, 0, 0, 0, 0],
...    [0, 1, 0, 0, 0, 0]])
>>> np.allclose(N, N_ex)
True

and

>>> N = sha3(1/2, 1/2)
>>> N_ex = np.array([
...    [0, 0, 1/2, 0, 1/2, 0],
...    [0, 0, 0, 1/2, 0, 1/2]])
>>> np.allclose(N, N_ex)
True
femutil.sha4(x, y)[source]

Shape functions for a 4-noded quad element

Parameters:
  • x (float) – x coordinate for a point within the element.
  • y (float) – y coordinate for a point within the element.
Returns:

N – Array of interpolation functions.
Return type:

Numpy array

Examples

We can check evaluating at two different points, namely (0, 0) and (1, 1). Thus

>>> N = sha4(0, 0)
>>> N_ex = np.array([
...    [1/4, 0, 1/4, 0, 1/4, 0, 1/4, 0],
...    [0, 1/4, 0, 1/4, 0, 1/4, 0, 1/4]])
>>> np.allclose(N, N_ex)
True

and

>>> N = sha4(1, 1)
>>> N_ex = np.array([
...    [0, 0, 0, 0, 1, 0, 0, 0],
...    [0, 0, 0, 0, 0, 1, 0, 0]])
>>> np.allclose(N, N_ex)
True
femutil.sha6(x, y)[source]

Shape functions for a 6-noded triangular element

Parameters:
  • x (float) – x coordinate for a point within the element.
  • y (float) – y coordinate for a point within the element.
Returns:

N – Array of interpolation functions.
Return type:

Numpy array

Examples

We can check evaluating at two different points, namely (0, 0) and (0.5, 0.5). Thus

>>> N = sha6(0, 0)
>>> N_ex = np.array([
...    [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
...    [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
>>> np.allclose(N, N_ex)
True

and

>>> N = sha6(1/2, 1/2)
>>> N_ex = np.array([
...     [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
...     [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0]])
>>> np.allclose(N, N_ex)
True
femutil.stdm3NT(r, s, coord)[source]

Strain-displacement interpolator B for a 3-noded triang element

Parameters:
  • r (float) – r component in the natural space.
  • s (float) – s component in the natural space.
  • coord (ndarray) – Coordinates of the nodes of the element (3, 2).
Returns:

  • det (float) – Determinant evaluated at (r, s).
  • B (ndarray) – Strain-displacement interpolator evaluated at (r, s).
femutil.stdm4NQ(r, s, coord)[source]

Strain-displacement interpolator B for a 4-noded quad element

Parameters:
  • r (float) – r component in the natural space.
  • s (float) – s component in the natural space.
  • coord (ndarray) – Coordinates of the nodes of the element (4, 2).
Returns:

  • ddet (float) – Determinant evaluated at (r, s).
  • B (ndarray) – Strain-displacement interpolator evaluated at (r, s).
femutil.stdm6NT(r, s, coord)[source]

Strain-displacement interpolator B for a 6-noded triang element

Parameters:
  • r (float) – r component in the natural space.
  • s (float) – s component in the natural space.
  • coord (ndarray) – Coordinates of the nodes of the element (6, 2).
Returns:

  • ddet (float) – Determinant evaluated at (r, s).
  • B (ndarray) – Strain-displacement interpolator evaluated at (r, s).
femutil.str_el3(coord, ul)[source]

Compute the strains at each element integration point

This one is used for 3-noded triangular elements.

Parameters:
  • coord (ndarray) – Coordinates of the nodes of the element (nn, 2).
  • ul (ndarray) – Array with displacements for the element.
Returns:

  • epsGT (ndarray) – Strain components for the Gauss points.
  • xl (ndarray) – Configuration of the Gauss points after deformation.
femutil.str_el4(coord, ul)[source]

Compute the strains at each element integration point

This one is used for 4-noded quadrilateral elements.

Parameters:
  • coord (ndarray) – Coordinates of the nodes of the element (4, 2).
  • ul (ndarray) – Array with displacements for the element.
Returns:

  • epsGT (ndarray) – Strain components for the Gauss points.
  • xl (ndarray) – Configuration of the Gauss points after deformation.
femutil.str_el6(coord, ul)[source]

Compute the strains at each element integration point

This one is used for 6-noded triangular elements.

Parameters:
  • coord (ndarray) – Coordinates of the nodes of the element (6, 2).
  • ul (ndarray) – Array with displacements for the element.
Returns:

  • epsGT (ndarray) – Strain components for the Gauss points.
  • xl (ndarray) – Configuration of the Gauss points after deformation.
femutil.umat(nu, E)[source]

2D Elasticity consitutive matrix in plane stress

For plane strain use effective properties.

Parameters:
  • nu (float) – Poisson coefficient (-1, 0.5).
  • E (float) – Young modulus (>0).
Returns:

C – Constitutive tensor in Voigt notation.
Return type:

ndarray

Examples

>>> C = umat(1/3, 8/3)
>>> C_ex = np.array([
...    [3, 1, 0],
...    [1, 3, 0],
...    [0, 0, 1]])
>>> np.allclose(C, C_ex)
True

Numeric integration routines

Weights and coordinates for Gauss-Legendre quadrature [1]_. The values for triangles is presented in section 5.5 of Bathe book [2].

References

[1]Wikipedia contributors. “Gaussian quadrature.” Wikipedia, The Free Encyclopedia, 2 Nov. 2015. Web. 25 Dec. 2015. url: https://en.wikipedia.org/wiki/Gaussian_quadrature
[2]Bathe, Klaus-Jürgen. Finite element procedures. Prentice Hall, Pearson Education, 2006.
gaussutil.gpoints2x2()[source]

Gauss points for a 2 by 2 grid

Returns:
  • xw (ndarray) – Weights for the Gauss-Legendre quadrature.
  • xp (ndarray) – Points for the Gauss-Legendre quadrature.
gaussutil.gpoints3()[source]

Gauss points for a triangle element (3 points)

Returns:
  • xw (ndarray) – Weights for the Gauss-Legendre quadrature.
  • xp (ndarray) – Points for the Gauss-Legendre quadrature.
gaussutil.gpoints7()[source]

Gauss points for a triangle (7 points)

Returns:
  • xw (ndarray) – Weights for the Gauss-Legendre quadrature.
  • xp (ndarray) – Points for the Gauss-Legendre quadrature.

Preprocessor subroutines

This module contains functions to preprocess the input files to compute a Finite Element Analysis.

preprocesor.boundary_conditions(cells, cell_data, phy_lin, nodes_array, bc_x, bc_y)[source]

Impose nodal point boundary conditions as required by SolidsPy

Parameters:
  • cell (dictionary) – Dictionary created by meshio with cells information.
  • cell_data (dictionary) – Dictionary created by meshio with cells data information.
  • phy_lin (int) – Physical line where BCs are to be imposed.
  • nodes_array (int) – Array with the nodal data and to be modified by BCs.
  • bc_y (bc_x,) –
    Boundary condition flag along the x and y direction:
    • -1: restrained
    • 0: free
Returns:

nodes_array – Array with the nodal data after imposing BCs according to SolidsPy.
Return type:

int
preprocesor.echomod(nodes, mats, elements, loads, folder='')[source]

Create echoes of the model input files

preprocesor.ele_writer(cells, cell_data, ele_tag, phy_sur, ele_type, mat_tag, nini)[source]

Extracts a subset of elements from a complete mesh according to the physical surface phy_sur and writes down the proper fields into an elements array.

Parameters:
  • cell (dictionary) – Dictionary created by meshio with cells information.
  • cell_data (dictionary) – Dictionary created by meshio with cells data information.
  • ele_tag (string) – Element type according to meshio convention, e.g., quad9 or line3.
  • phy_sur (int) – Physical surface for the subset.
  • ele_type (int) – Element type.
  • mat_tag (int) – Material profile for the subset.
  • ndof (int) – Number of degrees of freedom for the elements.
  • nnode (int) – Number of nodes for the element.
  • nini (int) – Element id for the first element in the set.
Returns:

  • nf (int) – Element id for the last element in the set
  • els_array (int) – Elemental data.
preprocesor.initial_params()[source]

Read initial parameters for the simulation

The parameters to be read are:

  • folder: location of the input files.
  • name: name for the output files (if echo is True).
  • echo: echo output files.
preprocesor.loading(cells, cell_data, phy_lin, P_x, P_y)[source]

Impose nodal boundary conditions as required by SolidsPy

Parameters:
  • cell (dictionary) – Dictionary created by meshio with cells information.
  • cell_data (dictionary) – Dictionary created by meshio with cells data information.
  • phy_lin (int) – Physical line where BCs are to be imposed.
  • nodes_array (int) – Array with the nodal data and to be modified by BCs.
  • P_y (P_x,) – Load components in x and y directions.
Returns:

nodes_array – Array with the nodal data after imposing BCs according to SolidsPy.
Return type:

int
preprocesor.node_writer(points, point_data)[source]

Write nodal data as required by SolidsPy

Parameters:
  • points (dictionary) – Nodal points
  • point_data (dictionary) – Physical data associatted to the nodes.
Returns:

nodes_array – Array with the nodal data according to SolidsPy.
Return type:

ndarray (int)
preprocesor.readin(folder='')[source]

Read the input files

preprocesor.rect_grid(length, height, nx, ny, eletype=None)[source]

Generate a structured mesh for a rectangle

The coordinates of the nodes will be defined in the domain [-length/2, length/2] x [-height/2, height/2].

Parameters:
  • length (float) – Length of the domain.
  • height (gloat) – Height of the domain.
  • nx (int) – Number of elements in the x direction.
  • ny (int) – Number of elements in the y direction.
  • eletype (None) – It does nothing right now.
Returns:

  • x (ndarray (float)) – x-coordinates for the nodes.
  • y (ndarray (float)) – y-coordinates for the nodes.
  • els (ndarray) – Array with element data.

Examples

>>> x, y, els = rect_grid(2, 2, 2, 2)
>>> x
array([-1.,  0.,  1., -1.,  0.,  1., -1.,  0.,  1.])
>>> y
array([-1., -1., -1.,  0.,  0.,  0.,  1.,  1.,  1.])
>>> els
array([[0, 1, 0, 0, 1, 4, 3],
       [1, 1, 0, 1, 2, 5, 4],
       [2, 1, 0, 3, 4, 7, 6],
       [3, 1, 0, 4, 5, 8, 7]])

Postprocessor subroutines

This module contains functions to postprocess results.

postprocesor.complete_disp(IBC, nodes, UG)[source]

Fill the displacement vectors with imposed and computed values.

IBC : ndarray (int)
IBC (Indicator of Boundary Conditions) indicates if the nodes has any type of boundary conditions applied to it.
UG : ndarray (float)
Array with the computed displacements.
nodes : ndarray (float)
Array with number and nodes coordinates
Returns:UC – Array with the displacements.
Return type:(nnodes, 2) ndarray (float)
postprocesor.eigvals(mat, tol=1e-06)[source]

Eigenvalues and eigenvectors for a 2x2 symmetric matrix/tensor

Parameters:
  • mat (ndarray) – Symmetric matrix.
  • tol (float (optional)) – Tolerance for considering a matrix diagonal.
Returns:

  • eig1 (float) – First eigenvalue.
  • eig2 (float) – Second eigenvalue.
  • vec1 (ndarray) – First eigenvector.
  • vec2 (ndarray) – Second eigenvector

Examples

>>> mat = np.array([[5, 6],
...              [6, 9]])
>>> eig1, eig2, vec1, vec2 =  eigvals(mat)
>>> np.allclose(eig1, 7 + 2*np.sqrt(10))
True
>>> np.allclose(eig2, 7 - 2*np.sqrt(10))
True
>>> np.allclose(vec1, np.array([-0.584710284663765, -0.8112421851755609]))
True
>>> np.allclose(vec2, np.array([-0.8112421851755609,0.584710284663765]))
True
postprocesor.energy(disp, stiff)[source]

Computes the potential energy for the current solution.

Parameters:
  • disp (ndarray (float)) – Array with the computed displacements.
  • stiff (ndarray (float)) – Global stiffness matrix.
Returns:

el_energy – Total energy in the system. \(-\frac{1}{2} U^T K U\)
Return type:

scalar (float)
postprocesor.fields_plot(elements, nodes, disp, E_nodes=None, S_nodes=None)[source]

Plot contours for displacements, strains and stresses

Parameters:
  • nodes (ndarray (float)) –
    Array with number and nodes coordinates:
    number coordX coordY BCX BCY
  • elements (ndarray (int)) – Array with the node number for the nodes that correspond to each element.
  • disp (ndarray (float)) – Array with the displacements.
  • E_nodes (ndarray (float)) – Array with strain field in the nodes.
  • S_nodes (ndarray (float)) – Array with stress field in the nodes.
postprocesor.mesh2tri(nodes, elements)[source]

Generate a matplotlib.tri.Triangulation object from the mesh

Parameters:
  • nodes (ndarray (float)) –
    Array with number and nodes coordinates:
    number coordX coordY BCX BCY
  • elements (ndarray (int)) – Array with the node number for the nodes that correspond to each element.
Returns:

tri – An unstructured triangular grid consisting of npoints points and ntri triangles.
Return type:

Triangulation
postprocesor.plot_node_field(field, nodes, elements, plt_type='contourf', levels=12, savefigs=False, title=None, figtitle=None, filename=None)[source]

Plot the nodal displacement using a triangulation

Parameters:
  • field (ndarray (float)) – Array with the field to be plotted. The number of columns determine the number of plots.
  • nodes (ndarray (float)) – Array with number and nodes coordinates number coordX coordY BCX BCY
  • elements (ndarray (int)) – Array with the node number for the nodes that correspond to each element.
  • plt_type (string (optional)) – Plot the field as one of the options: pcolor or contourf.
  • levels (int (optional)) – Number of levels to be used in contourf.
  • savefigs (bool (optional)) – Allow to save the figure.
  • title (Tuple of strings (optional)) – Titles of the plots. If not provided the plots will not have a title.
  • figtitle (Tuple of strings (optional)) – Titles of the plotting windows. If not provided the windows will not have a title.
  • filename (Tuple of strings (optional)) – Filenames to save the figures. Only used when savefigs=True. If not provided the name of the figures would be “outputk.pdf”, where k is the number of the column.
postprocesor.plot_truss(nodes, elements, mats, stresses=None, max_val=4, min_val=0.5, savefigs=False, title=None, figtitle=None, filename=None)[source]

Plot a truss and encodes the stresses in a colormap

Parameters:
  • UC ((nnodes, 2) ndarray (float)) – Array with the displacements.
  • nodes (ndarray (float)) – Array with number and nodes coordinates number coordX coordY BCX BCY
  • elements (ndarray (int)) – Array with the node number for the nodes that correspond to each element.
  • mats (ndarray (float)) – Array with material profiles.
  • loads (ndarray (float)) – Array with loads.
  • tol (float (optional)) – Minimum difference between cross-section areas of the members to be considered different.
  • savefigs (bool (optional)) – Allow to save the figure.
  • title (Tuple of strings (optional)) – Titles of the plots. If not provided the plots will not have a title.
  • figtitle (Tuple of strings (optional)) – Titles of the plotting windows. If not provided the windows will not have a title.
  • filename (Tuple of strings (optional)) – Filenames to save the figures. Only used when savefigs=True. If not provided the name of the figures would be “outputk.pdf”, where k is the number of the column.
postprocesor.principal_dirs(field)[source]

Compute the principal directions of a tensor field

Parameters:field (ndarray) – Tensor field. The tensor is written as “vector” using Voigt notation.
Returns:
  • eigs1 (ndarray) – Array with the first eigenvalues.
  • eigs2 (ndarray) – Array with the second eigenvalues.
  • vecs1 (ndarray) – Array with the first eigenvectors.
  • vecs2 (ndarray) – Array with the Second eigenvector.
postprocesor.strain_nodes(nodes, elements, mats, UC)[source]

Compute averaged strains and stresses at nodes

First, the variable is extrapolated from the Gauss point to nodes for each element. Then, these values are averaged according to the number of element that share that node. The theory for this technique can be found in [1]_.

Parameters:
  • nodes (ndarray (float)) – Array with nodes coordinates.
  • elements (ndarray (int)) – Array with the node number for the nodes that correspond to each element.
  • mats (ndarray (float)) – Array with material profiles.
  • UC (ndarray (float)) – Array with the displacements. This one contains both, the computed and imposed values.
Returns:

E_nodes – Strains evaluated at the nodes.
Return type:

ndarray

References

[1]O.C. Zienkiewicz and J.Z. Zhu, The Superconvergent patch recovery and a posteriori error estimators. Part 1. The recovery technique, Int. J. Numer. Methods Eng., 33, 1331-1364 (1992).
postprocesor.stress_truss(nodes, elements, mats, disp)[source]

Compute axial stresses in truss members

Parameters:
  • nodes (ndarray (float)) – Array with nodes coordinates.
  • elements (ndarray (int)) – Array with the node number for the nodes that correspond to each element.
  • mats (ndarray (float)) – Array with material profiles.
  • disp (ndarray (float)) – Array with the displacements. This one contains both, the computed and imposed values.
Returns:

stress – Stresses for each member on the truss
Return type:

ndarray

Examples

The following examples are taken from [1]_. In all the examples \(A=1\), \(E=1\).

>>> import assemutil as ass
>>> import solidspy.solutil as sol

>>> def fem_sol(nodes, elements, mats, loads):
...     DME , IBC , neq = ass.DME(nodes, elements)
...     KG = ass.assembler(elements, mats, nodes, neq, DME)
...     RHSG = ass.loadasem(loads, IBC, neq)
...     UG = sol.static_sol(KG, RHSG)
...     UC = complete_disp(IBC, nodes, UG)
...     return UC

Exercise 3.3-18

The axial stresses in this example are

\[[\sigma] = \left[\frac{1}{2},\frac{\sqrt{3}}{4},\frac{1}{4}\right]\]
>>> nodes = np.array([
... [0, 0.0,  0.0, 0, -1],
... [1, -1.0,  0.0, -1, -1],
... [2,  -np.cos(np.pi/6),  -np.sin(np.pi/6),  -1,  -1],
... [3,  -np.cos(np.pi/3),  -np.sin(np.pi/3),  -1,  -1]])
>>> mats = np.array([[1.0, 1.0]])
>>> elements = np.array([
... [0, 6, 0, 1, 0],
... [1, 6, 0, 2, 0],
... [2, 6, 0, 3, 0]])
>>> loads = np.array([[0, 1.0, 0]])
>>> disp = fem_sol(nodes, elements, mats, loads)
>>> stress = stress_truss(nodes, elements, mats, disp)
>>> stress_exact = np.array([1/2, np.sqrt(3)/4, 1/4])
>>> np.allclose(stress_exact, stress)
True

Exercise 3.3-19

The axial stresses in this example are

\[[\sigma] = \left[\frac{1}{\sqrt{2}+2}, \frac{\sqrt{2}}{\sqrt{2}+1}, \frac{1}{\sqrt{2}+2}\right]\]
>>> nodes = np.array([
...    [0, 0.0,  0.0, 0, 0],
...    [1, -1.0,  -1.0, -1, -1],
...    [2,  0.0,  -1.0,  -1,  -1],
...    [3,  1.0, -1.0,  -1,  -1]])
>>> mats = np.array([[1.0, 1.0]])
>>> elements = np.array([
...    [0, 6, 0, 1, 0],
...    [1, 6, 0, 2, 0],
...    [2, 6, 0, 3, 0]])
>>> loads = np.array([[0, 0, 1]])
>>> disp = fem_sol(nodes, elements, mats, loads)
>>> stress = stress_truss(nodes, elements, mats, disp)
>>> stress_exact = np.array([
...     1/(np.sqrt(2) + 2),
...     np.sqrt(2)/(np.sqrt(2) + 1),
...     1/(np.sqrt(2) + 2)])
>>> np.allclose(stress_exact, stress)
True

Exercise 3.3-22

The axial stresses in this example are

\[[\sigma] =\left[\frac{1}{3 \sqrt{2}},\frac{5}{12}, \frac{1}{2^{\frac{3}{2}}}, \frac{1}{12}, -\frac{1}{3 \sqrt{2}}\right]\]
>>> cathetus = np.cos(np.pi/4)
>>> nodes = np.array([
...    [0, 0.0,  0.0, 0, 0],
...    [1, -1.0,  0.0, -1, -1],
...    [2,  -cathetus,  cathetus,  -1,  -1],
...    [3,  0.0, 1.0,  -1,  -1],
...    [4,  cathetus, cathetus,  -1,  -1],
...    [5,  1.0, 0.0,  -1,  -1]])
>>> mats = np.array([[1.0, 1.0]])
>>> elements = np.array([
...    [0, 6, 0, 1, 0],
...    [1, 6, 0, 2, 0],
...    [2, 6, 0, 3, 0],
...    [3, 6, 0, 4, 0],
...    [4, 6, 0, 5, 0]])
>>> loads = np.array([[0, cathetus, -cathetus]])
>>> disp = fem_sol(nodes, elements, mats, loads)
>>> stress = stress_truss(nodes, elements, mats, disp)
>>> stress_exact = np.array([
...     1/(3*np.sqrt(2)),
...     5/12,
...     1/2**(3/2),
...     1/12,
...     -1/(3*np.sqrt(2))])
>>> np.allclose(stress_exact, stress)
True

References

[1]William Weaver and James Gere. Matrix Analysis of Framed Structures. Third Edition, Van Nostrand Reinhold, New York (1990).
postprocesor.tri_plot(tri, field, title='', levels=12, savefigs=False, plt_type='contourf', filename='solution_plot.pdf')[source]

Plot contours over triangulation

Parameters:
  • tri (ndarray (float)) – Array with number and nodes coordinates: number coordX coordY BCX BCY
  • field (ndarray (float)) – Array with data to be plotted for each node.
  • title (string (optional)) – Title of the plot.
  • levels (int (optional)) – Number of levels to be used in contourf.
  • savefigs (bool (optional)) – Allow to save the figure.
  • plt_type (string (optional)) – Plot the field as one of the options: pcolor or contourf
  • filename (string (optional)) – Filename to save the figures.

Solver routines

Utilities for solution of FEM systems

solutil.static_sol(mat, rhs)[source]

Solve a static problem [mat]{u_sol} = {rhs}

Parameters:
  • mat (array) – Array with the system of equations. It can be stored in dense or sparse scheme.
  • rhs (array) – Array with right-hand-side of the system of equations.
Returns:

u_sol – Solution of the system of equations.
Return type:

array

Element subroutines

Each UEL subroutine computes the local stiffness matrix for a given finite element.

New elements can be added by including additional subroutines.

uelutil.uel3ntrian(coord, enu, Emod)[source]

Triangular element with 3 nodes

Parameters:
  • coord (ndarray) – Coordinates for the nodes of the element (3, 2).
  • enu (float) – Poisson coefficient (-1, 0.5).
  • Emod (float) – Young modulus (>0).
Returns:

kl – Local stiffness matrix for the element (6, 6).
Return type:

ndarray

Examples

>>> coord = np.array([
...         [0, 0],
...         [1, 0],
...         [0, 1]])
>>> stiff = uel3ntrian(coord, 1/3, 8/3)
>>> stiff_ex = 1/2 * np.array([
...            [4, 2, -3, -1, -1, -1],
...            [2, 4, -1, -1, -1, -3],
...            [-3, -1, 3, 0, 0, 1],
...            [-1, -1, 0, 1, 1, 0],
...            [-1, -1, 0, 1, 1, 0],
...            [-1, -3, 1, 0, 0, 3]])
>>> np.allclose(stiff, stiff_ex)
True
uelutil.uel4nquad(coord, enu, Emod)[source]

Quadrilateral element with 4 nodes

Parameters:
  • coord (ndarray) – Coordinates for the nodes of the element (4, 2).
  • enu (float) – Poisson coefficient (-1, 0.5).
  • Emod (float) – Young modulus (>0).
Returns:

kl – Local stiffness matrix for the element (8, 8).
Return type:

ndarray

Examples

>>> coord = np.array([[-1, -1], [1, -1], [1, 1], [-1, 1]])
>>> stiff = uel4nquad(coord, 1/3, 8/3)
>>> stiff_ex = 1/6 * np.array([
...             [ 8,  3, -5,  0, -4, -3,  1,  0],
...             [ 3,  8,  0,  1, -3, -4,  0, -5],
...             [-5,  0,  8, -3,  1,  0, -4,  3],
...             [ 0,  1, -3,  8,  0, -5,  3, -4],
...             [-4, -3,  1,  0,  8,  3, -5,  0],
...             [-3, -4,  0, -5,  3,  8,  0,  1],
...             [ 1,  0, -4,  3, -5,  0,  8, -3],
...             [ 0, -5,  3, -4,  0,  1, -3,  8]])
>>> np.allclose(stiff, stiff_ex)
True
uelutil.uel6ntrian(coord, enu, Emod)[source]

Triangular element with 6 nodes

Parameters:
  • coord (ndarray) – Coordinates for the nodes of the element (6, 2).
  • enu (float) – Poisson coefficient (-1, 0.5).
  • Emod (float) – Young modulus (>0).
Returns:

kl – Local stiffness matrix for the element (12, 12).
Return type:

ndarray

Examples

>>> coord = np.array([
...         [0, 0],
...         [1, 0],
...         [0, 1],
...         [0.5, 0],
...         [0.5, 0.5],
...         [0, 0.5]])
>>> stiff = uel6ntrian(coord,1/3, 8/3)
>>> stiff_ex = 1/6 * np.array([
...            [12, 6, 3, 1, 1, 1, -12, -4, 0, 0, -4, -4],
...            [6, 12, 1, 1, 1, 3, -4, -4, 0, 0, -4, -12],
...            [3, 1, 9, 0, 0, -1, -12, -4, 0, 4, 0, 0],
...            [1, 1, 0, 3, -1, 0, -4, -4, 4, 0, 0, 0],
...            [1, 1, 0, -1, 3, 0, 0, 0, 0, 4, -4, -4],
...            [1, 3, -1, 0, 0, 9, 0, 0, 4, 0, -4, -12],
...            [-12, -4, -12, -4, 0, 0, 32, 8, -8, -8, 0, 8],
...            [-4, -4, -4, -4, 0, 0, 8, 32, -8, -24, 8, 0],
...            [0, 0, 0, 4, 0, 4, -8, -8, 32, 8, -24, -8],
...            [0, 0, 4, 0, 4, 0, -8, -24, 8, 32, -8, -8],
...            [-4, -4, 0, 0, -4, -4, 0, 8, -24, -8, 32, 8],
...            [-4, -12, 0, 0, -4, -12, 8, 0, -8, -8, 8, 32]])
>>> np.allclose(stiff, stiff_ex)
True
uelutil.uelbeam2DU(coord, I, Emod)[source]
2D-2-noded beam element
without axial deformation
Parameters:
  • coord (ndarray) – Coordinates for the nodes of the element (2, 2).
  • A (float) – Cross section area.
  • Emod (float) – Young modulus (>0).
Returns:

kl – Local stiffness matrix for the element (4, 4).
Return type:

ndarray
uelutil.uelspring(coord, enu, Emod)[source]

1D-2-noded Spring element

Parameters:
  • coord (ndarray) – Coordinates for the nodes of the element (2, 2).
  • enu (float) – Fictitious parameter.
  • Emod (float) – Stiffness coefficient (>0).
Returns:

kl – Local stiffness matrix for the element (4, 4).
Return type:

ndarray

Examples

>>> coord = np.array([
...         [0, 0],
...         [1, 0]])
>>> stiff = uelspring(coord, 1/3, 8/3)
>>> stiff_ex = 8/3 * np.array([
...    [1, 0, -1, 0],
...    [0, 0, 0, 0],
...    [-1, 0, 1, 0],
...    [0, 0, 0, 0]])
>>> np.allclose(stiff, stiff_ex)
True
uelutil.ueltruss2D(coord, A, Emod)[source]

2D-2-noded truss element

Parameters:
  • coord (ndarray) – Coordinates for the nodes of the element (2, 2).
  • A (float) – Cross section area.
  • Emod (float) – Young modulus (>0).
Returns:

kl – Local stiffness matrix for the element (4, 4).
Return type:

ndarray

Examples

>>> coord = np.array([
...         [0, 0],
...         [1, 0]])
>>> stiff = ueltruss2D(coord, 1.0 , 1.0)
>>> stiff_ex =  np.array([
...    [1, 0, -1, 0],
...    [0, 0, 0, 0],
...    [-1, 0, 1, 0],
...    [0, 0, 0, 0]])
>>> np.allclose(stiff, stiff_ex)
True