vector-field-demo

Ansh for defining vector fields¶

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 how to define vector field on a mesh.

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"}

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
- boundary:
  - line: 0
  - tetra: 0
  - triangle: 6238
- shell:
  - line: 0
  - tetra: 41120
  - triangle: 0
- core:
  - line: 0
  - tetra: 22577
  - triangle: 0
- 0:
  - line: 78
  - tetra: 0
  - triangle: 2820
Scale: 1.0
Info: 2 items
- name: shell
- about: one of the unit test meshes

In [6]:

f = ansh.VectorField(
    mesh=m, element_type=ansh.element_types.TetP1(), value=[42.0, 42.0, 42.0]
)

In [7]:

Out[7]:

Ansh VectorField

Mesh: 11669 points, 72833 cells, 4 subregions
- 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
  - boundary:
    
    line: 0
    tetra: 0
    triangle: 6238
  - shell:
    
    line: 0
    tetra: 41120
    triangle: 0
  - core:
    
    line: 0
    tetra: 22577
    triangle: 0
  - 0:
    
    line: 78
    tetra: 0
    triangle: 2820
- Scale: 1.0
- Info: 2 items
  - name: shell
  - about: one of the unit test meshes
Element: ElementTetP1
Value: min=42.0, max=42.0
array([[42., 42., 42.], [42., 42., 42.], [42., 42., 42.], ..., [42., 42., 42.], [42., 42., 42.], [42., 42., 42.]], shape=(11669, 3))

Basic vector field properties¶

Mesh: the Ansh mesh on which the field is defined.

In [8]:

f.mesh

Out[8]:

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
- boundary:
  - line: 0
  - tetra: 0
  - triangle: 6238
- shell:
  - line: 0
  - tetra: 41120
  - triangle: 0
- core:
  - line: 0
  - tetra: 22577
  - triangle: 0
- 0:
  - line: 78
  - tetra: 0
  - triangle: 2820
Scale: 1.0
Info: 2 items
- name: shell
- about: one of the unit test meshes

In [9]:

f.mesh is m

Out[9]:

True

Element type: the skfem element type used to approximate the field.

In [10]:

f.element_type

Out[10]:

<skfem.element.element_vector.ElementVector at 0x7f84aaf83260>

Note: ansh.element_types is a convenience named tuple that exposes all the skfem element types to Ansh. To see all the available element types, hit Tab after ansh.element_types. in a jupyter notebook (or an ipython shell). The auto-complete should also work with popular Python LSPs in your preferred IDE.

In [11]:

f.element_type.elem

Out[11]:

<skfem.element.element_tet.element_tet_p1.ElementTetP1 at 0x7f84aaff3a70>

Value: the value of the vector field on degree of freedom (DOF) points. For example, for polynomial basis functions of degree one (like TetP1), DOFs will coincide with the mesh nodes.

In [12]:

f.value

Out[12]:

array([[42., 42., 42.],
       [42., 42., 42.],
       [42., 42., 42.],
       ...,
       [42., 42., 42.],
       [42., 42., 42.],
       [42., 42., 42.]], shape=(11669, 3))

In [13]:

f.mesh.points.shape[0]

Out[13]:

Additional vector field properties and methods¶

Degree of freedom locations: the points at which the value of the vector field approximation is allowed to change. The value of field is interpolated in rest of the geometry using the basis functions.

In [14]:

f.dof_locations

Out[14]:

