mesh-demo

Ansh for defining unstructured meshes¶

ansh is a Python package to define and visualise finite-element scalar and vector fields. It is a thin wrapper around meshio, pyvista, and skfem. meshio is used for reading and writing meshes. skfem is used for finite-element field definition. pyvista is used from plotting both meshes and fields.

The notebook shows mesh specific functionality of Ansh.

In [1]:

import ansh
import numpy as np
import pyvista as pv

In [2]:

pv.set_jupyter_backend("html")  # This is only for website; remove when running locally

In [3]:

m = ansh.Mesh.from_file("../tests/meshes/med/omega.med")

In [4]:

m.info = {"name": "shell", "about": "one of the unit test meshes"}

Basic mesh properties¶

In [5]:

Out[5]:

Ansh Mesh

Bounds: [-2.999, -2.998, -3.0] - [2.998, 2.997, 3.0]
Points: 11669
Cells: 72833
- line: 78
- tetra: 63697
- triangle: 9058
Subregions: 4
- core:
  - line: 0
  - tetra: 22577
  - triangle: 0
- shell:
  - line: 0
  - tetra: 41120
  - triangle: 0
- boundary:
  - line: 0
  - tetra: 0
  - triangle: 6238
- 0:
  - line: 78
  - tetra: 0
  - triangle: 2820
Scale: 1.0
Info: 2 items
- name: shell
- about: one of the unit test meshes

Meaning of basic mesh attributes:¶

Bounds: mesh bounds give an idea about the size of the mesh.

Points: points of the mesh in 3D.

In [6]:

m.points

Out[6]:

array([[ 1.22464680e-16, -2.99951957e-32, -2.00000000e+00],
       [ 1.22464680e-16, -2.99951957e-32,  2.00000000e+00],
       [ 1.83697020e-16, -4.49927935e-32, -3.00000000e+00],
       ...,
       [-1.42971150e+00,  2.56445340e+00,  3.66067141e-01],
       [-1.25520072e+00,  2.24344656e+00,  7.50555066e-01],
       [-1.37227559e+00,  2.28844921e+00,  1.26925109e+00]],
      shape=(11669, 3))

Cells: cell blocks or list of cell blocks that has information about cell type and connectivity.

In [7]:

m.cells

Out[7]:

[CellBlock(type: line, cells: 78),
 CellBlock(type: tetra, cells: 63697),
 CellBlock(type: triangle, cells: 9058)]

In [8]:

m.cells[1]

Out[8]:

CellBlock(type: tetra, cells: 63697)

In [9]:

m.cells[1].type

Out[9]:

'tetra'

In [10]:

m.cells[1].connectivity

Out[10]:

array([[  515,  6864,   438,   422],
       [  241,   317,  4890,   240],
       [ 7249,   901,  4894,  7247],
       ...,
       [ 3842,  3695,  9014, 11668],
       [ 9014, 11668,  3695, 11620],
       [ 9014,  3842, 11668, 11620]], shape=(63697, 4))

Subregions: dictionary defining subregions of the mesh with region name as keys and cell sets as values. Cell sets are indices of the cells for each entry in Mesh.cells list, that belong to the subregion.

In [11]:

m.subregions

Out[11]:

{'core': [array([], dtype=int64),
  array([    0,     1,     2, ..., 22574, 22575, 22576], shape=(22577,)),
  array([], dtype=int64)],
 'shell': [array([], dtype=int64),
  array([22577, 22578, 22579, ..., 63694, 63695, 63696], shape=(41120,)),
  array([], dtype=int64)],
 'boundary': [array([], dtype=int64),
  array([], dtype=int64),
  array([2820, 2821, 2822, ..., 9055, 9056, 9057], shape=(6238,))],
 '0': [array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
         17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
         34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
         51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67,
         68, 69, 70, 71, 72, 73, 74, 75, 76, 77]),
  array([], dtype=int64),
  array([   0,    1,    2, ..., 2817, 2818, 2819], shape=(2820,))]}

