Ansh for defining scalar 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 scalar 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]:
m
Out[5]:
Ansh Mesh
- Bounds: [-2.999, -2.998, -3.0] - [2.998, 2.997, 3.0]
- Points: 11669
-
- line: 78
- tetra: 63697
- triangle: 9058
-
- shell:
- line: 0
- tetra: 41120
- triangle: 0
- core:
- line: 0
- tetra: 22577
- triangle: 0
- boundary:
- line: 0
- tetra: 0
- triangle: 6238
- 0:
- line: 78
- tetra: 0
- triangle: 2820
- shell:
- Scale: 1.0
-
- name: shell
- about: one of the unit test meshes
In [6]:
f = ansh.ScalarField(mesh=m, element_type=ansh.element_types.TetP1(), value=42.0)
In [7]:
f
Out[7]:
Ansh ScalarField
-
- Bounds: [-2.999, -2.998, -3.0] - [2.998, 2.997, 3.0]
- Points: 11669
-
- line: 78
- tetra: 63697
- triangle: 9058
-
- shell:
- line: 0
- tetra: 41120
- triangle: 0
- core:
- line: 0
- tetra: 22577
- triangle: 0
- boundary:
- line: 0
- tetra: 0
- triangle: 6238
- 0:
- line: 78
- tetra: 0
- triangle: 2820
- shell:
- Scale: 1.0
-
- name: shell
- about: one of the unit test meshes
- Element: ElementTetP1
-
array([42., 42., 42., ..., 42., 42., 42.], shape=(11669,))
Basic scalar 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
-
- line: 78
- tetra: 63697
- triangle: 9058
-
- shell:
- line: 0
- tetra: 41120
- triangle: 0
- core:
- line: 0
- tetra: 22577
- triangle: 0
- boundary:
- line: 0
- tetra: 0
- triangle: 6238
- 0:
- line: 78
- tetra: 0
- triangle: 2820
- shell:
- Scale: 1.0
-
- name: shell
- about: one of the unit test meshes
In [9]:
f.mesh is m
Out[9]:
True
- Element type: the
skfemelement type used to approximate the field.
In [10]:
f.element_type
Out[10]:
<skfem.element.element_tet.element_tet_p1.ElementTetP1 at 0x7f63465f7d40>
Note:
ansh.element_typesis a convenience named tuple that exposes all theskfemelement types to Ansh. To see all the available element types, hitTabafteransh.element_types.in a jupyter notebook (or anipythonshell). The auto-complete should also work with popular Python LSPs in your preferred IDE.
- Value: the value of the scalar 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 [11]:
f.value
Out[11]:
array([42., 42., 42., ..., 42., 42., 42.], shape=(11669,))
In [12]:
f.mesh.points.shape[0]
Out[12]:
11669
Additional scalar field properties and methods¶
- Degree of freedom locations: the points at which the value of the scalar field approximation is allowed to change. The value of field is interpolated in rest of the geometry using the basis functions.
In [13]:
f.dof_locations
Out[13]:
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 scalar field approximations can be obtained by calling the field over the coordinates of one or more points inside the mesh.
In [14]:
f([0, 0, 0])
Out[14]:
array([42.])
In [15]:
f(np.array([[0, 2.0, -1.0], [0, 0, 0], [-2, -1, 1.9]]))
Out[15]:
array([42., 42., 42.])
- To
skfem: create askfem.DiscreteFieldfromansh.ScalaFieldfor further processing such as matrix assembly. It is also possible to just extractskfembasis.
In [16]:
f.to_skfem()
Out[16]:
<skfem DiscreteField object> Quadrature points per element: 4 Number of elements: 63697 Order: 0 Attributes: grad
In [17]:
f.skfem_basis
Out[17]:
<skfem CellBasis(MeshTet1, ElementTetP1) object> Number of elements: 63697 Number of DOFs: 11669 Size: 32612864 B
- To PyVista: add the value of the field to the
pyvista.UnstructuredGridobtained throughansh.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 toTetP1is calculated internally to add the values topyvista.UnstructuredGrid.
In [18]:
f.to_pyvista()
Out[18]:
UnstructuredGrid (0x7f633d652ec0) 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 scalar field using
pyvista. The convenience function creates a movable clip widget which can be utilised to inspect the value of the field within the geometry.
In [19]:
f.plot(cmap="RdBu")
It is also possible to plot the value of the field for specific subregions.
In [20]:
f.plot(subregions=["core"], cmap="RdBu")
- Set value: it is possible to set the value of the scalar field as a callable. The callable is given the DOF locations and it is expected to output a 1D Numpy array of values corresponding to each DOF location.
In [21]:
def norm_val(dofs):
return np.linalg.norm(dofs, axis=-1)
In [22]:
f.value = norm_val
In [23]:
f.plot(cmap="hot")
- Set value per subregion: it is possible to update the value for a given subregion.
In [24]:
f.set_subregion_value(name="boundary", value=0.0)
Note:
set_subregion_valuemethod also accepts callable similar tovaluesetter.
In [25]:
f.plot(cmap="hot")
Scalar field read write¶
- Write scalar 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
.vtufile can be opened in ParaView for interactive visualisation.
In [26]:
f.to_file("scalar.h5")
In [27]:
!h5glance --attrs scalar.h5
scalar.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: 'ScalarField' └value [float64: 11669]
In [28]:
f.to_file("scalar.vtu")
- Read scalar field: one can read an HDF5 file created using the method above or a
.vtufile.
In [29]:
ansh.ScalarField.from_file("scalar.h5")
Out[29]:
Ansh ScalarField
-
- Bounds: [-2.999, -2.998, -3.0] - [2.998, 2.997, 3.0]
- Points: 11669
-
- line: 78
- tetra: 63697
- triangle: 9058
-
- 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
- 0:
- Scale: 1.0
-
- about: one of the unit test meshes
- name: shell
- Element: ElementTetP1
-
array([2. , 2. , 0. , ..., 2.95880067, 2.67804299, 2.95484998], shape=(11669,))
In [30]:
ansh.ScalarField.from_file("scalar.vtu")
Out[30]:
Ansh ScalarField
-
- Bounds: [-2.999, -2.998, -3.0] - [2.998, 2.997, 3.0]
- Points: 11669
-
- line: 78
- tetra: 63697
- triangle: 9058
-
- shell:
- line: 0
- tetra: 41120
- triangle: 0
- core:
- line: 0
- tetra: 22577
- triangle: 0
- boundary:
- line: 0
- tetra: 0
- triangle: 6238
- 0:
- line: 78
- tetra: 0
- triangle: 2820
- shell:
- Scale: 1.0
-
- AnshScale: [1.]
- name: ['shell']
- about: ['one of the unit test meshes']
- Element: ElementTetP1
-
pyvista_ndarray([2. , 2. , 0. , ..., 2.95880067, 2.67804299, 2.95484998], shape=(11669,))
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 [31]:
ho_f = ansh.ScalarField(f.mesh, ansh.element_types.TetP2(), f.value[0])
ho_f
Out[31]:
Ansh ScalarField
-
- Bounds: [-2.999, -2.998, -3.0] - [2.998, 2.997, 3.0]
- Points: 11669
-
- line: 78
- tetra: 63697
- triangle: 9058
-
- shell:
- line: 0
- tetra: 41120
- triangle: 0
- core:
- line: 0
- tetra: 22577
- triangle: 0
- boundary:
- line: 0
- tetra: 0
- triangle: 6238
- 0:
- line: 78
- tetra: 0
- triangle: 2820
- shell:
- Scale: 1.0
-
- name: shell
- about: one of the unit test meshes
- Element: ElementTetP2
-
array([2., 2., 2., ..., 2., 2., 2.], shape=(90153,))
In [32]:
ho_f.to_file("scalar_ho.vtu")
In [33]:
ansh.ScalarField.from_file("scalar_ho.vtu")
Out[33]:
Ansh ScalarField
-
- Bounds: [-2.999, -2.998, -3.0] - [2.998, 2.997, 3.0]
- Points: 11669
-
- line: 78
- tetra: 63697
- triangle: 9058
-
- shell:
- line: 0
- tetra: 41120
- triangle: 0
- core:
- line: 0
- tetra: 22577
- triangle: 0
- boundary:
- line: 0
- tetra: 0
- triangle: 6238
- 0:
- line: 78
- tetra: 0
- triangle: 2820
- shell:
- Scale: 1.0
-
- AnshScale: [1.]
- name: ['shell']
- about: ['one of the unit test meshes']
- Element: ElementTetP1
-
pyvista_ndarray([2., 2., 2., ..., 2., 2., 2.], shape=(11669,))
Note: dues to projection from P2 to P1 element type, there is a slight change in the min and max value of the field