Repository: AppliedMechanics-EAFIT/SolidsPy (language: python)

Issue:
fix numpy depracated

fix shape_tri6 function femutil

add github-action to CI/CD PyPi



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml
new file mode 100644
index 0000000..f4c4273
--- /dev/null
+++ b/.github/workflows/python-publish.yml
@@ -0,0 +1,153 @@
+name: SolidsPy CI/CD Pipeline
+
+on:
+  push:
+    branches:
+      - main
+      - develop_v2.0
+    tags:
+      - 'v*.*.*'
+  pull_request:
+    branches:
+      - main
+      - develop_v2.0
+
+concurrency:
+  group: ci-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  test:
+    name: Test on Python ${{ matrix.python-version }}
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: ['3.7', '3.8', '3.9', '3.10']
+
+    steps:
+      - name: Checkout code
+        uses: actions/checkout@v4
+
+      - name: Set up Python ${{ matrix.python-version }}
+        uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install build
+          pip install -e .[dev]
+
+      - name: Run tests with coverage
+        run: |
+          pytest --cov=solidspy tests/
+
+      - name: Generate coverage report
+        run: |
+          coverage xml
+
+      - name: Upload coverage report
+        uses: actions/upload-artifact@v4
+        with:
+          name: coverage-report-${{ matrix.python-version }}
+          path: coverage.xml
+
+  release:
+    name: Release Process
+    needs: test
+    runs-on: ubuntu-latest
+    if: startsWith(github.ref, 'refs/tags/v')
+
+    steps:
+      - name: Checkout code
+        uses: actions/checkout@v4
+
+      - name: Extract version from tag
+        run: echo "VERSION=${GITHUB_REF#refs/tags/v}" >> $GITHUB_ENV
+
+      - name: Validate VERSION Format
+        run: |
+          if [[ ! "${{ env.VERSION }}" =~ ^[0-9]+\.[0-9]+\.[0-9]+$ ]]; then
+            echo "Invalid version format: ${{ env.VERSION }}"
+            exit 1
+          fi
+
+      - name: Update __version__ in __init__.py
+        run: |
+          sed -i 's|__version__ = ".*"|__version__ = "${{ env.VERSION }}"|' src/solidspy/__init__.py
+
+      - name: Update version in pyproject.toml
+        run: |
+          sed -i 's/version = ".*"/version = "${{ env.VERSION }}"/' pyproject.toml
+
+      - name: Update version in conda-recipe/meta.yaml
+        run: |
+          sed -i 's|{% set version = ".*" %}|{% set version = "'${{ env.VERSION }}'" %}|' conda-recipe/meta.yaml
+
+      - name: Update release in docs/conf.py
+        run: |
+          sed -i 's|release = ".*"|release = "${{ env.VERSION }}"|' docs/conf.py
+
+      - name: Configure Git
+        run: |
+          git config user.name "github-actions[bot]"
+          git config user.email "github-actions[bot]@users.noreply.github.com"
+
+      - name: Commit version update
+        run: |
+          git add .
+          git commit -m "Update version to ${{ env.VERSION }}" || echo "No changes to commit"
+
+      - name: Push changes (Force)
+        uses: ad-m/github-push-action@v0.8.0
+        with:
+          github_token: ${{ secrets.GITHUB_TOKEN }}
+          branch: develop_v2.0
+          directory: .
+          force: true
+
+      # Build and Publish Steps
+      - name: Set up Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: '3.8'
+
+      - name: Install build tools
+        run: |
+          python -m pip install --upgrade pip
+          pip install build twine
+
+      - name: Build Python package
+        run: python -m build
+
+      - name: Publish package to TestPyPI
+        env:
+          TWINE_USERNAME: '__token__'
+          TWINE_PASSWORD: ${{ secrets.TEST_PYPI_API_TOKEN }}
+        run: |
+          python -m twine upload --repository-url https://test.pypi.org/legacy/ dist/* --verbose
+
+      - name: Set up Miniconda
+        uses: conda-incubator/setup-miniconda@v2
+        with:
+          auto-update-conda: true
+          python-version: '3.10'
+          channels: conda-forge
+          use-mamba: false
+
+      - name: Install conda-build and anaconda-client via conda
+        run: |
+          conda install -y conda-build anaconda-client
+          conda list anaconda-client
+
+      - name: Build Conda package
+        run: conda build conda-recipe --output-folder dist -c conda-forge
+
+      - name: Upload Conda package
+        env:
+          ANACONDA_API_TOKEN: ${{ secrets.ANACONDA_API_TOKEN }}
+        run: |
+          for pkg in dist/*/*.tar.bz2; do
+            conda run -n base anaconda upload "$pkg" --force
+          done
\ No newline at end of file
diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml
new file mode 100644
index 0000000..b3f7385
--- /dev/null
+++ b/conda-recipe/meta.yaml
@@ -0,0 +1,47 @@
+{% set name = "solidspy" %}
+{% set version = "1.4.9" %}
+
+package:
+  name: {{ name|lower }}
+  version: {{ version }}
+
+source:
+  git_url: https://github.com/kssgarcia/SolidsPy.git
+  git_tag: v{{ version }}
+
+build:
+  number: 0
+  script: python -m pip install . --no-deps -vv
+
+requirements:
+  host:
+    - python
+    - pip
+  run:
+    - python
+    - numpy
+    - scipy
+    - matplotlib
+    - easygui
+    - meshio >=5.0
+
+test:
+  imports:
+    - solidspy
+
+about:
+  home: "https://github.com/yourusername/SolidsPyKevin"
+  license: "MIT"
+  license_file: "LICENSE"
+  summary: "A Python package for finite element analysis."
+  description: |
+    SolidsPyKevin is a Python library for performing finite element analysis.
+  doc_url: "https://github.com/yourusername/SolidsPyKevin#readme"
+  dev_url: "https://github.com/yourusername/SolidsPyKevin"
+
+extra:
+  recipe-maintainers:
+    - kssgarcia
+
+channels:
+  - conda-forge
\ No newline at end of file
diff --git a/docs/conf.py b/docs/conf.py
index 5e70db7..d3ba85d 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -21,7 +21,7 @@
 import sys
 sys.path.insert(0, os.path.abspath('../solidspy'))
 
-import solidspy
+import src
 
 
 # -- General configuration ------------------------------------------------
@@ -65,7 +65,7 @@
 # The short X.Y version.
 version = '1.0'
 # The full version, including alpha/beta/rc tags.
-release = '1.0.16'
+release = "1.4.9",
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
diff --git a/docs/tutorials/old_format.rst b/docs/tutorials/old_format.rst
index 142da5b..28fc30b 100644
--- a/docs/tutorials/old_format.rst
+++ b/docs/tutorials/old_format.rst
@@ -4,7 +4,7 @@ Using the old format format for SolidsPy
 :Author: Nicolás Guarín-Zapata
 :Date: September 2019
 
-Until version 1.0.16 the data for a simulation was stored in 4 files
+Until version 1.2.8 the data for a simulation was stored in 4 files
 
 - ``nodes.txt``
 
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..5208340
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,55 @@
+[build-system]
+requires = ["setuptools >= 68.0.0", "wheel"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "SolidsPyKevin"
+version = "1.4.9"
+description = "A simple Finite Element program"
+readme = {file = "README.rst", content-type = "text/x-rst"}
+license = {file = "LICENSE"}
+authors = [
+    {name = "Juan Gomez", email = "jgomezc1@eafit.edu.co"},
+    {name = "Nicolas Guarin-Zapata", email = "nicoguarin@gmail.com"}
+]
+maintainers = [
+    {name = "Juan Gomez", email = "jgomezc1@eafit.edu.co"}
+]
+requires-python = ">=3.7"
+classifiers = [
+    "Development Status :: 4 - Beta",
+    "Intended Audience :: Education",
+    "Topic :: Scientific/Engineering",
+    "License :: OSI Approved :: MIT License",
+    "Programming Language :: Python :: 3",
+    "Programming Language :: Python :: 3.7",
+    "Programming Language :: Python :: 3.8",
+    "Programming Language :: Python :: 3.9",
+    "Programming Language :: Python :: 3.10"
+]
+keywords = ["finite-elements", "fem", "scientific-computing"]
+
+dependencies = [
+    "numpy",
+    "scipy",
+    "matplotlib",
+    "easygui",
+    "meshio==3.0"
+]
+
+[project.optional-dependencies]
+dev = [
+    "pytest",
+    "coverage",
+    "pytest-cov",
+    "black",
+    "flake8"
+]
+
+[project.urls]
+Homepage = "https://github.com/AppliedMechanics-EAFIT/SolidsPy"
+Repository = "https://github.com/AppliedMechanics-EAFIT/SolidsPy"
+BugTracker = "https://github.com/AppliedMechanics-EAFIT/SolidsPy/issues"
+
+[tool.setuptools.packages.find]
+where = ["src"]
diff --git a/requirements-dev.txt b/requirements-dev.txt
deleted file mode 100644
index 6d455de..0000000
--- a/requirements-dev.txt
+++ /dev/null
@@ -1,3 +0,0 @@
--r requirements.txt
-pytest
-coverage
diff --git a/requirements-docs.txt b/requirements-docs.txt
deleted file mode 100644
index efcd1c2..0000000
--- a/requirements-docs.txt
+++ /dev/null
@@ -1,2 +0,0 @@
--r requirements.txt
-
diff --git a/requirements.txt b/requirements.txt
deleted file mode 100644
index f8980ef..0000000
--- a/requirements.txt
+++ /dev/null
@@ -1,5 +0,0 @@
-numpy
-scipy
-matplotlib
-easygui
-meshio==3.0
diff --git a/setup.py b/setup.py
deleted file mode 100755
index 9587ea8..0000000
--- a/setup.py
+++ /dev/null
@@ -1,87 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-"""A simple Finite Element program
-
-https://github.com/AppliedMechanics-EAFIT/SolidsPy
-"""
-
-# Always prefer setuptools over distutils
-from setuptools import setup, find_packages
-# To use a consistent encoding
-from codecs import open
-from os import path
-
-here = path.abspath(path.dirname(__file__))
-
-# Get the long description from the README file
-with open(path.join(here, 'README.rst'), encoding='utf-8') as f:
-    long_description = f.read()
-
-requirements = ['numpy',
-                'scipy',
-                'matplotlib',
-                'easygui',
-                'meshio==3.0']
-
-setup(
-    name='solidspy',
-
-    version='1.0.16',
-
-    description='A simple Finite Element program',
-    long_description=long_description,
-    long_description_content_type="text/x-rst",
-
-    # The project's main homepage.
-    url='https://github.com/AppliedMechanics-EAFIT/SolidsPy',
-
-
-    # Author details
-    author='Juan Gomez <jgomezc1@eafit.edu.co>, Nicolas Guarin-Zapata <nicoguarin@gmail.com>',
-    author_email='jgomezc1@eafit.edu.co',
-
-    license='MIT',
-
-    # See https://pypi.python.org/pypi?%3Aaction=list_classifiers
-    classifiers=[
-        # How mature is this project? Common values are
-        #   3 - Alpha
-        #   4 - Beta
-        #   5 - Production/Stable
-        'Development Status :: 4 - Beta',
-
-        # Indicate who your project is intended for
-        'Intended Audience :: Education',
-        'Topic :: Scientific/Engineering',
-
-        # Pick your license as you wish (should match "license" above)
-        'License :: OSI Approved :: MIT License',
-
-        # Specify the Python versions you support here. In particular, ensure
-        # that you indicate whether you support Python 2, Python 3 or both.
-        'Programming Language :: Python :: 3',
-    ],
-
-    keywords='finite-elements fem scientific-computing',
-
-
-    # You can just specify the packages manually here if your project is
-    # simple. Or you can use find_packages().
-    packages=find_packages(exclude=['contrib', 'docs', 'tests']),
-
-    # List run-time dependencies here.  These will be installed by pip when
-    # your project is installed. For an analysis of "install_requires" vs pip's
-    # requirements files see:
-    # https://packaging.python.org/en/latest/requirements.html
-    install_requires=requirements,
-
-
-    # To provide executable scripts, use entry points in preference to the
-    # "scripts" keyword. Entry points provide cross-platform support and allow
-    # pip to create the appropriate form of executable for the target platform.
-    entry_points={
-        'console_scripts': [
-            'sample=sample:main',
-        ],
-    },
-)
diff --git a/solidspy/__init__.py b/src/solidspy/__init__.py
similarity index 97%
rename from solidspy/__init__.py
rename to src/solidspy/__init__.py
index 51e9e5b..50f10c9 100755
--- a/solidspy/__init__.py
+++ b/src/solidspy/__init__.py
@@ -11,7 +11,7 @@
            "uelutil",
            "solids_GUI"]
 
-__version__ = "1.0.16"
+__version__ = "1.4.9"
 
 
 __citation__ = """@software{solidspy,
diff --git a/solidspy/__main__.py b/src/solidspy/__main__.py
similarity index 100%
rename from solidspy/__main__.py
rename to src/solidspy/__main__.py
diff --git a/solidspy/assemutil.py b/src/solidspy/assemutil.py
similarity index 72%
rename from solidspy/assemutil.py
rename to src/solidspy/assemutil.py
index a56b4bd..d9a92a9 100755
--- a/solidspy/assemutil.py
+++ b/src/solidspy/assemutil.py
@@ -8,12 +8,16 @@
 
 """
 import numpy as np
-from scipy.sparse import coo_matrix
+from numpy.typing import NDArray
+from scipy.sparse import coo_matrix, csr_matrix, spmatrix
+from typing import Callable, Optional, Tuple
 import solidspy.uelutil as ue
 import solidspy.femutil as fem
 
 
-def eqcounter(cons, ndof_node=2):
+def eqcounter(
+    cons: NDArray[np.int_], ndof_node: int = 2
+) -> Tuple[int, NDArray[np.int_]]:
     """Count active equations
 
     Creates boundary conditions array bc_array
@@ -23,6 +27,9 @@ def eqcounter(cons, ndof_node=2):
     cons : ndarray.
       Array with constraints for each node.
 
+    ndof_node : int, optional
+      Number of degrees of freedom per node. Default is 2.
+
     Returns
     -------
     neq : int
@@ -44,7 +51,13 @@ def eqcounter(cons, ndof_node=2):
     return neq, bc_array
 
 
-def DME(cons, elements, ndof_node=2, ndof_el_max=18, ndof_el=None):
+def DME(
+    cons: NDArray[np.int_],
+    elements: NDArray[np.int_],
+    ndof_node: int = 2,
+    ndof_el_max: int = 18,
+    ndof_el: Optional[Callable[[int], int]] = None,
+) -> Tuple[NDArray[np.int_], NDArray[np.int_], int]:
     """Create assembly array operator
 
     Count active equations, create boundary conditions array ``bc_array``
@@ -62,7 +75,7 @@ def DME(cons, elements, ndof_node=2, ndof_el_max=18, ndof_el=None):
       Number of maximum degrees of freedom per element. By default it is
       18.
     ndof_el : callable, optional
-      Function that return number of degrees of freedom for elements. It
+      Function that returns the number of degrees of freedom for elements. It
       is needed for user elements.
 
     Returns
@@ -76,7 +89,7 @@ def DME(cons, elements, ndof_node=2, ndof_el_max=18, ndof_el=None):
 
     """
     nels = elements.shape[0]
-    assem_op = np.zeros([nels, ndof_el_max], dtype=np.integer)
+    assem_op = np.zeros([nels, ndof_el_max], dtype=int)
     neq, bc_array = eqcounter(cons, ndof_node=ndof_node)
     for ele in range(nels):
         iet = elements[ele, 1]
@@ -88,7 +101,7 @@ def DME(cons, elements, ndof_node=2, ndof_el_max=18, ndof_el=None):
     return assem_op, bc_array, neq
 
 
-def ele_fun(eletype):
+def ele_fun(eletype: int) -> Callable[[NDArray[np.float64], NDArray[np.float64]], Tuple[NDArray[np.float64], NDArray[np.float64]]]:
     """Return the function for the element type given
 
     Parameters
@@ -109,14 +122,23 @@ def ele_fun(eletype):
         5: ue.spring,
         6: ue.truss2D,
         7: ue.beam2DU,
-        8: ue.beam2D}
+        8: ue.beam2D,
+    }
     try:
         return elem_id[eletype]
-    except:
+    except KeyError:
         raise ValueError("You entered an invalid type of element.")
 
 
-def retriever(elements, mats, nodes, ele, uel=None):
+def retriever(
+    elements: NDArray[np.int_],
+    mats: NDArray[np.float64],
+    nodes: NDArray[np.float64],
+    ele: int,
+    uel: Optional[
+        Callable[[NDArray[np.float64], NDArray[np.float64]], Tuple[NDArray[np.float64], NDArray[np.float64]]]
+    ] = None,
+) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """Computes the elemental stiffness matrix of element ``ele``
 
     Parameters
@@ -129,6 +151,8 @@ def retriever(elements, mats, nodes, ele, uel=None):
       Array with the nodal numbers and coordinates.
     ele : int.
       Identifier of the element to be assembled.
+    uel : callable, optional
+      User element function that returns stiffness and mass matrices.
 
     Returns
     -------
@@ -146,7 +170,17 @@ def retriever(elements, mats, nodes, ele, uel=None):
     return kloc, mloc
 
 
-def assembler(elements, mats, nodes, neq, assem_op, sparse=True, uel=None):
+def assembler(
+    elements: NDArray[np.int_],
+    mats: NDArray[np.float64],
+    nodes: NDArray[np.float64],
+    neq: int,
+    assem_op: NDArray[np.int_],
+    sparse: bool = True,
+    uel: Optional[
+        Callable[[NDArray[np.float64], NDArray[np.float64]], Tuple[NDArray[np.float64], NDArray[np.float64]]]
+    ] = None,
+) -> Tuple[spmatrix, spmatrix]:
     """Assembles the global stiffness matrix
 
     Parameters
@@ -169,12 +203,10 @@ def assembler(elements, mats, nodes, neq, assem_op, sparse=True, uel=None):
 
     Returns
     -------
-    kglob : ndarray (float)
-      Array with the global stiffness matrix. It might be
-      dense or sparse, depending on the value of _sparse_
-    mglob : ndarray (float)
-      Array with the global mass matrix. It might be
-      dense or sparse, depending on the value of _sparse_
+    kglob : sparse matrix (float)
+      Array with the global stiffness matrix. It is sparse in CSR format.
+    mglob : sparse matrix (float)
+      Array with the global mass matrix. It is sparse in CSR format.
 
     """
     if sparse:
@@ -187,7 +219,16 @@ def assembler(elements, mats, nodes, neq, assem_op, sparse=True, uel=None):
     return kglob, mglob
 
 
-def dense_assem(elements, mats, nodes, neq, assem_op, uel=None):
+def dense_assem(
+    elements: NDArray[np.int_],
+    mats: NDArray[np.float64],
+    nodes: NDArray[np.float64],
+    neq: int,
+    assem_op: NDArray[np.int_],
+    uel: Optional[
+        Callable[[NDArray[np.float64], NDArray[np.float64]], Tuple[NDArray[np.float64], NDArray[np.float64]]]
+    ] = None,
+) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """
     Assembles the global stiffness matrix
     using a dense storing scheme
@@ -204,7 +245,7 @@ def dense_assem(elements, mats, nodes, neq, assem_op, uel=None):
       Assembly operator.
     neq : int
       Number of active equations in the system.
-    uel : callable function (optional)
+    uel : callable, optional
       Python function that returns the local stiffness matrix.
 
     Returns
@@ -212,6 +253,9 @@ def dense_assem(elements, mats, nodes, neq, assem_op, uel=None):
     kglob : ndarray (float)
       Array with the global stiffness matrix in a dense numpy
       array.
+    mglob : ndarray (float)
+      Array with the global mass matrix in a dense numpy
+      array.
 
     """
     kglob = np.zeros((neq, neq))
@@ -233,13 +277,22 @@ def dense_assem(elements, mats, nodes, neq, assem_op, uel=None):
     return kglob, mglob
 
 
-def sparse_assem(elements, mats, nodes, neq, assem_op, uel=None):
+def sparse_assem(
+    elements: NDArray[np.int_],
+    mats: NDArray[np.float64],
+    nodes: NDArray[np.float64],
+    neq: int,
+    assem_op: NDArray[np.int_],
+    uel: Optional[
+        Callable[[NDArray[np.float64], NDArray[np.float64]], Tuple[NDArray[np.float64], NDArray[np.float64]]]
+    ] = None,
+) -> Tuple[csr_matrix, csr_matrix]:
     """
     Assembles the global stiffness matrix
     using a sparse storing scheme
 
     The scheme used to assemble is COOrdinate list (COO), and
-    it converted to Compressed Sparse Row (CSR) afterward
+    it is converted to Compressed Sparse Row (CSR) afterward
     for the solution phase [1]_.
 
     Parameters
@@ -254,7 +307,7 @@ def sparse_assem(elements, mats, nodes, neq, assem_op, uel=None):
       Assembly operator.
     neq : int
       Number of active equations in the system.
-    uel : callable function (optional)
+    uel : callable, optional
       Python function that returns the local stiffness matrix.
 
     Returns
@@ -262,6 +315,9 @@ def sparse_assem(elements, mats, nodes, neq, assem_op, uel=None):
     kglob : sparse matrix (float)
       Array with the global stiffness matrix in a sparse
       Compressed Sparse Row (CSR) format.
+    mglob : sparse matrix (float)
+      Array with the global mass matrix in a sparse
+      Compressed Sparse Row (CSR) format.
 
     References
     ----------
@@ -295,7 +351,12 @@ def sparse_assem(elements, mats, nodes, neq, assem_op, uel=None):
     return stiff, mass
 
 
-def loadasem(loads, bc_array, neq, ndof_node=2):
+def loadasem(
+    loads: NDArray[np.float64],
+    bc_array: NDArray[np.int_],
+    neq: int,
+    ndof_node: int = 2,
+) -> NDArray[np.float64]:
     """Assembles the global Right Hand Side Vector
 
     Parameters
@@ -307,6 +368,8 @@ def loadasem(loads, bc_array, neq, ndof_node=2):
     neq : int
       Number of equations in the system after removing the nodes
       with imposed displacements.
+    ndof_node : int, optional
+      Number of degrees of freedom per node. By default it is 2.
 
     Returns
     -------
diff --git a/solidspy/blochutil.py b/src/solidspy/blochutil.py
similarity index 86%
rename from solidspy/blochutil.py
rename to src/solidspy/blochutil.py
index a05f591..a853f6a 100644
--- a/solidspy/blochutil.py
+++ b/src/solidspy/blochutil.py
@@ -3,27 +3,31 @@
 """
 Bloch analysis utilities
 ------------------------
-Routines used in Bloch Analisis.
+Routines used in Bloch Analysis.
 
 """
 import numpy as np
-from scipy.sparse import coo_matrix
+from scipy.sparse import coo_matrix, csr_matrix
+from typing import List, Dict
 
 
-#%% Auxiliar
-def reassign_dof_bloch(ndofs, dofs_ima, dofs_ref):
+#%% Auxiliary
+def reassign_dof_bloch(
+    ndofs: int,
+    dofs_ima: List[int],
+    dofs_ref: List[int]
+) -> List[int]:
     """Find dof number after applying boundary conditions
 
     Parameters
     ----------
     ndofs : int
         Number of degrees of freedom.
-    dofs_ima : list, int
+    dofs_ima : list of int
         List with image degrees of freedom.
-    dofs_ref : list, int
+    dofs_ref : list of int
         List with reference degrees of freedom.
 
-
     Examples
     --------
     
@@ -41,7 +45,7 @@ def reassign_dof_bloch(ndofs, dofs_ima, dofs_ref):
 
     Returns
     -------
-    new_num : list, int
+    new_num : list of int
         List with new numbering for the dof after applying Bloch analysis.
 
     """
@@ -64,9 +68,18 @@ def reassign_dof_bloch(ndofs, dofs_ima, dofs_ref):
 
 
 #%% Sparse matrices
-def bloch_transform_mat(wavevector, coords, ndofs, nodes_ref, nodes_ima,
-                        dofs_ref, dofs_ima, nodes_ref_dof, nodes_ima_dof,
-                        new_num):
+def bloch_transform_mat(
+    wavevector: np.ndarray,
+    coords: np.ndarray,
+    ndofs: int,
+    nodes_ref: List[int],
+    nodes_ima: List[int],
+    dofs_ref: List[int],
+    dofs_ima: List[int],
+    nodes_ref_dof: Dict[int, List[int]],
+    nodes_ima_dof: Dict[int, List[int]],
+    new_num: List[int]
+) -> csr_matrix:
     """Form transformation matrices for Bloch conditions
 
     Parameters
@@ -77,20 +90,19 @@ def bloch_transform_mat(wavevector, coords, ndofs, nodes_ref, nodes_ima,
         Coordinates array.
     ndofs : int
         Number of degrees of freedom.
-    nodes_ref : list, int
+    nodes_ref : list of int
         List with reference nodes.
-    nodes_ima : list, int
+    nodes_ima : list of int
         List with image nodes.
-    dofs_ref : list, int
+    dofs_ref : list of int
         List with reference degrees of freedom.        
-    dofs_ima : list, int
+    dofs_ima : list of int
         List with image degrees of freedom.
-    nodes_ref_dof : dictionary
+    nodes_ref_dof : dict of int to list of int
         Dictionary that maps nodes to dof for reference nodes.
-    nodes_ima_dof : dictionary
+    nodes_ima_dof : dict of int to list of int
         Dictionary that maps nodes to dof for image nodes.
 
-
     Examples
     --------
     
diff --git a/solidspy/femutil.py b/src/solidspy/femutil.py
similarity index 71%
rename from solidspy/femutil.py
rename to src/solidspy/femutil.py
index d13dd1b..91b0066 100755
--- a/solidspy/femutil.py
+++ b/src/solidspy/femutil.py
@@ -1,11 +1,11 @@
- # -*- coding: utf-8 -*-
+# -*- coding: utf-8 -*-
 """
 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.
@@ -13,7 +13,6 @@
 
 The notation used is similar to the one used by Bathe [1]_.
 
-
 References
 ----------
 .. [1] Bathe, Klaus-Jürgen. Finite element procedures. Prentice Hall,
@@ -22,17 +21,19 @@
 """
 import numpy as np
 import solidspy.gaussutil as gau
+from typing import Tuple, Callable
+from numpy.typing import NDArray
 
 
-def eletype(eletype):
+def eletype(eletype: int) -> Tuple[int, int, int]:
     """Assigns number to degrees of freedom
 
-    According to iet assigns number of degrees of freedom, number of
+    According to `eletype` assigns number of degrees of freedom, number of
     nodes and minimum required number of integration points.
 
     Parameters
     ----------
-    eletype :  int
+    eletype : int
       Type of element. These are:
         1. 4 node bilinear quadrilateral.
         2. 6 node quadratic triangle.
@@ -51,6 +52,11 @@ def eletype(eletype):
     ngpts : int
       Number of Gauss points for the selected element.
 
+    Raises
+    ------
+    ValueError
+        If an invalid element type is provided.
+
     """
     elem_id = {
         1: (8, 4, 4),
@@ -63,14 +69,14 @@ def eletype(eletype):
         8: (6, 2, 3)}
     try:
         return elem_id[eletype]
-    except:
+    except KeyError:
         raise ValueError("You entered an invalid type of element.")
 
 
 #%% Shape functions and derivatives
 
 # Triangles
-def shape_tri3(r, s):
+def shape_tri3(r: float, s: float) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """
     Shape functions and derivatives for a linear element
 
@@ -83,9 +89,9 @@ def shape_tri3(r, s):
 
     Returns
     -------
-    N : ndarray (float)
+    N : NDArray[np.float64]
         Array with the shape functions evaluated at the point (r, s).
-    dNdr : ndarray (float)
+    dNdr : NDArray[np.float64]
         Array with the derivative of the shape functions evaluated at
         the point (r, s).
 
@@ -96,17 +102,19 @@ def shape_tri3(r, s):
 
     >>> N, _ = shape_tri3(0, 0)
     >>> N_ex = np.array([
-    ...    [1, 0, 0, 0, 0, 0],
-    ...    [0, 1, 0, 0, 0, 0]])
+    ...    [1, 0, 0],
+    ...    [0, 1, 0],
+    ...    [0, 0, 1]])
     >>> np.allclose(N, N_ex)
     True
 
     and
 
-    >>> N, _ = shape_tri3(1/2, 1/2)
+    >>> N, _ = shape_tri3(0.5, 0.5)
     >>> N_ex = np.array([
-    ...    [0, 0, 1/2, 0, 1/2, 0],
-    ...    [0, 0, 0, 1/2, 0, 1/2]])
+    ...    [0, 0.5, 0.5],
+    ...    [0, 1, 0],
+    ...    [0, 0, 1]])
     >>> np.allclose(N, N_ex)
     True
 
@@ -118,9 +126,9 @@ def shape_tri3(r, s):
     return N, dNdr
 
 
-def shape_tri6(r, s):
+def shape_tri6(r: float, s: float) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """
-    Shape functions and derivatives for a quadratic element
+    Shape functions and derivatives for a quadratic triangular element.
 
     Parameters
     ----------
@@ -131,26 +139,63 @@ def shape_tri6(r, s):
 
     Returns
     -------
-    N : ndarray (float)
-        Array with the shape functions evaluated at the point (r, s).
-    dNdr : ndarray (float)
-        Array with the derivative of the shape functions evaluated at
-        the point (r, s).
+    N : NDArray[np.float64]
+        Shape functions evaluated at the point (r, s).
+    dNdrds : NDArray[np.float64]
+        Derivatives of the shape functions with respect to r and s.
+
     """
-    N = np.array(
-        [(1 - r - s) - 2*r*(1 - r - s) - 2*s*(1 - r - s),
-         r - 2*r*(1 - r - s) - 2*r*s,
-         s - 2*r*s - 2*s*(1-r-s),
-         4*r*(1 - r - s),
-         4*r*s,
-         4*s*(1 - r - s)])
-    dNdr = np.array([
-        [4*r + 4*s - 3, 4*r - 1, 0, -8*r - 4*s + 4, 4*s, -4*s],
-        [4*r + 4*s - 3, 0, 4*s - 1, -4*r, 4*r, -4*r - 8*s + 4]])
-    return N, dNdr
+    # Convert inputs to NumPy arrays
+    r = np.asarray(r)
+    s = np.asarray(s)
+    
+    # Check if inputs are scalar
+    scalar_input = False
+    if r.ndim == 0 and s.ndim == 0:
+        r = r.reshape(1)
+        s = s.reshape(1)
+        scalar_input = True
+
+    # Compute shape functions N1 to N6
+    N = np.vstack((
+        (1 - r - s) - 2 * r * (1 - r - s) - 2 * s * (1 - r - s),
+        r - 2 * r * (1 - r - s) - 2 * r * s,
+        s - 2 * r * s - 2 * s * (1 - r - s),
+        4 * r * (1 - r - s),
+        4 * r * s,
+        4 * s * (1 - r - s)
+    )).T  # Shape: (n_points, 6)
+
+    # Compute derivatives with respect to r
+    dNdr = np.vstack((
+        4 * r + 4 * s - 3,
+        4 * r - 1,
+        np.zeros_like(r),
+        -8 * r - 4 * s + 4,
+        4 * s,
+        -4 * s
+    )).T  # Shape: (n_points, 6)
+
+    # Compute derivatives with respect to s
+    dNds = np.vstack((
+        4 * r + 4 * s - 3,
+        np.zeros_like(r),
+        4 * s - 1,
+        -4 * r,
+        4 * r,
+        -4 * r - 8 * s + 4
+    )).T  # Shape: (n_points, 6)
+
+    # Stack derivatives into a single array
+    dNdrds = np.stack((dNdr, dNds), axis=1)  # Shape: (n_points, 2, 6)
+
+    if scalar_input:
+        return N[0], dNdrds[0]
+    else:
+        return N, dNdrds
 
 # Quadrilaterals
-def shape_quad4(r, s):
+def shape_quad4(r: float, s: float) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """
     Shape functions and derivatives for a bilinear element
 
@@ -163,9 +208,9 @@ def shape_quad4(r, s):
 
     Returns
     -------
-    N : ndarray (float)
+    N : NDArray[np.float64]
         Array with the shape functions evaluated at the point (r, s).
-    dNdr : ndarray (float)
+    dNdr : NDArray[np.float64]
         Array with the derivative of the shape functions evaluated at
         the point (r, s).
 
@@ -176,8 +221,8 @@ def shape_quad4(r, s):
 
     >>> N, _ = shape_quad4(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]])
+    ...    [0.25, 0.0, 0.25, 0.0, 0.25, 0.0, 0.25, 0.0],
+    ...    [0.0, 0.25, 0.0, 0.25, 0.0, 0.25, 0.0, 0.25]])
     >>> np.allclose(N, N_ex)
     True
 
@@ -185,8 +230,8 @@ def shape_quad4(r, s):
 
     >>> N, _ = shape_quad4(1, 1)
     >>> N_ex = np.array([
-    ...    [0, 0, 0, 0, 1, 0, 0, 0],
-    ...    [0, 0, 0, 0, 0, 1, 0, 0]])
+    ...    [0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0],
+    ...    [0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0]])
     >>> np.allclose(N, N_ex)
     True
 
@@ -202,7 +247,7 @@ def shape_quad4(r, s):
     return N, dNdr
 
 
-def shape_quad9(r, s):
+def shape_quad9(r: float, s: float) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """
     Shape functions and derivatives for a biquadratic element
 
@@ -215,9 +260,9 @@ def shape_quad9(r, s):
 
     Returns
     -------
-    N : ndarray (float)
+    N : NDArray[np.float64]
         Array with the shape functions evaluated at the point (r, s).
-    dNdr : ndarray (float)
+    dNdr : NDArray[np.float64]
         Array with the derivative of the shape functions evaluated at
         the point (r, s).
     """
@@ -253,10 +298,9 @@ def shape_quad9(r, s):
     return N, dNdr
 
 
-
-def shape_quad8(r, s):
+def shape_quad8(r: float, s: float) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """
-    Shape functions and derivatives for a 8-noded serendipity element
+    Shape functions and derivatives for an 8-noded serendipity element
 
     Parameters
     ----------
@@ -267,12 +311,13 @@ def shape_quad8(r, s):
 
     Returns
     -------
-    N : ndarray (float)
+    N : NDArray[np.float64]
         Array with the shape functions evaluated at the point (r, s).
-    dNdr : ndarray (float)
+    dNdr : NDArray[np.float64]
         Array with the derivative of the shape functions evaluated at
         the point (r, s).
     """
+
     N = 0.25*np.array(
         [(1.0 - r)*(1.0 - s) - (1.0 - r)*(1.0 - s**2) - (1.0 - s)*(1.0 - r**2),
          (1.0 + r)*(1.0 - s) - (1.0 + r)*(1.0 - s**2) - (1.0 - s)*(1.0 - r**2),
@@ -302,8 +347,9 @@ def shape_quad8(r, s):
     return N, dNdr
 
 
-## 3D elements
-def shape_tet4(r, s, t):
+
+# 3D elements
+def shape_tet4(r: float, s: float, t: float) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """
     Shape functions and derivatives for a linear tetrahedron
 
@@ -318,9 +364,9 @@ def shape_tet4(r, s, t):
 
     Returns
     -------
-    N : ndarray (float)
+    N : NDArray[np.float64]
         Array with the shape functions evaluated at the point (r, s, t).
-    dNdr : ndarray (float)
+    dNdr : NDArray[np.float64]
         Array with the derivative of the shape functions evaluated at
         the point (r, s, t).
     """
@@ -332,7 +378,7 @@ def shape_tet4(r, s, t):
     return N, dNdr
 
 
-def shape_hex8(r, s, t):
+def shape_hex8(r: float, s: float, t: float) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """
     Shape functions and derivatives for a trilinear element
 
@@ -347,9 +393,9 @@ def shape_hex8(r, s, t):
 
     Returns
     -------
-    N : ndarray (float)
+    N : NDArray[np.float64]
         Array with the shape functions evaluated at the point (r, s, t).
-    dNdr : ndarray (float)
+    dNdr : NDArray[np.float64]
         Array with the derivative of the shape functions evaluated at
         the point (r, s, t).
     """
@@ -375,7 +421,12 @@ def shape_hex8(r, s, t):
 
 
 #%% Derivative matrices
-def elast_diff_2d(r, s, coord, element):
+def elast_diff_2d(
+    r: float,
+    s: float,
+    coord: NDArray[np.float64],
+    element: Callable[[float, float], Tuple[NDArray[np.float64], NDArray[np.float64]]]
+) -> Tuple[NDArray[np.float64], NDArray[np.float64], float]:
     """
     Interpolation matrices for elements for plane elasticity
 
@@ -385,19 +436,22 @@ def elast_diff_2d(r, s, coord, element):
         Horizontal coordinate of the evaluation point.
     s : float
         Vertical coordinate of the evaluation point.
-    coord : ndarray (float)
+    coord : NDArray[np.float64]
         Coordinates of the element.
+    element : Callable
+        Shape function generator for the element.
 
     Returns
     -------
-    H : ndarray (float)
+    H : NDArray[np.float64]
         Array with the shape functions evaluated at the point (r, s)
         for each degree of freedom.
-    B : ndarray (float)
+    B : NDArray[np.float64]
         Array with the displacement to strain matrix evaluated
         at the point (r, s).
     det : float
         Determinant of the Jacobian.
+
     """
     N, dNdr = element(r, s)
     det, jaco_inv = jacoper(dNdr, coord)
@@ -413,9 +467,14 @@ def elast_diff_2d(r, s, coord, element):
     return H, B, det
 
 
-def elast_diff_axi(r, s, coord, element):
+def elast_diff_axi(
+    r: float,
+    s: float,
+    coord: NDArray[np.float64],
+    element: Callable[[float, float], Tuple[NDArray[np.float64], NDArray[np.float64]]]
+) -> Tuple[NDArray[np.float64], NDArray[np.float64], float]:
     """
-    Interpolation matrices for elements for axisymetric elasticity
+    Interpolation matrices for elements for axisymmetric elasticity
 
     Parameters
     ----------
@@ -423,19 +482,27 @@ def elast_diff_axi(r, s, coord, element):
         Horizontal coordinate of the evaluation point.
     s : float
         Vertical coordinate of the evaluation point.
-    coord : ndarray (float)
+    coord : NDArray[np.float64]
         Coordinates of the element.
+    element : Callable
+        Shape function generator for the element.
 
     Returns
     -------
-    H : ndarray (float)
+    H : NDArray[np.float64]
         Array with the shape functions evaluated at the point (r, s)
         for each degree of freedom.
-    B : ndarray (float)
+    B : NDArray[np.float64]
         Array with the displacement to strain matrix evaluated
         at the point (r, s).
     det : float
         Determinant of the Jacobian.
+
+    Raises
+    ------
+    ValueError
+        If the radial coordinate `x` is negative.
+
     """
     N, dNdr = element(r, s)
     x = N.dot(coord[:, 0])
@@ -456,11 +523,15 @@ def elast_diff_axi(r, s, coord, element):
     return H, B, det
 
 
-## 3D elements
-def elast_diff_3d(r, s, t, coord, interp=shape_hex8):
+def elast_diff_3d(
+    r: float,
+    s: float,
+    t: float,
+    coord: NDArray[np.float64],
+    interp: Callable[[float, float, float], Tuple[NDArray[np.float64], NDArray[np.float64]]] = shape_hex8
+) -> Tuple[NDArray[np.float64], NDArray[np.float64], float]:
     """
-    Interpolation matrices for a trilinear element for
-    elasticity
+    Interpolation matrices for a trilinear element for elasticity
 
     Parameters
     ----------
@@ -470,19 +541,22 @@ def elast_diff_3d(r, s, t, coord, interp=shape_hex8):
         Horizontal coordinate of the evaluation point.
     t : float
         Vertical coordinate of the evaluation point.
-    coord : ndarray (float)
+    coord : NDArray[np.float64]
         Coordinates of the element.
+    interp : Callable, optional
+        Shape function generator for the element. Defaults to `shape_hex8`.
 
     Returns
     -------
-    H : ndarray (float)
+    H : NDArray[np.float64]
         Array with the shape functions evaluated at the point (r, s, t)
         for each degree of freedom.
-    B : ndarray (float)
+    B : NDArray[np.float64]
         Array with the displacement to strain matrix evaluated
         at the point (r, s, t).
     det : float
         Determinant of the Jacobian.
+
     """
     N, dNdr = interp(r, s, t)
     det, jaco_inv = jacoper(dNdr, coord)
@@ -505,23 +579,33 @@ def elast_diff_3d(r, s, t, coord, interp=shape_hex8):
 
 
 #%%
-def jacoper(dNdr, coord):
+def jacoper(
+    dNdr: NDArray[np.float64],
+    coord: NDArray[np.float64]
+) -> Tuple[float, NDArray[np.float64]]:
     """
     Compute the Jacobian of the transformation evaluated at
     the Gauss point
 
     Parameters
     ----------
-    dNdr : ndarray
+    dNdr : NDArray[np.float64]
       Derivatives of the interpolation function with respect to the
       natural coordinates.
-    coord : ndarray
+    coord : NDArray[np.float64]
       Coordinates of the nodes of the element (nnodes, ndim).
 
     Returns
     -------
-    jaco_inv : ndarray (ndim, ndim)
-      Jacobian of the transformation evaluated at a point.
+    det : float
+      Determinant of the Jacobian.
+    jaco_inv : NDArray[np.float64]
+      Inverse of the Jacobian matrix.
+
+    Raises
+    ------
+    ValueError
+        If the Jacobian determinant is close to zero or negative.
 
     """
     jaco = dNdr @ coord
@@ -537,14 +621,14 @@ def jacoper(dNdr, coord):
 
 
 #%% Material routines
-def umat(params):
-    """2D Elasticity consitutive matrix in plane stress
+def umat(params: Tuple[float, float]) -> NDArray[np.float64]:
+    """2D Elasticity constitutive matrix in plane stress
 
     For plane strain use effective properties.
 
     Parameters
     ----------
-    params : tuple
+    params : Tuple[float, float]
         Material parameters in the following order:
 
             E : float
@@ -554,13 +638,13 @@ def umat(params):
 
     Returns
     -------
-    C : ndarray
+    C : NDArray[np.float64]
       Constitutive tensor in Voigt notation.
 
     Examples
     --------
 
-    >>> params = 8/3, 1/3
+    >>> params = (8/3, 1/3)
     >>> C = umat(params)
     >>> C_ex = np.array([
     ...    [3, 1, 0],
@@ -583,14 +667,14 @@ def umat(params):
     return C
 
 
-def elast_mat_2d(params):
-    """2D Elasticity consitutive matrix in plain stress
+def elast_mat_2d(params: Tuple[float, float]) -> NDArray[np.float64]:
+    """2D Elasticity constitutive matrix in plane stress
 
-    For plane stress use effective properties.
+    For plane strain use effective properties.
 
     Parameters
     ----------
-    params : tuple
+    params : Tuple[float, float]
         Material parameters in the following order:
 
             E : float
@@ -600,7 +684,7 @@ def elast_mat_2d(params):
 
     Returns
     -------
-    C : ndarray
+    C : NDArray[np.float64]
       Constitutive tensor in Voigt notation.
 
     """
@@ -615,12 +699,12 @@ def elast_mat_2d(params):
     return C
 
 
-def elast_mat_axi(params):
-    """Elasticity consitutive matrix for axisymmetric problems.
+def elast_mat_axi(params: Tuple[float, float]) -> NDArray[np.float64]:
+    """Elasticity constitutive matrix for axisymmetric problems.
 
     Parameters
     ----------
-    params : tuple
+    params : Tuple[float, float]
         Material parameters in the following order:
 
             E : float
@@ -630,7 +714,7 @@ def elast_mat_axi(params):
 
     Returns
     -------
-    C : ndarray
+    C : NDArray[np.float64]
       Constitutive tensor in Voigt notation.
 
     """
@@ -645,12 +729,12 @@ def elast_mat_axi(params):
     return C
 
 
-def elast_mat(params):
-    """3D Elasticity consitutive matrix.
+def elast_mat(params: Tuple[float, float]) -> NDArray[np.float64]:
+    """3D Elasticity constitutive matrix.
 
     Parameters
     ----------
-    params : tuple
+    params : Tuple[float, float]
         Material parameters in the following order:
 
             E : float
@@ -660,7 +744,7 @@ def elast_mat(params):
 
     Returns
     -------
-    C : ndarray
+    C : NDArray[np.float64]
         Constitutive tensor in Voigt notation.
 
     """
@@ -677,25 +761,27 @@ def elast_mat(params):
     return C
 
 
-
 #%% Elemental strains
-def str_el3(coord, ul):
+def str_el3(
+    coord: NDArray[np.float64],
+    ul: NDArray[np.float64]
+) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """Compute the strains at each element integration point
 
     This one is used for 3-noded triangular elements.
 
     Parameters
     ----------
-    coord : ndarray
+    coord : NDArray[np.float64]
       Coordinates of the nodes of the element (nn, 2).
-    ul : ndarray
+    ul : NDArray[np.float64]
       Array with displacements for the element.
 
     Returns
     -------
-    epsGT : ndarray
+    epsGT : NDArray[np.float64]
       Strain components for the Gauss points.
-    xl : ndarray
+    xl : NDArray[np.float64]
       Configuration of the Gauss points after deformation.
 
     """
@@ -711,23 +797,26 @@ def str_el3(coord, ul):
     return epsG.T, xl
 
 
-def str_el6(coord, ul):
+def str_el6(
+    coord: NDArray[np.float64],
+    ul: NDArray[np.float64]
+) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """Compute the strains at each element integration point
 
     This one is used for 6-noded triangular elements.
 
     Parameters
     ----------
-    coord : ndarray
+    coord : NDArray[np.float64]
       Coordinates of the nodes of the element (6, 2).
-    ul : ndarray
+    ul : NDArray[np.float64]
       Array with displacements for the element.
 
     Returns
     -------
-    epsGT : ndarray
+    epsGT : NDArray[np.float64]
       Strain components for the Gauss points.
-    xl : ndarray
+    xl : NDArray[np.float64]
       Configuration of the Gauss points after deformation.
 
     """
@@ -743,23 +832,26 @@ def str_el6(coord, ul):
     return epsG.T, xl
 
 
-def str_el4(coord, ul):
+def str_el4(
+    coord: NDArray[np.float64],
+    ul: NDArray[np.float64]
+) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """Compute the strains at each element integration point
 
     This one is used for 4-noded quadrilateral elements.
 
     Parameters
     ----------
-    coord : ndarray
+    coord : NDArray[np.float64]
       Coordinates of the nodes of the element (4, 2).
-    ul : ndarray
+    ul : NDArray[np.float64]
       Array with displacements for the element.
 
     Returns
     -------
-    epsGT : ndarray
+    epsGT : NDArray[np.float64]
       Strain components for the Gauss points.
-    xl : ndarray
+    xl : NDArray[np.float64]
       Configuration of the Gauss points after deformation.
 
     """
diff --git a/solidspy/gaussutil.py b/src/solidspy/gaussutil.py
similarity index 95%
rename from solidspy/gaussutil.py
rename to src/solidspy/gaussutil.py
index b959ea8..69219e4 100755
--- a/solidspy/gaussutil.py
+++ b/src/solidspy/gaussutil.py
@@ -3,22 +3,25 @@
 Numeric integration routines
 ----------------------------
 Weights and coordinates for Gauss-Legendre quadrature [WGAUSS]_. The
-values for triangles is presented in section 5.5 of Bathe book [BATHE]_.
+values for triangles are presented in section 5.5 of Bathe's book [BATHE]_.
 
 References
 ----------
 .. [WGAUSS] Wikipedia contributors. "Gaussian quadrature." Wikipedia,
    The Free Encyclopedia, 2 Nov.  2015. Web. 25 Dec. 2015.
-   url: https://en.wikipedia.org/wiki/Gaussian_quadrature
-.. [BATHE] Bathe, Klaus-Jürgen. Finite element procedures. Prentice Hall,
-   Pearson Education, 2006.
+   URL: https://en.wikipedia.org/wiki/Gaussian_quadrature
+.. [BATHE] Bathe, Klaus-Jürgen. Finite Element Procedures.
+   Prentice Hall, Pearson Education, 2006.
 """
+
 from itertools import product
 import numpy as np
+from typing import Tuple, List
+from numpy.typing import NDArray
 
 
 #%% General
-def gauss_1d(npts):
+def gauss_1d(npts: int) -> Tuple[List[float], List[float]]:
     """Return Gauss points and weights for Gauss quadrature in 1D
 
     Parameters
@@ -28,10 +31,15 @@ def gauss_1d(npts):
 
     Returns
     -------
-    wts : ndarray
-      Weights for the Gauss-Legendre quadrature.
-    pts : ndarray
+    pts : List[float]
       Points for the Gauss-Legendre quadrature.
+    wts : List[float]
+      Weights for the Gauss-Legendre quadrature.
+
+    Raises
+    ------
+    ValueError
+        If the number of points is not between 2 and 10.
     """
     if npts == 2:
         pts = [-0.577350269189625764, 0.577350269189625764]
@@ -103,7 +111,7 @@ def gauss_1d(npts):
     return pts, wts
 
 
-def gauss_nd(npts, ndim=2):
+def gauss_nd(npts: int, ndim: int = 2) -> Tuple[NDArray[np.float64], List[float]]:
     """
     Return Gauss points and weights for Gauss quadrature in
     an ND hypercube.
@@ -111,14 +119,21 @@ def gauss_nd(npts, ndim=2):
     Parameters
     ----------
     npts : int
-      Number of quadrature points.
+      Number of quadrature points per dimension.
+    ndim : int, optional
+      Number of dimensions (default is 2).
 
     Returns
     -------
-    nd_wts : ndarray
-      Weights for the Gauss-Legendre quadrature.
-    nd_pts : ndarray
+    nd_pts : NDArray[np.float64]
       Points for the Gauss-Legendre quadrature.
+    nd_wts : List[float]
+      Weights for the Gauss-Legendre quadrature.
+
+    Raises
+    ------
+    ValueError
+        If the number of points is not between 2 and 10.
     """
     pts, wts = gauss_1d(npts)
     nd_pts = np.array(list(product(pts, repeat=ndim)))
@@ -127,20 +142,30 @@ def gauss_nd(npts, ndim=2):
     return nd_pts, nd_wts
 
 
-def gauss_tri(order=2):
+def gauss_tri(order: int = 2) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """
-    Gauss points and weights for a triangle up to order 7
+    Gauss points and weights for a triangle up to order 7.
+
+    Parameters
+    ----------
+    order : int, optional
+      Order of the quadrature rule (default is 2).
 
     Returns
     -------
-    pts : ndarray
+    pts : NDArray[np.float64]
       Points for the Gauss-Legendre quadrature.
-    wts : ndarray
+    wts : NDArray[np.float64]
       Weights for the Gauss-Legendre quadrature.
 
+    Raises
+    ------
+    ValueError
+        If the order is not between 1 and 7.
+
     References
     ----------
-    .. [BATHE] Bathe, Klaus-Jürgen. Finite element procedures.
+    .. [BATHE] Bathe, Klaus-Jürgen. Finite Element Procedures.
        Prentice Hall, Pearson Education, 2006.
     """
     if order == 1:
@@ -192,17 +217,27 @@ def gauss_tri(order=2):
     return pts, wts
 
 
-def gauss_tet(order=3):
+def gauss_tet(order: int = 3) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """
-    Gauss points and weights for a tetrahedron up to order 10
+    Gauss points and weights for a tetrahedron up to order 10.
+
+    Parameters
+    ----------
+    order : int, optional
+      Order of the quadrature rule (default is 3).
 
     Returns
     -------
-    pts : ndarray
+    pts : NDArray[np.float64]
       Points for the Gauss-Legendre quadrature.
-    wts : ndarray
+    wts : NDArray[np.float64]
       Weights for the Gauss-Legendre quadrature.
 
+    Raises
+    ------
+    ValueError
+        If the order is not between 1 and 10.
+
     References
     ----------
     .. [WITHE] Witherden, F. D., & Vincent, P. E. (2015).
@@ -213,7 +248,6 @@ def gauss_tet(order=3):
        Matthieu Ancellin, & Darius Arnold. (2020, May 5).
        nschloe/quadpy v0.14.11 (Version v0.14.11). Zenodo.
        http://doi.org/10.5281/zenodo.3786435
-
     """
     if order == 1:
         pts = np.array([[0.25, 0.25, 0.25]])
diff --git a/solidspy/postprocesor.py b/src/solidspy/postprocesor.py
similarity index 85%
rename from solidspy/postprocesor.py
rename to src/solidspy/postprocesor.py
index 3a68c13..4fc9efd 100755
--- a/solidspy/postprocesor.py
+++ b/src/solidspy/postprocesor.py
@@ -11,6 +11,8 @@
 import solidspy.uelutil as uel
 import matplotlib.pyplot as plt
 from matplotlib.tri import Triangulation
+from typing import Optional, Tuple, List
+from numpy.typing import NDArray
 
 # Set plotting defaults
 gray = '#757575'
@@ -27,7 +29,13 @@
 
 
 #%% Plotting routines
-def fields_plot(elements, nodes, disp, E_nodes=None, S_nodes=None):
+def fields_plot(
+    elements: NDArray[np.int_],
+    nodes: NDArray[np.float64],
+    disp: NDArray[np.float64],
+    E_nodes: Optional[NDArray[np.float64]] = None,
+    S_nodes: Optional[NDArray[np.float64]] = None
+) -> None:
     """Plot contours for displacements, strains and stresses
 
     Parameters
@@ -40,9 +48,9 @@ def fields_plot(elements, nodes, disp, E_nodes=None, S_nodes=None):
         to each element.
     disp : ndarray (float)
         Array with the displacements.
-    E_nodes : ndarray (float)
+    E_nodes : ndarray (float), optional
         Array with strain field in the nodes.
-    S_nodes : ndarray (float)
+    S_nodes : ndarray (float), optional
         Array with stress field in the nodes.
 
     """
@@ -76,9 +84,15 @@ def fields_plot(elements, nodes, disp, E_nodes=None, S_nodes=None):
 
 
 
-
-def tri_plot(tri, field, title="", levels=12, savefigs=False,
-             plt_type="contourf", filename="solution_plot.pdf"):
+def tri_plot(
+    tri: Triangulation,
+    field: NDArray[np.float64],
+    title: str = "",
+    levels: int = 12,
+    savefigs: bool = False,
+    plt_type: str = "contourf",
+    filename: str = "solution_plot.pdf"
+) -> None:
     """Plot contours over triangulation
 
     Parameters
@@ -88,16 +102,16 @@ def tri_plot(tri, field, title="", levels=12, savefigs=False,
         `number coordX coordY BCX BCY`
     field : ndarray (float)
         Array with data to be plotted for each node.
-    title : string (optional)
+    title : string, optional
         Title of the plot.
-    levels : int (optional)
+    levels : int, optional
         Number of levels to be used in ``contourf``.
-    savefigs : bool (optional)
+    savefigs : bool, optional
         Allow to save the figure.
-    plt_type : string (optional)
+    plt_type : string, optional
         Plot the field as one of the options: ``pcolor`` or
         ``contourf``
-    filename : string (optional)
+    filename : string, optional
         Filename to save the figures.
     """
     if plt_type == "pcolor":
@@ -112,9 +126,17 @@ def tri_plot(tri, field, title="", levels=12, savefigs=False,
         plt.savefig(filename)
 
 
-def plot_node_field(field, nodes, elements, plt_type="contourf", levels=12,
-                    savefigs=False, title=None, figtitle=None,
-                    filename=None):
+def plot_node_field(
+    field: NDArray[np.float64],
+    nodes: NDArray[np.float64],
+    elements: NDArray[np.int_],
+    plt_type: str = "contourf",
+    levels: int = 12,
+    savefigs: bool = False,
+    title: Optional[List[str]] = None,
+    figtitle: Optional[List[str]] = None,
+    filename: Optional[List[str]] = None
+) -> None:
     """Plot the nodal displacement using a triangulation
 
     Parameters
@@ -126,22 +148,21 @@ def plot_node_field(field, nodes, elements, plt_type="contourf", levels=12,
         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)
+        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)
+    levels : int, optional
         Number of levels to be used in ``contourf``.
-    savefigs : bool (optional)
+    savefigs : bool, optional
         Allow to save the figure.
-    title : Tuple of strings (optional)
+    title : Tuple of strings, optional
         Titles of the plots. If not provided the plots will not have
         a title.
-    figtitle : Tuple of strings (optional)
+    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)
+    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.
@@ -172,9 +193,18 @@ def plot_node_field(field, nodes, elements, plt_type="contourf", levels=12,
             plt.savefig(filename[cont])
 
 
-def plot_truss(nodes, elements, mats, stresses=None, max_val=4,
-               min_val=0.5, savefigs=False, title=None, figtitle=None,
-               filename=None):
+def plot_truss(
+    nodes: NDArray[np.float64],
+    elements: NDArray[np.int_],
+    mats: NDArray[np.float64],
+    stresses: Optional[NDArray[np.float64]] = None,
+    max_val: float = 4,
+    min_val: float = 0.5,
+    savefigs: bool = False,
+    title: Optional[str] = None,
+    figtitle: Optional[str] = None,
+    filename: Optional[str] = None
+) -> None:
     """Plot a truss and encodes the stresses in a colormap
 
     Parameters
@@ -191,22 +221,21 @@ def plot_truss(nodes, elements, mats, stresses=None, max_val=4,
         Array with material profiles.
     loads : ndarray (float)
         Array with loads.
-    tol : float (optional)
+    tol : float, optional
         Minimum difference between cross-section areas of the members
         to be considered different.
-    savefigs : bool (optional)
+    savefigs : bool, optional
         Allow to save the figure.
-    title : Tuple of strings (optional)
+    title : Tuple of strings, optional
         Titles of the plots. If not provided the plots will not have
         a title.
-    figtitle : Tuple of strings (optional)
+    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)
+    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.
-
     """
     min_area = mats[:, 1].min()
     max_area = mats[:, 1].max()
@@ -242,7 +271,10 @@ def plot_truss(nodes, elements, mats, stresses=None, max_val=4,
 
 
 #%% Auxiliar functions for plotting
-def mesh2tri(nodes, elements):
+def mesh2tri(
+    nodes: NDArray[np.float64],
+    elements: NDArray[np.int_]
+) -> Triangulation:
     """Generate a  matplotlib.tri.Triangulation object from the mesh
 
     Parameters
@@ -290,7 +322,12 @@ def mesh2tri(nodes, elements):
 
 
 #%% Auxiliar variables computation
-def complete_disp(bc_array, nodes, sol, ndof_node=2):
+def complete_disp(
+    bc_array: NDArray[np.int_],
+    nodes: NDArray[np.float64],
+    sol: NDArray[np.float64],
+    ndof_node: int = 2
+) -> NDArray[np.float64]:
     """
     Fill the displacement vectors with imposed and computed values.
 
@@ -311,7 +348,7 @@ def complete_disp(bc_array, nodes, sol, ndof_node=2):
 
     """
     nnodes = nodes.shape[0]
-    sol_complete = np.zeros([nnodes, ndof_node], dtype=np.float)
+    sol_complete = np.zeros([nnodes, ndof_node], dtype=float)
     for row in range(nnodes):
         for col in range(ndof_node):
             cons = bc_array[row, col]
@@ -322,7 +359,12 @@ def complete_disp(bc_array, nodes, sol, ndof_node=2):
     return sol_complete
 
 
-def strain_nodes(nodes, elements, mats, sol_complete):
+def strain_nodes(
+    nodes: NDArray[np.float64],
+    elements: NDArray[np.int_],
+    mats: NDArray[np.float64],
+    sol_complete: NDArray[np.float64]
+) -> Tuple[NDArray[np.float64], NDArray[np.float64]]:
     """Compute averaged strains and stresses at nodes
 
     First, the variable is extrapolated from the Gauss
@@ -347,14 +389,15 @@ def strain_nodes(nodes, elements, mats, sol_complete):
     -------
     E_nodes : ndarray
         Strains evaluated at the nodes.
-
+    S_nodes : ndarray
+        Stresses evaluated at the nodes.
 
     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).
+       recovery and a posteriori error estimators. Part 1. The
+       recovery technique, Int. J. Numer. Methods Eng., 33,
+       1331-1364 (1992).
     """
     nelems = elements.shape[0]
     nnodes = nodes.shape[0]
@@ -369,7 +412,7 @@ def strain_nodes(nodes, elements, mats, sol_complete):
     IELCON = elements[:, 3:]
 
     for el in range(nelems):
-        young, poisson = mats[np.int(elements[el, 2]), :]
+        young, poisson = mats[int(elements[el, 2]), :]
         shear = young/(2*(1 + poisson))
         fact1 = young/(1 - poisson**2)
         fact2 = poisson*young/(1 - poisson**2)
@@ -401,7 +444,12 @@ def strain_nodes(nodes, elements, mats, sol_complete):
     return E_nodes, S_nodes
 
 
-def stress_truss(nodes, elements, mats, disp):
+def stress_truss(
+    nodes: NDArray[np.float64],
+    elements: NDArray[np.int_],
+    mats: NDArray[np.float64],
+    disp: NDArray[np.float64]
+) -> NDArray[np.float64]:
     r"""Compute axial stresses in truss members
 
     Parameters
@@ -536,8 +584,8 @@ def stress_truss(nodes, elements, mats, disp):
     References
     ----------
     .. [1] William Weaver and James Gere. Matrix Analysis
-        of Framed Structures. Third Edition, Van Nostrand
-        Reinhold, New York (1990).
+       of Framed Structures. Third Edition, Van Nostrand
+       Reinhold, New York (1990).
 
     """
     neles = elements.shape[0]
@@ -556,14 +604,17 @@ def stress_truss(nodes, elements, mats, disp):
     return stress
 
 
-def eigvals(mat, tol=1e-6):
+def eigvals(
+    mat: NDArray[np.float64],
+    tol: float = 1e-6
+) -> Tuple[float, float, NDArray[np.float64], NDArray[np.float64]]:
     """Eigenvalues and eigenvectors for a 2x2 symmetric matrix/tensor
 
     Parameters
     ----------
     mat : ndarray
         Symmetric matrix.
-    tol : float (optional)
+    tol : float, optional
         Tolerance for considering a matrix diagonal.
 
     Returns
@@ -575,7 +626,7 @@ def eigvals(mat, tol=1e-6):
     vec1 : ndarray
         First eigenvector.
     vec2 : ndarray
-        Second eigenvector
+        Second eigenvector.
 
     Examples
     --------
@@ -617,7 +668,9 @@ def eigvals(mat, tol=1e-6):
     return eig1, eig2, vec1, vec2
 
 
-def principal_dirs(field):
+def principal_dirs(
+    field: NDArray[np.float64]
+) -> Tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]:
     """Compute the principal directions of a tensor field
 
     Parameters
@@ -656,7 +709,10 @@ def principal_dirs(field):
     return eigs1, eigs2, vecs1, vecs2
 
 
-def energy(disp, stiff):
+def energy(
+    disp: NDArray[np.float64],
+    stiff: NDArray[np.float64]
+) -> float:
     r"""
     Computes the potential energy for the current solution.
 
@@ -674,7 +730,7 @@ def energy(disp, stiff):
 
     """
     force = stiff.dot(disp)
-    el_energy = -0.5*force.dot(disp)
+    el_energy = -0.5 * force.dot(disp)
     return el_energy
 
 
diff --git a/solidspy/preprocesor.py b/src/solidspy/preprocesor.py
similarity index 81%
rename from solidspy/preprocesor.py
rename to src/solidspy/preprocesor.py
index 1c5b1c6..33918ea 100755
--- a/solidspy/preprocesor.py
+++ b/src/solidspy/preprocesor.py
@@ -9,19 +9,25 @@
 """
 import sys
 import numpy as np
+from typing import Tuple, Dict, Any, Optional
 
-
-def readin(folder=""):
+def readin(folder: str = "") -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
     """Read the input files"""
     nodes = np.loadtxt(folder + 'nodes.txt', ndmin=2)
     mats = np.loadtxt(folder + 'mater.txt', ndmin=2)
-    elements = np.loadtxt(folder + 'eles.txt', ndmin=2, dtype=np.int)
+    elements = np.loadtxt(folder + 'eles.txt', ndmin=2, dtype=int)
     loads = np.loadtxt(folder + 'loads.txt', ndmin=2)
 
     return nodes, mats, elements, loads
 
 
-def echomod(nodes, mats, elements, loads, folder=""):
+def echomod(
+    nodes: np.ndarray,
+    mats: np.ndarray,
+    elements: np.ndarray,
+    loads: np.ndarray,
+    folder: str = ""
+) -> None:
     """Create echoes of the model input files"""
     np.savetxt(folder + "KNODES.txt", nodes, fmt='%5.2f', delimiter=' ')
     np.savetxt(folder + "KMATES.txt", mats, fmt='%5.2f', delimiter=' ')
@@ -29,7 +35,7 @@ def echomod(nodes, mats, elements, loads, folder=""):
     np.savetxt(folder + "KLOADS.txt", loads, fmt='%5.2f', delimiter=' ')
 
 
-def initial_params():
+def initial_params() -> str:
     """Read initial parameters for the simulation
 
     The parameters to be read are:
@@ -58,7 +64,15 @@ def initial_params():
     return folder
 
 
-def ele_writer(cells, cell_data, ele_tag, phy_sur,  ele_type, mat_tag, nini):
+def ele_writer(
+    cells: Dict[str, np.ndarray],
+    cell_data: Dict[str, Dict[str, Any]],
+    ele_tag: str,
+    phy_sur: int,
+    ele_type: int,
+    mat_tag: int,
+    nini: int
+) -> Tuple[int, np.ndarray]:
     """
     Extracts a subset of elements from a complete mesh according to the
     physical surface  phy_sur and writes down the proper fields into an
@@ -102,16 +116,16 @@ def ele_writer(cells, cell_data, ele_tag, phy_sur,  ele_type, mat_tag, nini):
     phy_surface = cell_data[ele_tag]['gmsh:physical']
     ele_id = [cont for cont, _ in enumerate(phy_surface[:])
               if phy_surface[cont] == phy_sur]
-    els_array = np.zeros([len(ele_id) , 3 + nnode], dtype=int)
-    els_array[: , 0] = range(nini , len(ele_id) + nini )
-    els_array[: , 1] = ele_type
-    els_array[: , 2] = mat_tag
-    els_array[: , 3::] = eles[ele_id, :]
+    els_array = np.zeros([len(ele_id), 3 + nnode], dtype=int)
+    els_array[:, 0] = range(nini, len(ele_id) + nini)
+    els_array[:, 1] = ele_type
+    els_array[:, 2] = mat_tag
+    els_array[:, 3::] = eles[ele_id, :]
     nf = nini + len(ele_id)
-    return nf , els_array
+    return nf, els_array
 
 
-def node_writer(points , point_data):
+def node_writer(points: np.ndarray, point_data: Dict[str, Any]) -> np.ndarray:
     """Write nodal data as required by SolidsPy
 
     Parameters
@@ -133,7 +147,14 @@ def node_writer(points , point_data):
     return nodes_array
 
 
-def boundary_conditions(cells, cell_data, phy_lin, nodes_array, bc_x, bc_y):
+def boundary_conditions(
+    cells: Dict[str, np.ndarray],
+    cell_data: Dict[str, Dict[str, Any]],
+    phy_lin: int,
+    nodes_array: np.ndarray,
+    bc_x: int,
+    bc_y: int
+) -> np.ndarray:
     """Impose nodal point boundary conditions as required by SolidsPy
 
     Parameters
@@ -171,7 +192,13 @@ def boundary_conditions(cells, cell_data, phy_lin, nodes_array, bc_x, bc_y):
     return nodes_array
 
 
-def loading(cells, cell_data, phy_lin, P_x, P_y):
+def loading(
+    cells: Dict[str, np.ndarray],
+    cell_data: Dict[str, Dict[str, Any]],
+    phy_lin: int,
+    P_x: float,
+    P_y: float
+) -> np.ndarray:
     """Impose nodal boundary conditions as required by SolidsPy
 
     Parameters
@@ -205,12 +232,18 @@ def loading(cells, cell_data, phy_lin, P_x, P_y):
     ncargas = len(nodes_carga)
     cargas = np.zeros((ncargas, 3))
     cargas[:, 0] = nodes_carga
-    cargas[:, 1] = P_x/ncargas
-    cargas[:, 2] = P_y/ncargas
+    cargas[:, 1] = P_x / ncargas
+    cargas[:, 2] = P_y / ncargas
     return cargas
 
 
-def rect_grid(length, height, nx, ny, eletype=None):
+def rect_grid(
+    length: float,
+    height: float,
+    nx: int,
+    ny: int,
+    eletype: Optional[Any] = None
+) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
     """Generate a structured mesh for a rectangle
 
     The coordinates of the nodes will be defined in the
@@ -253,13 +286,13 @@ def rect_grid(length, height, nx, ny, eletype=None):
            [3, 1, 0, 4, 5, 8, 7]])
 
     """
-    y, x = np.mgrid[-height/2:height/2:(ny + 1)*1j,
-                    -length/2:length/2:(nx + 1)*1j]
-    els = np.zeros((nx*ny, 7), dtype=int)
+    y, x = np.mgrid[-height / 2:height / 2:(ny + 1) * 1j,
+                    -length / 2:length / 2:(nx + 1) * 1j]
+    els = np.zeros((nx * ny, 7), dtype=int)
     els[:, 1] = 1
     for row in range(ny):
         for col in range(nx):
-            cont = row*nx + col
+            cont = row * nx + col
             els[cont, 0] = cont
             els[cont, 3:7] = [cont + row, cont + row + 1,
                               cont + row + nx + 2, cont + row + nx + 1]
diff --git a/solidspy/solids_GUI.py b/src/solidspy/solids_GUI.py
similarity index 73%
rename from solidspy/solids_GUI.py
rename to src/solidspy/solids_GUI.py
index ea72120..f0337d3 100755
--- a/solidspy/solids_GUI.py
+++ b/src/solidspy/solids_GUI.py
@@ -19,33 +19,43 @@
 import solidspy.postprocesor as pos
 import solidspy.assemutil as ass
 import solidspy.solutil as sol
+from typing import Optional, Tuple, Dict, Any, Union
 
 
-def solids_GUI(plot_contours=True, compute_strains=False, folder=None):
+def solids_GUI(
+    plot_contours: bool = True,
+    compute_strains: bool = False,
+    folder: Optional[str] = None
+) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray, np.ndarray]]:
     """
     Run a complete workflow for a Finite Element Analysis
 
     Parameters
     ----------
-    plot_contours : Bool (optional)
+    plot_contours : bool, optional
         Boolean variable to plot contours of the computed variables.
         By default it is True.
-    compute_strains : Bool (optional)
+    compute_strains : bool, optional
         Boolean variable to compute Strains and Stresses at nodes.
         By default it is False.
-    folder : string (optional)
+    folder : Optional[str], optional
         String with the path to the input files. If not provided
-        it would ask for it in a pop-up window.
+        it will 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.
-
+    Union[np.ndarray, Tuple[np.ndarray, np.ndarray, np.ndarray]]
+        - If `compute_strains` is False:
+            UC : ndarray (nnodes, 2)
+                Displacements at nodes.
+        - If `compute_strains` is True:
+            Tuple containing:
+                UC : ndarray (nnodes, 2)
+                    Displacements at nodes.
+                E_nodes : ndarray (nnodes, 3)
+                    Strains at nodes.
+                S_nodes : ndarray (nnodes, 3)
+                    Stresses at nodes.
     """
     if folder is None:
         folder = pre.initial_params()
@@ -92,9 +102,12 @@ def solids_GUI(plot_contours=True, compute_strains=False, folder=None):
         return disp_complete
 
 
-
-def solids_auto(data, plot_contours=True, compute_strains=False,
-                verbose=True):
+def solids_auto(
+    data: Dict[str, Any],
+    plot_contours: bool = True,
+    compute_strains: bool = False,
+    verbose: bool = True
+) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray, np.ndarray]]:
     """
     Run a complete workflow for a Finite Element Analysis
 
@@ -103,22 +116,30 @@ def solids_auto(data, plot_contours=True, compute_strains=False,
     data : dict
         Simulation data composed of nodes, constrains, elements,
         materials and loads.
-    plot_contours : Bool (optional)
+    plot_contours : bool, optional
         Boolean variable to plot contours of the computed variables.
         By default it is True.
-    compute_strains : Bool (optional)
+    compute_strains : bool, optional
         Boolean variable to compute Strains and Stresses at nodes.
         By default it is False.
+    verbose : bool, optional
+        If True, prints detailed information during execution.
+        By default it is True.
 
     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.
-
+    Union[np.ndarray, Tuple[np.ndarray, np.ndarray, np.ndarray]]
+        - If `compute_strains` is False:
+            UC : ndarray (nnodes, 2)
+                Displacements at nodes.
+        - If `compute_strains` is True:
+            Tuple containing:
+                UC : ndarray (nnodes, 2)
+                    Displacements at nodes.
+                E_nodes : ndarray (nnodes, 3)
+                    Strains at nodes.
+                S_nodes : ndarray (nnodes, 3)
+                    Stresses at nodes.
     """
     # Retrieving data
     nodes = data["nodes"]
@@ -127,7 +148,6 @@ def solids_auto(data, plot_contours=True, compute_strains=False,
     mats = data["mats"]
     loads = data["loads"]
 
-
     # Pre-processing
     assem_op, bc_array, neq = ass.DME(cons, elements)
     if verbose:
diff --git a/solidspy/solutil.py b/src/solidspy/solutil.py
similarity index 71%
rename from solidspy/solutil.py
rename to src/solidspy/solutil.py
index 2861d08..224490a 100755
--- a/solidspy/solutil.py
+++ b/src/solidspy/solutil.py
@@ -7,31 +7,31 @@
 """
 from numpy import ndarray
 from numpy.linalg import solve
-from scipy.sparse.csr import csr_matrix
+from scipy.sparse import csr_matrix
 from scipy.sparse.linalg import spsolve
+from typing import Union
 
 
-def static_sol(mat, rhs):
+def static_sol(mat: Union[csr_matrix, ndarray], rhs: ndarray) -> ndarray:
     """Solve a static problem [mat]{u_sol} = {rhs}
     
     Parameters
     ----------
-    mat : array
+    mat : Union[csr_matrix, ndarray]
         Array with the system of equations. It can be stored in
         dense or sparse scheme.
-    rhs : array
+    rhs : ndarray
         Array with right-hand-side of the system of equations.
     
     Returns
     -------
-    u_sol : array
+    u_sol : ndarray
         Solution of the system of equations.
 
     Raises
     ------
-    
-        
-
+    TypeError
+        If `mat` is neither a `numpy.ndarray` nor a `scipy.sparse.csr_matrix`.
     """
     if type(mat) is csr_matrix:
         u_sol = spsolve(mat, rhs)
diff --git a/solidspy/uelutil.py b/src/solidspy/uelutil.py
similarity index 86%
rename from solidspy/uelutil.py
rename to src/solidspy/uelutil.py
index b438582..2620982 100755
--- a/solidspy/uelutil.py
+++ b/src/solidspy/uelutil.py
@@ -11,19 +11,24 @@
 import numpy as np
 import solidspy.femutil as fem
 import solidspy.gaussutil as gau
+from numpy import ndarray
+from typing import Tuple, Sequence
 
 
 #%% Continuum elements
 
 # Triangles
-def elast_tri3(coord, params):
+def elast_tri3(
+    coord: ndarray, 
+    params: Sequence[float]
+) -> Tuple[ndarray, ndarray]:
     """Triangular element with 3 nodes
 
     Parameters
     ----------
     coord : ndarray
       Coordinates for the nodes of the element (3, 2).
-    params : tuple
+    params : Sequence[float]
       Material parameters in the following order:
 
           young : float
@@ -48,7 +53,7 @@ def elast_tri3(coord, params):
     ...         [1, 0],
     ...         [0, 1]])
     >>> params = [8/3, 1/3]
-    >>> stiff, mass = uel3ntrian(coord, params)
+    >>> stiff, mass = elast_tri3(coord, params)
     >>> stiff_ex = 1/2 * np.array([
     ...            [4, 2, -3, -1, -1, -1],
     ...            [2, 4, -1, -1, -1, -3],
@@ -77,14 +82,17 @@ def elast_tri3(coord, params):
     return stiff_mat, mass_mat
 
 
-def elast_tri6(coord, params):
+def elast_tri6(
+    coord: ndarray, 
+    params: Sequence[float]
+) -> Tuple[ndarray, ndarray]:
     """Triangular element with 6 nodes
 
     Parameters
     ----------
     coord : ndarray
         Coordinates for the nodes of the element (6, 2).
-    params : tuple
+    params : Sequence[float]
         Material parameters in the following order:
 
             young : float
@@ -112,7 +120,7 @@ def elast_tri6(coord, params):
     ...         [0.5, 0.5],
     ...         [0, 0.5]])
     >>> params = [8/3, 1/3]
-    >>> stiff, mass = uel6ntrian(coord, params)
+    >>> stiff, mass = elast_tri6(coord, params)
     >>> 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],
@@ -148,14 +156,17 @@ def elast_tri6(coord, params):
 
 
 # Quadrilaterals
-def elast_quad4(coord, params):
+def elast_quad4(
+    coord: ndarray, 
+    params: Sequence[float]
+) -> Tuple[ndarray, ndarray]:
     """Quadrilateral element with 4 nodes
 
     Parameters
     ----------
     coord : ndarray
         Coordinates for the nodes of the element (4, 2).
-    params : tuple
+    params : Sequence[float]
         Material parameters in the following order:
 
             young : float
@@ -177,7 +188,7 @@ def elast_quad4(coord, params):
 
     >>> coord = np.array([[-1, -1], [1, -1], [1, 1], [-1, 1]])
     >>> params = [8/3, 1/3, 1]
-    >>> stiff, mass = uel4nquad(coord, params)
+    >>> stiff, mass = elast_quad4(coord, params)
     >>> stiff_ex = 1/6 * np.array([
     ...             [ 8,  3, -5,  0, -4, -3,  1,  0],
     ...             [ 3,  8,  0,  1, -3, -4,  0, -5],
@@ -215,27 +226,30 @@ def elast_quad4(coord, params):
         H, B, det = fem.elast_diff_2d(r, s, coord, fem.shape_quad4)
         factor = det * gwts[cont]
         stiff_mat += factor * (B.T @ C @ B)
-        mass_mat += dens*factor* (H.T @ H)
+        mass_mat += dens * factor * (H.T @ H)
     return stiff_mat, mass_mat
 
 
-def elast_quad9(coord, params):
+def elast_quad9(
+    coord: ndarray, 
+    params: Sequence[float]
+) -> Tuple[ndarray, ndarray]:
     """
     Quadrilateral element with 9 nodes for classic elasticity
     under plane-strain
 
     Parameters
     ----------
-    coord : coord
+    coord : ndarray
         Coordinates of the element.
-    params : list
+    params : Sequence[float]
         List with material parameters in the following order:
         [Young modulus, Poisson coefficient, density].
 
     Returns
     -------
     stiff_mat : ndarray (float)
-        Local stifness matrix.
+        Local stiffness matrix.
     mass_mat : ndarray (float)
         Local mass matrix.
     """
@@ -256,26 +270,30 @@ def elast_quad9(coord, params):
     return stiff_mat, mass_mat
 
 
-def elast_quad8(coord, params):
+def elast_quad8(
+    coord: ndarray, 
+    params: Sequence[float]
+) -> Tuple[ndarray, ndarray]:
     """
-    Quadrilateral element with 9 nodes for classic elasticity
+    Quadrilateral element with 8 nodes for classic elasticity
     under plane-strain
 
     Parameters
     ----------
-    coord : coord
+    coord : ndarray
         Coordinates of the element.
-    params : list
+    params : Sequence[float]
         List with material parameters in the following order:
         [Young modulus, Poisson coefficient, density].
 
     Returns
     -------
     stiff_mat : ndarray (float)
-        Local stifness matrix.
+        Local stiffness matrix.
     mass_mat : ndarray (float)
         Local mass matrix.
     """
+    
     E, nu, rho = params
     C = fem.umat((E, nu))
     stiff_mat = np.zeros((16, 16))
@@ -287,28 +305,31 @@ def elast_quad8(coord, params):
         H, B, det = fem.elast_diff_2d(r, s, coord, fem.shape_quad8)
         factor = det * gwts[cont]
         stiff_mat += factor * (B.T @ C @ B)
-        mass_mat += rho*factor * (H.T @ H)
+        mass_mat += rho * factor * (H.T @ H)
     return stiff_mat, mass_mat
 
 
 ## Axisymmetric
-def elast_axi_quad9(coord, params):
+def elast_axi_quad9(
+    coord: ndarray, 
+    params: Sequence[float]
+) -> Tuple[ndarray, ndarray]:
     """
     Quadrilateral element with 9 nodes for classic elasticity
     under axisymmetric conditions
 
     Parameters
     ----------
-    coord : coord
+    coord : ndarray
         Coordinates of the element.
-    params : list
+    params : Sequence[float]
         List with material parameters in the following order:
         [Young modulus, Poisson coefficient, density].
 
     Returns
     -------
     stiff_mat : ndarray (float)
-        Local stifness matrix.
+        Local stiffness matrix.
     mass_mat : ndarray (float)
         Local mass matrix.
     """
@@ -324,25 +345,29 @@ def elast_axi_quad9(coord, params):
         H, B, det = fem.elast_diff_axi(r, s, coord, fem.shape_quad9)
         factor = det * gwts[cont]
         stiff_mat += factor * (B.T @ C @ B)
-        mass_mat += rho*factor * (H.T @ H)
+        mass_mat += rho * factor * (H.T @ H)
     return stiff_mat, mass_mat
 
+
 ## 3D elements
-def elast_tet4(coord, params):
+def elast_tet4(
+    coord: ndarray, 
+    params: Sequence[float]
+) -> Tuple[ndarray, ndarray]:
     """Tetraedral element with 4 nodes for classic elasticity
 
     Parameters
     ----------
-    coord : coord
+    coord : ndarray
         Coordinates of the element.
-    params : list
+    params : Sequence[float]
         List with material parameters in the following order:
         [Young modulus, Poisson coefficient, density].
 
     Returns
     -------
     stiff_mat : ndarray (float)
-        Local stifness matrix.
+        Local stiffness matrix.
     mass_mat : ndarray (float)
         Local mass matrix.
     """
@@ -358,25 +383,28 @@ def elast_tet4(coord, params):
         H, B, det = fem.elast_diff_3d(r, s, t, coord, fem.shape_tet4)
         factor = det * gwts[cont] / 6
         stiff_mat += factor * (B.T @ C @ B)
-        mass_mat += rho*factor * (H.T @ H)
+        mass_mat += rho * factor * (H.T @ H)
     return stiff_mat, mass_mat
 
 
-def elast_hex8(coord, params):
+def elast_hex8(
+    coord: ndarray, 
+    params: Sequence[float]
+) -> Tuple[ndarray, ndarray]:
     """Hexaedral element with 8 nodes for classic elasticity
 
     Parameters
     ----------
-    coord : coord
+    coord : ndarray
         Coordinates of the element.
-    params : list
+    params : Sequence[float]
         List with material parameters in the following order:
         [Young modulus, Poisson coefficient, density].
 
     Returns
     -------
     stiff_mat : ndarray (float)
-        Local stifness matrix.
+        Local stiffness matrix.
     mass_mat : ndarray (float)
         Local mass matrix.
     """
@@ -392,12 +420,15 @@ def elast_hex8(coord, params):
         H, B, det = fem.elast_diff_3d(r, s, t, coord)
         factor = det * gwts[cont]
         stiff_mat += factor * (B.T @ C @ B)
-        mass_mat += rho*factor * (H.T @ H)
+        mass_mat += rho * factor * (H.T @ H)
     return stiff_mat, mass_mat
 
 
 #%% Structural elements
-def spring(coord, stiff):
+def spring(
+    coord: ndarray, 
+    stiff: float
+) -> Tuple[ndarray, ndarray]:
     """1D 2-noded Spring element
 
     Parameters
@@ -422,7 +453,8 @@ def spring(coord, stiff):
     >>> coord = np.array([
     ...         [0, 0],
     ...         [1, 0]])
-    >>> stiff, mass = uelspring(coord, 8/3)
+    >>> stiff = 8/3
+    >>> stiff, mass = spring(coord, stiff)
     >>> stiff_ex = 8/3 * np.array([
     ...    [1, 0, -1, 0],
     ...    [0, 0, 0, 0],
@@ -445,14 +477,17 @@ def spring(coord, stiff):
     return stiff_mat, np.zeros(4)
 
 
-def truss2D(coord, params):
+def truss2D(
+    coord: ndarray, 
+    params: Sequence[float]
+) -> Tuple[ndarray, ndarray]:
     """2D 2-noded truss element
 
     Parameters
     ----------
     coord : ndarray
       Coordinates for the nodes of the element (2, 2).
-    params : tuple
+    params : Sequence[float]
       Element parameters in the following order:
 
           young : float
@@ -476,7 +511,7 @@ def truss2D(coord, params):
     ...         [0, 0],
     ...         [1, 0]])
     >>> params = [1.0 , 1.0]
-    >>> stiff, mass = ueltruss2D(coord, params)
+    >>> stiff, mass = truss2D(coord, params)
     >>> stiff_ex =  np.array([
     ...    [1, 0, -1, 0],
     ...    [0, 0, 0, 0],
@@ -494,7 +529,7 @@ def truss2D(coord, params):
         [nx, ny, 0, 0],
         [0, 0, nx, ny]])
     young, area  = params[:2]
-    stiff = area*young/length
+    stiff = area * young / length
     stiff_mat = stiff * np.array([
         [1, -1],
         [-1, 1]])
@@ -504,7 +539,7 @@ def truss2D(coord, params):
     else:
         dens = params[-1]
     mass = area * length * dens
-    mass_mat = mass/6*np.array([
+    mass_mat = mass / 6 * np.array([
             [2, 0, 1, 0],
             [0, 2, 0, 1],
             [1, 0, 2, 0],
@@ -512,14 +547,17 @@ def truss2D(coord, params):
     return stiff_mat, mass_mat
 
 
-def beam2DU(coord, params):
+def beam2DU(
+    coord: ndarray, 
+    params: Sequence[float]
+) -> Tuple[ndarray, ndarray]:
     """2D 2-noded beam element without axial deformation
 
     Parameters
     ----------
     coord : ndarray
         Coordinates for the nodes of the element (2, 2).
-    params : tuple
+    params : Sequence[float]
         Element parameters in the following order:
     
             young : float
@@ -548,7 +586,7 @@ def beam2DU(coord, params):
         [0, 0, 0, ny, nx, 0],
         [0, 0, 0, 0, 0, 1]])
     young, area_moment = params[:2]
-    bending_stiff = area_moment*young/L**3
+    bending_stiff = area_moment * young / L**3
     stiff_mat = bending_stiff * np.array([
         [12, 6*L, -12, 6*L],
         [6*L, 4*L**2, -6*L, 2*L**2],
@@ -560,7 +598,7 @@ def beam2DU(coord, params):
     else:
         dens, area = params[2:]
     mass = area * L * dens
-    mass_mat = mass/420*np.array([
+    mass_mat = mass / 420 * np.array([
             [156, 22*L, 54, -13*L],
             [22*L, 4*L**2, 13*L, -3*L**2],
             [54, 13*L, 156, -22*L],
@@ -570,14 +608,17 @@ def beam2DU(coord, params):
     return stiff_mat, mass_mat
 
 
-def beam2D(coord, params):
+def beam2D(
+    coord: ndarray, 
+    params: Sequence[float]
+) -> Tuple[ndarray, ndarray]:
     """2D 2-noded beam element with axial deformation
 
     Parameters
     ----------
     coord : ndarray
       Coordinates for the nodes of the element (2, 2).
-    params : tuple
+    params : Sequence[float]
         Element parameters in the following order:
     
             young : float
@@ -609,8 +650,8 @@ def beam2D(coord, params):
         [0, 0, 0, 0, 0, 1]])
     young, area_moment, area = params[:3]
     bending_stiff = area_moment * young
-    ratio = area/area_moment
-    stiff_mat = bending_stiff/L**3 * np.array([
+    ratio = area / area_moment
+    stiff_mat = bending_stiff / L**3 * np.array([
         [ratio*L**2, 0, 0, -ratio*L**2, 0, 0],
         [0, 12, 6*L, 0, -12, 6*L],
         [0, 6*L, 4*L**2, 0, -6*L, 2*L**2],
@@ -623,7 +664,7 @@ def beam2D(coord, params):
     else:
         dens = params[3:]
     mass = area * L * dens
-    mass_mat = mass/420*np.array([
+    mass_mat = mass / 420 * np.array([
             [140, 0, 0, 70, 0, 0],
             [0, 156, 22*L, 0, 54, -13*L],
             [0, 22*L, 4*L**2, 0, 13*L, -3*L**2],


Test patch:
diff --git a/tests/test_assemutil.py b/tests/test_assemutil.py
index b96f956..fc005a8 100755
--- a/tests/test_assemutil.py
+++ b/tests/test_assemutil.py
@@ -30,7 +30,7 @@ def test_sparse_assem():
         [0, 1, 2, 3, 8, 9, 6, 7],
         [2, 3, 4, 5, 10, 11, 8, 9],
         [6, 7, 8, 9, 14, 15, 12, 13],
-        [8, 9, 10, 11, 16, 17, 14, 15]], dtype=np.int)
+        [8, 9, 10, 11, 16, 17, 14, 15]], dtype=int)
 
     neq = 18
     K_ass, _ = ass.sparse_assem(elements, mats, nodes, neq, assem_op)
@@ -80,7 +80,7 @@ def test_dense_assem():
         [0, 1, 2, 3, 8, 9, 6, 7],
         [2, 3, 4, 5, 10, 11, 8, 9],
         [6, 7, 8, 9, 14, 15, 12, 13],
-        [8, 9, 10, 11, 16, 17, 14, 15]], dtype=np.int)
+        [8, 9, 10, 11, 16, 17, 14, 15]], dtype=int)
 
     neq = 18
     K_ass, _ = ass.dense_assem(elements, mats, nodes, neq, assem_op)
@@ -124,7 +124,7 @@ def uel_ones(elcoord, params):
         [1, 2],
         [2, 2]])
     cons = np.zeros((9, 2))
-    elements = np.ones((4, 7), dtype=np.int)
+    elements = np.ones((4, 7), dtype=int)
     elements[:, 0] = range(0, 4)
     elements[:, 2] = 0
     elements[:, 3:] = np.array([