In [12]:

m.subregions["core"]

Out[12]:

[array([], dtype=int64),
 array([    0,     1,     2, ..., 22574, 22575, 22576], shape=(22577,)),
 array([], dtype=int64)]

Scale: the scale or a multiplying factor for the mesh points.

Info: additional metadata attached to the mesh object.

In [13]:

m.info

Out[13]:

{'name': 'shell', 'about': 'one of the unit test meshes'}

Additional mesh properties and methods¶

Cell centres: list of cell centre coordinate arrays corresponding to each cell block entry in Mesh.cells

In [14]:

type(m.cell_centres)

Out[14]:

list

In [15]:

len(m.cell_centres)

Out[15]:

In [16]:

m.cells[1]

Out[16]:

CellBlock(type: tetra, cells: 63697)

In [17]:

m.cell_centres[1]

Out[17]:

array([[ 1.55419836, -0.89677233,  0.73276849],
       [ 1.71388184,  0.54067495,  0.76760184],
       [ 0.89402151,  1.59210337,  0.24942144],
       ...,
       [-1.29613864,  2.29507797,  1.29634942],
       [-1.29087363,  2.25669716,  1.29889778],
       [-1.26516341,  2.29243455,  1.2626589 ]], shape=(63697, 3))

to PyVista: returns a pyvista.UnstructuredGrid. If mesh subregions are defined, they stored as cell data in the form of boolean arrays that can act as masks for selecting the cells. Additionally, the mesh information is attached as Field Data.

In [18]:

m.to_pyvista()

Out[18]:

UnstructuredGrid (0x7f01b8ad6200)
  N Cells:    72833
  N Points:   11669
  X Bounds:   -2.999e+00, 2.998e+00
  Y Bounds:   -2.998e+00, 2.997e+00
  Z Bounds:   -3.000e+00, 3.000e+00
  N Arrays:   8

In [19]:

m.to_pyvista()["core"]

Out[19]:

pyvista_ndarray([False, False, False, ..., False, False, False],
                shape=(72833,))

In [20]:

m.to_pyvista().field_data

Out[20]:

pyvista DataSetAttributes
Association     : NONE
Contains arrays :
    AnshMesh                bool       (1,)
    AnshScale               float64    (1,)
    name                    str        "shell"
    about                   str        "one of the unit..."

Note: the Field Data shows an extra entry AnshMesh set to True. This key marks the UnstructuredGrid as created by Ansh.

Plot: convenience method to plot the mesh based on PyVista's plotting functionality.

In [21]:

m.plot(subregions=["core"])

Add subregion: method to add a subregion to the mesh by passing indices of the cells for corresponding cell blocks or selector functions the return a boolian mask. The selector functions ingest cell centres corresponding to each cell block.

In [22]:

def select_unit_cube(centres):
    # Here centres is an array of cell centres for a cell block
    return (
        ((centres[:, 0] < 1.0) & (centres[:, 0] > -1.0))  # x-axis
        & ((centres[:, 1] < 1.0) & (centres[:, 1] > -1.0))  # y-axis
        & ((centres[:, 2] < 1.0) & (centres[:, 2] > -1.0))  # z-axis
    )

In [23]:

m.add_subregion(
    subregion_name="cube",
    value=[np.array([], dtype=np.int_), select_unit_cube, np.array([], dtype=np.int_)],
)

In [24]:

m.plot(subregions=["cube"], cmap="viridis")

Note: since the mesh was not aligned with the cube geometry, the selected cube subregion will show crinkle. Hence, it is always advisable to create subregions during mesh generation.

In [25]:

m.subregions.keys()

Out[25]:

dict_keys(['core', 'shell', 'boundary', '0', 'cube'])

Mesh conversion to meshio and skfem: convert ansh mesh to meshio and skfem mesh.

In [26]:

m.to_meshio()

Out[26]:

<meshio mesh object>
  Number of points: 11669
  Number of cells:
    line: 78
    tetra: 63697
    triangle: 9058
  Cell sets: core, shell, boundary, 0, cube

