Notebook 2: Variables

We can add discrete “variables” (unknowns associated with the mesh points) to a mesh, assign values to them and build expressions that sympy can understand, manipulate and simplify.

This notebook introduces the concept of MeshVariables in Underworld3. These are both data containers and sympy symbolic objects. We show you how to inspect a meshVariable, set the data values in the MeshVariable and visualise them, and build expressions that sympy can understand, manipulate and simplify.

import underworld3 as uw
import numpy as np
import sympy

mesh = uw.meshing.uw.meshing.CubedSphere(
    radiusOuter=1.0,
    radiusInner=0.547,
    numElements=8,
    simplex=True,
    verbose=False,
)

x,y,z = mesh.CoordinateSystem.X
r,th,phi = mesh.CoordinateSystem.R
r_vec = mesh.CoordinateSystem.unit_e_0
th_vec =  mesh.CoordinateSystem.unit_e_1
phi_vec =  mesh.CoordinateSystem.unit_e_2

This example shows how we can add a scalar field with a single value associated with each mesh node, and a vector field which has quadratic interpolation (points at the nodes plus interpolating points along mesh edges).

# mesh variable example / test

scalar_var = uw.discretisation.MeshVariable(
    varname="Radius",
    mesh=mesh, 
    vtype = uw.VarType.SCALAR,
    varsymbol=r"r"
)

vector_var = uw.discretisation.MeshVariable(
    varname="Vertical_Vec",
    mesh=mesh, 
    degree=2, #quadratic interpolation
    vtype = uw.VarType.VECTOR,
    varsymbol=r"\mathbf{v}",
)

To set values of the variable, we first have to unlock it using the access context manager, and then we can evaluate a function at the coordinates appropriate to fill up each variable. The locking is a PETSc requirement which is used to make sure data remains synchronised when we run code in parallel.

with mesh.access(scalar_var, vector_var):
    scalar_var.data[:,0] = uw.function.evaluate(r, scalar_var.coords)
    vector_var.data[:,:] = uw.function.evaluate(r_vec, vector_var.coords)

Variables are like most underworld and PETSc objects - they can be examined using their view() method. The information that you will see is split into the underworld representation (listed under MeshVariable) and the PETSc representation (listed under FE Data which also includes the numerical values).

scalar_var.view()

Class: <class ‘underworld3.discretisation._MeshVariable’>

MeshVariable:

symbol: \({ \hspace{ 0.05pt } {r} }\)

shape: \((1, 1)\)

degree: \(1\)

continuous: True

type: SCALAR

FE Data:

PETSc field id: \(0\)

PETSc field name: Radius

array([[1.        ],
       [1.        ],
       [1.        ],
       ...,
       [0.88636051],
       [0.83309475],
       [0.79595533]])

# Visualise it / them

import pyvista as pv
import underworld3.visualisation as vis

pvmesh = vis.mesh_to_pv_mesh(mesh)
pvmesh.point_data["z"] = vis.scalar_fn_to_pv_points(pvmesh, mesh.CoordinateSystem.X[2])
pvmesh.point_data["r"] = vis.scalar_fn_to_pv_points(pvmesh, scalar_var.sym[0])
pvmesh.point_data["V"] = vis.vector_fn_to_pv_points(pvmesh, (1-scalar_var.sym[0]) * vector_var.sym)

if mesh.dim==3:
    pvmesh_clipped = pvmesh.clip( normal='z', crinkle=False,origin=(0.0,0.0,0.01))

# pvmesh.plot(show_edges=True, show_scalar_bar=False)

pl = pv.Plotter(window_size=(750, 750))

pl.add_mesh(pvmesh_clipped, 
            show_edges=True,
            scalars="z", 
            opacity=0.6,
            show_scalar_bar=False)

pl.add_arrows(pvmesh.points, 
              pvmesh.point_data["V"],
              mag=0.25,
              opacity=0.6,
              color="Black",
              show_scalar_bar=False)

pl.export_html("html5/echidna_plot.html")

from IPython.display import IFrame
IFrame(src="html5/echidna_plot.html", width=600, height=400)

Interactive Image: Spherical shell mesh cut away to show radial arrows with length decreasing away from the centre.

More information

The meshVariable code is described API docs here.

mesh._mark_faces_inside_and_out()

points = np.array([(0.0,0.75,0.0), (0.0,0.0,0.1)])

mesh.get_closest_local_cells(points)

Closest points (2,)

array([2195,   -1])

mesh.test_if_points_in_cells( points, np.array([2195,
       2829]))

array([ True, False])

mesh.faces_inner_control_points.shape
# (4, 7615, 3)

(4, 7615, 3)