import mfem.ser as mfem
import ipywidgets as widgets
from glvis import glvis
from math import *
element_types = {
"2d": ["TRIANGLE", "QUADRILATERAL"],
"3d": ["TETRAHEDRON", "HEXAHEDRON"]
}
basis_map = {
"H1": mfem.H1_FECollection,
"L2": mfem.L2_FECollection
}
class callback_coeff(mfem.PyCoefficient):
def EvalValue(self, x):
return self.callback(x)
class callback_vcoeff(mfem.VectorPyCoefficient):
def EvalValue(self, x):
return self.callback(x)
# these have to be globals, if they are gc'd the runtime crashes (segfault)
coeff = fec = fespace = None
def generate(shape: tuple,
element_type: str,
order: int,
mesh_order: int,
surf: bool,
basis: mfem.FiniteElementCollection,
func: callable = None,
transform: callable = None):
global fec, fespace, coeff
if element_type in element_types["2d"]: shape = shape[:2]
mesh = mfem.mesh.Mesh(*shape, element_type)
if surf:
mesh.SetCurvature(mesh_order, True, 3)
else:
mesh.SetCurvature(mesh_order)
if transform is not None:
mesh_coeff = callback_vcoeff(mesh.SpaceDimension())
mesh_coeff.callback = transform
mesh.Transform(mesh_coeff)
fec = basis(order, mesh.Dimension())
fespace = mfem.FiniteElementSpace(mesh, fec)
if func is None: return mesh
x = mfem.GridFunction(fespace)
coeff = callback_coeff()
coeff.callback = func
x.ProjectCoefficient(coeff)
return (mesh, x)
# Widget parts
nx = widgets.IntSlider(description="nx:", min=1, value=8, max=50, step=1, continuous_update=False)
ny = widgets.IntSlider(description="ny:", min=1, value=8, max=50, step=1, continuous_update=False)
nz = widgets.IntSlider(description="nz:", min=1, value=8, max=50, step=1, continuous_update=False)
order = widgets.BoundedIntText(description='FEM Order:', value=2, min=1, step=1, continuous_update=False, layout={"width": "130px"})
mesh_order = widgets.BoundedIntText(description='Mesh Order:', value=2, min=1, step=1, continuous_update=False, layout={"width": "130px"})
surf = widgets.Checkbox(description='Surface Mesh', value=False)
element_type = widgets.Dropdown(
value="TRIANGLE",
options=[e for l in element_types.values() for e in l],
description='Mesh Elements:',
style={'description_width': 'initial'},
layout={"width": "240px"}
)
basis = widgets.Dropdown(
# value=basis_map.values()[0],
options=basis_map.keys(),
description="FEM Basis:",
style={'description_width': 'initial'},
layout={"width": "130px"}
)
func = widgets.Textarea(
disabled=False,
layout={"height": "56px", "width": "240px"}
)
transform = widgets.Textarea(
disabled=False,
layout={"height": "56px", "width": "240px"}
)
g = glvis(mfem.mesh.Mesh(nx.value, ny.value, element_type.value))
def toggle_dim(event=None):
if element_type.value in element_types["2d"]:
nz.disabled = True
func.placeholder = "FEM Function, e.g. x+y"
if surf == True:
transform.placeholder = "Mesh Transformation, e.g. (x,y,xy)"
else:
transform.placeholder = "Mesh Transformation, e.g. (x+y,x-y)"
else:
nz.disabled = False
func.placeholder = "FEM Function, e.g. x+y+z"
transform.placeholder = "Mesh Transformation, e.g. (x+y,x-y,z)"
# Setup handlers
def show(event=None):
# GridFunction callback
if func.value != "":
def gfcb(x):
local = {
"x": x[0],
"y": x[1],
"z": x[2] if len(x) > 2 else 0
}
return eval(func.value, None, local)
else:
gfcb = None
# Mesh callback
if transform.value != "":
def mcb(x):
local = {
"x": x[0],
"y": x[1],
"z": x[2] if (len(x) > 2 or surf == True) else 0
}
return eval(transform.value, None, local)
else:
mcb = None
g.plot(generate(
shape=(nx.value, ny.value, nz.value),
element_type=element_type.value,
order=order.value,
mesh_order=mesh_order.value,
surf=surf.value,
basis=basis_map[basis.value],
func=gfcb,
transform=mcb
))
button = widgets.Button(description="Update")
button.on_click(show)
# Figure updates on any change to nx or ny
nx.observe(show, names="value")
ny.observe(show, names="value")
nz.observe(show, names="value")
order.observe(show, names="value")
mesh_order.observe(show, names="value")
surf.observe(show, names="value")
element_type.observe(show, names="value")
basis.observe(show, names="value")
element_type.observe(toggle_dim, names="value")
centered = widgets.Layout(display='inline-flex',
align_items='center',
justify_content="center")
# Initial state
toggle_dim()
# Build widget
widgets.HBox([
g,
widgets.VBox([nx, ny, nz,
element_type,
mesh_order,
surf, transform,
basis, order, func,
button], layout=centered)
])
# Sample 2D Mesh Transformation
def spiral(x, y):
r = 1+x+4*y
phi = 4*pi*y
return (r*cos(phi), r*sin(phi))
# Sample 2D Surface Mesh Transformation
def breather(x, y):
u = 13.2*(2*x-1)
v = 37.4*(2*y-1)
b = 0.4
r = 1 - b*b
w = sqrt(r)
denom = b*((w*cosh(b*u))*(w*cosh(b*u)) + (b*sin(w*v))*(b*sin(w*v)))
return(-u + (2*r*cosh(b*u)*sinh(b*u)) / denom,
(2*w*cosh(b*u)*(-(w*cos(v)*cos(w*v)) - sin(v)*sin(w*v))) / denom,
(2*w*cosh(b*u)*(-(w*sin(v)*cos(w*v)) + cos(v)*sin(w*v))) / denom)
# Sample 3D Mesh Transformation
def toroid(x, y, z):
ns = 3 # Number of shifts
R = 1.0 # Major radius of the torus
r = 0.4 # Minor radius of the torus
phi = 2.0 * pi * z;
theta = phi * ns / 4;
u = sqrt(2) * (y - 0.5) * r;
v = sqrt(2) * (x - 0.5) * r;
return((R + u * cos(theta) + v * sin(theta)) * cos(phi),
(R + u * cos(theta) + v * sin(theta)) * sin(phi),
v * cos(theta) - u * sin(theta))
# Sample 2D Function
def wave2d(x,y):
r = sqrt((x-0.5)**2+(y-0.5)**2)
return cos(8*pi*r)*(1-r)**2
# Sample 3D Function
def wave3d(x, y, z):
r = sqrt((x-0.5)**2+(y-0.5)**2+(z-0.5)**2)
return cos(4*pi*r)*(1-r)