In [27]:

ansh.Mesh.from_meshio(m.to_meshio())

Out[27]:

Ansh Mesh

Bounds: [-2.999, -2.998, -3.0] - [2.998, 2.997, 3.0]
Points: 11669
Cells: 72833
- line: 78
- tetra: 63697
- triangle: 9058
Subregions: 5
- core:
  - line: 0
  - tetra: 22577
  - triangle: 0
- shell:
  - line: 0
  - tetra: 41120
  - triangle: 0
- boundary:
  - line: 0
  - tetra: 0
  - triangle: 6238
- 0:
  - line: 78
  - tetra: 0
  - triangle: 2820
- cube:
  - line: 0
  - tetra: 6640
  - triangle: 0
Scale: 1.0
Info: 2 items
- name: shell
- about: one of the unit test meshes

In [28]:

m.to_skfem()

Out[28]:

<skfem MeshTet1 object>
  Number of elements: 63697
  Number of vertices: 11669
  Number of nodes: 11669
  Named subdomains [# elements]: core [22577], shell [41120], cube [6640]
  Named boundaries [# facets]: boundary [6238], 0 [2820]

Mesh read write¶

Read mesh: mesh file reading is delegated to the powerful meshio. Hence, (theoretically) it is possible to read any file that meshio can read. However, ansh.Mesh sanitises the meshio.Mesh object for inter-operability with pyvista and skfem. So, there is a possibility of a read failure. Please open an issue if your mesh is not read properly. Additionally, one can read HDF5 (.h5) files if they were saved using ansh.Mesh.to_file().

In [29]:

read_mesh = ansh.Mesh.from_file("../tests/meshes/netgen/periodic_3d.vol")
read_mesh

Out[29]:

Ansh Mesh

Bounds: [0.0, 0.0, 0.0] - [1.0, 1.0, 1.0]
Points: 65
Cells: 346
- triangle: 108
- tetra: 166
- line: 72
Subregions: 6
- 1:
  - triangle: 20
  - tetra: 166
  - line: 12
- 2:
  - triangle: 16
  - tetra: 0
  - line: 12
- default:
  - triangle: 16
  - tetra: 0
  - line: 12
- 4:
  - triangle: 20
  - tetra: 0
  - line: 12
- 5:
  - triangle: 16
  - tetra: 0
  - line: 12
- outer:
  - triangle: 20
  - tetra: 0
  - line: 12
Scale: 1.0
Info: 2 items
- netgen:identifications: [[ 1 3 1] [ 2 5 1] [ 4 7 1] [ 6 8 1] [ 9 11 1] [10 12 1] [15 13 1] [16 14 1] [21 19 1] [22 20 1] [25 23 1] [26 24 1] [38 54 1] [39 55 1] [40 56 1]]
- netgen:identificationtypes: [[2]]

In [30]:

read_mesh.plot(subregions=["2", "4"])

Write mesh: ansh.Mesh can be written to an Abaqus .inp, Ansh HDF5 .h5, or a VTK unstructured grid .vtu format. The Abaqus and HDF5 format aligns closely with the outline of ansh.Mesh object. The VTK format requires a conversion of cell sets to cell data for the subregion definition. The aim of Ansh is not to convert from one mesh file format to another. Hence, if you are unhappy with the supported formats, please use meshio to convert to your preferred file format.

In [31]:

read_mesh.to_file("read.inp")

In [32]:

ansh.Mesh.from_file("read.inp")

Out[32]:

Ansh Mesh

Bounds: [0.0, 0.0, 0.0] - [1.0, 1.0, 1.0]
Points: 65
Cells: 346
- triangle: 108
- tetra: 166
- line: 72
Subregions: 6
- 1:
  - triangle: 0
  - tetra: 0
  - line: 12
- 2:
  - triangle: 0
  - tetra: 0
  - line: 12
- default:
  - triangle: 0
  - tetra: 0
  - line: 12
- 4:
  - triangle: 0
  - tetra: 0
  - line: 12
- 5:
  - triangle: 0
  - tetra: 0
  - line: 12
- outer:
  - triangle: 0
  - tetra: 0
  - line: 12
Scale: 1.0

In [33]:

read_mesh.to_file("read.h5")

In [34]:

!h5glance --attrs read.h5

read.h5
└mesh
  ├3 attributes:
  │ ├netgen:identifications: array [int64: (15, 3)]
  │ ├netgen:identificationtypes: array [int64: (1, 1)]
  │ └scale: np.float64(1.0)
  ├cells
  │ ├0
  │ │ └triangle  [uint32: 108 × 3]
  │ ├1
  │ │ └tetra [uint32: 166 × 4]
  │ └2
  │   └line  [uint32: 72 × 2]
  ├points    [float64: 65 × 3]
  └subregions
    ├1
    │ ├0 [int64: 20]
    │ ├1 [int64: 166]
    │ └2 [int64: 12]
    ├2
    │ ├0 [int64: 16]
    │ ├1 [int64: 0]
    │ └2 [int64: 12]
    ├4
    │ ├0 [int64: 20]
    │ ├1 [int64: 0]
    │ └2 [int64: 12]
    ├5
    │ ├0 [int64: 16]
    │ ├1 [int64: 0]
    │ └2 [int64: 12]
    ├default
    │ ├0 [int64: 16]
    │ ├1 [int64: 0]
    │ └2 [int64: 12]
    └outer
      ├0 [int64: 20]
      ├1 [int64: 0]
      └2 [int64: 12]

In [35]:

ansh.Mesh.from_file("read.h5")

Out[35]:

Ansh Mesh

Bounds: [0.0, 0.0, 0.0] - [1.0, 1.0, 1.0]
Points: 65
Cells: 346
- triangle: 108
- tetra: 166
- line: 72
Subregions: 6
- 1:
  - triangle: 20
  - tetra: 166
  - line: 12
- 2:
  - triangle: 16
  - tetra: 0
  - line: 12
- 4:
  - triangle: 20
  - tetra: 0
  - line: 12
- 5:
  - triangle: 16
  - tetra: 0
  - line: 12
- default:
  - triangle: 16
  - tetra: 0
  - line: 12
- outer:
  - triangle: 20
  - tetra: 0
  - line: 12
Scale: 1.0
Info: 2 items
- netgen:identifications: [[ 1 3 1] [ 2 5 1] [ 4 7 1] [ 6 8 1] [ 9 11 1] [10 12 1] [15 13 1] [16 14 1] [21 19 1] [22 20 1] [25 23 1] [26 24 1] [38 54 1] [39 55 1] [40 56 1]]
- netgen:identificationtypes: [[2]]

In [36]:

read_mesh.to_file("read.vtu")

In [37]:

ansh.Mesh.from_file("read.vtu")

Out[37]:

Ansh Mesh

Bounds: [0.0, 0.0, 0.0] - [1.0, 1.0, 1.0]
Points: 65
Cells: 346
- triangle: 108
- tetra: 166
- line: 72
Subregions: 6
- 1:
  - triangle: 20
  - tetra: 166
  - line: 12
- 2:
  - triangle: 16
  - tetra: 0
  - line: 12
- default:
  - triangle: 16
  - tetra: 0
  - line: 12
- 4:
  - triangle: 20
  - tetra: 0
  - line: 12
- 5:
  - triangle: 16
  - tetra: 0
  - line: 12
- outer:
  - triangle: 20
  - tetra: 0
  - line: 12
Scale: 1.0
Info: 3 items
- AnshScale: [1.]
- netgen:identifications: [[ 1 3 1] [ 2 5 1] [ 4 7 1] [ 6 8 1] [ 9 11 1] [10 12 1] [15 13 1] [16 14 1] [21 19 1] [22 20 1] [25 23 1] [26 24 1] [38 54 1] [39 55 1] [40 56 1]]
- netgen:identificationtypes: [2]