array([[ 1.24900090e-16,  0.00000000e+00, -2.00000000e+00],
       [ 1.11022302e-16,  0.00000000e+00,  2.00000000e+00],
       [ 1.83880688e-16,  0.00000000e+00, -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))

Value at any point: the value of the vector field approximations can be obtained by calling the field over the coordinates of one or more points inside the mesh.

In [15]:

f([0, 0, 0])

Out[15]:

array([[42., 42., 42.]])

In [16]:

f(np.array([[0, 2.0, -1.0], [0, 0, 0], [-2, -1, 1.9]]))

Out[16]:

array([[42., 42., 42.],
       [42., 42., 42.],
       [42., 42., 42.]])

To skfem: create a skfem.DiscreteField from ansh.VectorField for further processing such as matrix assembly. It is also possible to just extract skfem basis.

In [17]:

f.to_skfem()

Out[17]:

<skfem DiscreteField object>
  Quadrature points per element: 4
  Number of elements: 63697
  Order: 1
  Attributes: grad

In [18]:

f.skfem_basis

Out[18]:

<skfem CellBasis(MeshTet1, ElementVector) object>
  Number of elements: 63697
  Number of DOFs: 35007
  Size: 293515776 B

To PyVista: add the value of the field to the pyvista.UnstructuredGrid obtained through ansh.Mesh.to_pyvista().

Note: It is important to note that for creation of a PyVista object, the value of the field needs to be interpolated over the mesh points or cell centres. In this example, we do not have this problem since DOFs coincide with the mesh points. However, if the field definition involved element_type=ansh.element_types.TetP2(), a projection on to TetP1 is calculated internally to add values to the pyvista.UnstructuredGrid.

In [19]:

f.to_pyvista()

Out[19]:

UnstructuredGrid (0x7f84a2062b60)
  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:   9

Plot: plot the value of the vector field using pyvista. The convenience function creates a glyph map of the vector field with an optional colour_with to colour the glyphs using either x, y, or z components. Additionally, one can colour the glyphs using magnitude of the vector.

In [20]:

f.plot(colour_with="z", factor=3e-3, cmap="RdBu")

You can also plot the value of the field for specific subregions.

In [21]:

f.plot(subregions=["core"], colour_with="z", factor=3e-3, cmap="RdBu")

Set value: it is possible to set the value of the vector field as a callable. The callable is given the DOF locations and it is expected to output a two dimensional Numpy array of vector values corresponding to each DOF location.

In [22]:

def set_dof_val(dofs):
    return dofs

In [23]:

f.value = set_dof_val

In [24]:

f.plot(colour_with="z", factor=8e-2, cmap="RdBu")

Set value per subregion: it is possible to update the value for a given subregion.

In [25]:

f.set_subregion_value(name="shell", value=[0.0, 0.0, 0.0])

Note: set_subregion_value method also accepts a callable similar to value setter.

In [26]:

f.plot(colour_with="x", factor=8e-2, cmap="RdBu")

Vector field read write¶

Write vector field: at the moment, the scalar field can be written only as an HDF5 (.h5) or a VTK Unstructured Grid (.vtu) file. Please create an issue if you want to add support for other file formats.

Note: the .vtu file can be opened in ParaView for interactive visualisation.

In [27]:

f.to_file("vector.h5")

In [28]:

!h5glance --attrs vector.h5

vector.h5
├mesh
│ ├3 attributes:
│ │ ├about: 'one of the unit test meshes'
│ │ ├name: 'shell'
│ │ └scale: np.float64(1.0)
│ ├cells
│ │ ├0
│ │ │ └line  [int64: 78 × 2]
│ │ ├1
│ │ │ └tetra [int64: 63697 × 4]
│ │ └2
│ │   └triangle  [int64: 9058 × 3]
│ ├points    [float64: 11669 × 3]
│ └subregions
│   ├0
│   │ ├0 [int64: 78]
│   │ ├1 [int64: 0]
│   │ └2 [int64: 2820]
│   ├boundary
│   │ ├0 [int64: 0]
│   │ ├1 [int64: 0]
│   │ └2 [int64: 6238]
│   ├core
│   │ ├0 [int64: 0]
│   │ ├1 [int64: 22577]
│   │ └2 [int64: 0]
│   └shell
│     ├0 [int64: 0]
│     ├1 [int64: 41120]
│     └2 [int64: 0]
└field
  ├2 attributes:
  │ ├element_type: 'ElementTetP1'
  │ └type: 'VectorField'
  └value [float64: 11669 × 3]

In [29]:

f.to_file("vector.vtu")

Read vector field: one can read an HDF5 file created using the method above or a .vtu file.

In [30]:

ansh.VectorField.from_file("vector.h5")

Out[30]:

Ansh VectorField

Mesh: 11669 points, 72833 cells, 4 subregions
- 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
  - 0:
    
    line: 78
    tetra: 0
    triangle: 2820
  - boundary:
    
    line: 0
    tetra: 0
    triangle: 6238
  - core:
    
    line: 0
    tetra: 22577
    triangle: 0
  - shell:
    
    line: 0
    tetra: 41120
    triangle: 0
- Scale: 1.0
- Info: 2 items
  - about: one of the unit test meshes
  - name: shell
Element: ElementTetP1
Value: min=-1.8633949892534885, max=1.866760251091686
array([[0., 0., 0.], [0., 0., 0.], [0., 0., 0.], ..., [0., 0., 0.], [0., 0., 0.], [0., 0., 0.]], shape=(11669, 3))

In [31]:

ansh.VectorField.from_file("vector.vtu")

Out[31]:

Ansh VectorField

Mesh: 11669 points, 72833 cells, 4 subregions
- 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
  - boundary:
    
    line: 0
    tetra: 0
    triangle: 6238
  - shell:
    
    line: 0
    tetra: 41120
    triangle: 0
  - core:
    
    line: 0
    tetra: 22577
    triangle: 0
  - 0:
    
    line: 78
    tetra: 0
    triangle: 2820
- Scale: 1.0
- Info: 3 items
  - AnshScale: [1.]
  - name: ['shell']
  - about: ['one of the unit test meshes']
Element: ElementTetP1
Value: min=-1.8633949892534885, max=1.866760251091686
pyvista_ndarray([[0., 0., 0.], [0., 0., 0.], [0., 0., 0.], ..., [0., 0., 0.], [0., 0., 0.], [0., 0., 0.]], shape=(11669, 3))

Projection to P1 for non-linear elements with `.vtu`¶

As mentioned previously, if you save a field with a higher order element as a .vtu file, it will be projected to a P1 field due to VTK limitations.

In [32]:

ho_f = ansh.VectorField(f.mesh, ansh.element_types.TetP2(), [2., 2., 2.])
ho_f

Out[32]:

Ansh VectorField

Mesh: 11669 points, 72833 cells, 4 subregions
- 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
  - boundary:
    
    line: 0
    tetra: 0
    triangle: 6238
  - shell:
    
    line: 0
    tetra: 41120
    triangle: 0
  - core:
    
    line: 0
    tetra: 22577
    triangle: 0
  - 0:
    
    line: 78
    tetra: 0
    triangle: 2820
- Scale: 1.0
- Info: 2 items
  - name: shell
  - about: one of the unit test meshes
Element: ElementTetP2
Value: min=2.0, max=2.0
array([[2., 2., 2.], [2., 2., 2.], [2., 2., 2.], ..., [2., 2., 2.], [2., 2., 2.], [2., 2., 2.]], shape=(90153, 3))

In [33]:

ho_f.to_file("vector_ho.vtu")

In [34]:

ansh.VectorField.from_file("vector_ho.vtu")

Out[34]:

Ansh VectorField

Mesh: 11669 points, 72833 cells, 4 subregions
- 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
  - boundary:
    
    line: 0
    tetra: 0
    triangle: 6238
  - shell:
    
    line: 0
    tetra: 41120
    triangle: 0
  - core:
    
    line: 0
    tetra: 22577
    triangle: 0
  - 0:
    
    line: 78
    tetra: 0
    triangle: 2820
- Scale: 1.0
- Info: 3 items
  - AnshScale: [1.]
  - name: ['shell']
  - about: ['one of the unit test meshes']
Element: ElementTetP1
Value: min=1.9999999999999523, max=2.0000000000000346
pyvista_ndarray([[2., 2., 2.], [2., 2., 2.], [2., 2., 2.], ..., [2., 2., 2.], [2., 2., 2.], [2., 2., 2.]], shape=(11669, 3))

Note: dues to projection from P2 to P1 element type, there is a slight change in the min and max value of the field

Ansh for defining vector fields¶

Ansh Mesh

Ansh VectorField

Basic vector field properties¶

Ansh Mesh

Additional vector field properties and methods¶

Vector field read write¶

Ansh VectorField

Ansh VectorField

Projection to P1 for non-linear elements with .vtu¶

Ansh VectorField

Ansh VectorField

Projection to P1 for non-linear elements with `.vtu`¶