Polyfem Python Tutorial¶
This is a jupyter notebook. The "real" notebook can be found here.
Polyfem relies on 3 main objects:
Settingsthat contains the main settings such discretization order (e.g., $P_1$ or $P_2$), material parameters, formulation, etc.Problemthat describe the problem you want to solve, that is the boundary conditions and right-hand side. There are some predefined problems, such asDrivenCavity, or generic problems, such asGenericTensor.Solverthat is the actual FEM solver.
The usage of specific problems is indented for benchmarking, in general you want to use the GenericTensor for tensor-based PDEs (e.g., elasticity) or GenericScalar for scalar PDEs (e.g., Poisson).
A typical use of Polyfem is:
settings = polyfempy.Settings(
pde=polyfempy.PDEs.LinearElasticity, # or any other PDE
discr_order=2
)
# set necessary settings
# e.g. settings.discr_order = 2
problem = polyfempy.Problem() # or any other problem
# set problem related data
# e.g. problem.set_displacement(1, [0, 0], [True, False])
settings.problem = problem
#now we can create a solver and solve
solver = polyfempy.Solver()
solver.settings(settings)
solver.load_mesh_from_path(mesh_path, normalize_mesh=False)
solver.solve()
Note 1: for legacy reasons Polyfem always normalizes the mesh (i.e., rescale it to lay in the $[0,1]^d$ box, you can use normalize_mesh=False when loading to disable this feature.
Note 2: the solution $u(x)$ of a FEM solver are the coefficients $u_i$ you need to multiply the bases $\varphi_i(x)$ with: $$ u(x)=\sum u_i \varphi_i(x). $$ The coefficients $u_i$ are unrelated with the mesh vertices because of reordering of the nodes or high-order bases. For instance $P_2$ bases have additional nodes on the edges which do not exist in the mesh.
For this reason Polyfem uses a visualization mesh where the solution is sampled at the vertices. This mesh has two advantages:
- it solves the problem of nodes reordering and additional nodes in the same way
- it provides a "true" visualization for high order solution by densely sampling each element (a $P_2$ solution is a piecewise quadratic function which is visualized in a picewise linear fashion, thus the need of a dense element sampling).
To control the resolution of the visualization mesh use vismesh_rel_area while loading.
Examples¶
Some imports for plotting
import meshplot as mp
algebra
import numpy as np
and finally polyfempy
import polyfempy as pf
Utility¶
Creates a quad mesh of n_pts x n_pts in the form of a regular grid
def create_quad_mesh(n_pts):
extend = np.linspace(0,1,n_pts)
x, y = np.meshgrid(extend, extend, sparse=False, indexing='xy')
pts = np.column_stack((x.ravel(), y.ravel()))
faces = np.ndarray([(n_pts-1)**2, 4],dtype=np.int32)
index = 0
for i in range(n_pts-1):
for j in range(n_pts-1):
faces[index, :] = np.array([
j + i * n_pts,
j+1 + i * n_pts,
j+1 + (i+1) * n_pts,
j + (i+1) * n_pts
])
index = index + 1
return pts, faces
Plate hole¶
This is the python version of the plate with hole example explained here.
Set the mesh path
mesh_path = "plane_hole.obj"
create settings:
- Pick linear $P_1$ elements (if the mesh would be a quad it would be $Q_1$)
- We are use a linear material model
settings = pf.Settings(
discr_order=1,
pde=pf.PDEs.LinearElasticity
)
and choose Young's modulus and poisson ratio
settings.set_material_params("E", 210000)
settings.set_material_params("nu", 0.3)
Next we setup the problem
problem = pf.Problem()
sideset 1 has symetric boundary in $x$
problem.set_x_symmetric(1)
sideset 4 has symmetric boundary in $y$
problem.set_y_symmetric(4)
sideset 3 has a force of [100, 0] applied
problem.set_force(3, [100, 0])
fianally set the problem
settings.problem = problem
Solve! Note: we normalize the mesh to be in $[0,1]^2$
solver = pf.Solver()
solver.settings(settings)
solver.load_mesh_from_path(mesh_path, normalize_mesh=True)
solver.solve()
Get the solution
pts, tets, disp = solver.get_sampled_solution()
diplace the mesh
vertices = pts + disp
and get the stresses
mises, _ = solver.get_sampled_mises_avg()
finally plot with the above code
mp.plot(vertices, tets, mises, return_plot=True)
Note that we used get_sampled_mises_avg to get the Von Mises stresses because they depend on the Jacobian of the displacement which, becase we use $P_1$ elements, is piece-wise constant. To avoid that effect in get_sampled_mises_avg the mises are averaged per vertex weighted by the area of the triangles. If you want the "real" mises just call
mises = solver.get_sampled_mises()
mp.plot(vertices, tets, mises, return_plot=True)
Plate hole with High-Order Mesh¶
This is the same example as before, but we use wildmeshing to create a curved mesh.
import wildmeshing as wm
mesh_path = "plane_hole.svg"
v, f, nodes, F_nodes = wm.triangulate_svg(mesh_path, cut_outside=True)
Now we proceed as before
#create settings
settings = pf.Settings(
discr_order=1, #pick linear P_1 elements, even if the geometry is P_3
pde=pf.PDEs.LinearElasticity #Linear elasticity
)
#Material parameters
settings.set_material_params("E", 210000)
settings.set_material_params("nu", 0.3)
Next we setup the problem as before
problem = pf.Problem()
#sideset 1 is symmetric in x
problem.set_x_symmetric(1)
#sideset 4 is symmetric in y
problem.set_y_symmetric(4)
#sideset 3 has a force of [100, 0] applied
problem.set_force(3, [100, 0])
fianally set the problem
settings.problem = problem
Create the solver and load the high-order mesh, the only difference with respect to before
solver = pf.Solver()
solver.settings(settings)
solver.set_high_order_mesh(v, f, nodes, F_nodes, normalize_mesh=True)
And finally, solve!
solver.solve()
Get and plot the solution, same as before
pts, tets, disp = solver.get_sampled_solution()
#diplace the mesh
vertices = pts + disp
#get the stresses
mises, _ = solver.get_sampled_mises_avg()
#plot
mp.plot(vertices, tets, mises, return_plot=True)
Torsion¶
Non-linear example. We want to torque a 3D bar around the $z$ direction.
The example is really similar as the one just above.
Sets the mesh, create a solver, and load the mesh.
In this case the mesh has already the correct size. We also choose a coarse visualization mesh
mesh_path = "square_beam_h.HYBRID"
solver = pf.Solver()
solver.load_mesh_from_path(mesh_path, normalize_mesh=False, vismesh_rel_area=0.001)
We want to use the default sideset marking, top of the mesh is 5 and bottom is 6.
Let's verify. We first extract the sidesets: p are some point, t triangles, and s the sidesets from 1 to 6
p, t, s = solver.get_boundary_sidesets()
Now we can plot it
tmp = np.zeros_like(s)
tmp[s==5] = 1
tmp[s==6] = 1
mp.plot(p, t, tmp, return_plot=True)
Now we create the settings, as before.
Note: It is an hex-mesh so we are using $Q_1$.
Differently from before we want a non-linear material model: NeoHookean.
To avoid ambiguities in the rotation, we want to do 5 steps of incremental loading.
settings = pf.Settings(
discr_order=1,
pde=pf.PDEs.NonLinearElasticity,
nl_solver_rhs_steps=5
)
Choose Young's modulus and Poisson's ratio, as before
settings.set_material_params("E", 200)
settings.set_material_params("nu", 0.35)
Now we setup problem with fixed sideset (5), rotating sideset (6), ahlf a tour along the $z$-axis.
problem = pf.Torsion(
fixed_boundary=5, #sideset 5 is fixed
turning_boundary=6, #sideset 6 rotates
n_turns = 0.5, #by half a tour
axis_coordiante=2, #around the z-axis, 2
)
and set the problem and solve
settings.problem = problem
#solving!
solver.settings(settings)
solver.solve()
takes approx 1 min because it is a complicated non-linear problem!
Get solution and stesses as before
Since we want to show only the surface there is no need of getting the whole volume, so we set boundary_only to True
pts, tets, disp = solver.get_sampled_solution(boundary_only=True)
vertices = pts + disp
mises, _ = solver.get_sampled_mises_avg(boundary_only=True)
and plot the 3d result!
mp.plot(vertices, tets, mises, shading={"flat":True}, return_plot=True)
Fluid simulation¶
Create the mesh using the utility function
pts, faces = create_quad_mesh(50)
create settings, pick linear $Q_2$ elements for velocity and $Q_1$ for pressure and select stokes as material model.
settings = pf.Settings(
discr_order=2,
pressure_discr_order=1,
pde=pf.PDEs.Stokes
)
Set the viscosity of the fluid
settings.set_material_params("viscosity", 1)
The default solver do not support mixed formulation, we need to choose UmfPackLU
settings.set_advanced_option("solver_type", "Eigen::UmfPackLU")
We use the standard Driven Cavity problem
problem = pf.DrivenCavity()
we set the problem
settings.problem = problem
We create the solver and set the settings
solver = pf.Solver()
solver.settings(settings)
This time we are using pts and faces instead of loading from the disk
solver.set_mesh(pts, faces, vismesh_rel_area=0.001)
Solve!
solver.solve()
We now get the solution and the pressure
pts, tris, velocity = solver.get_sampled_solution()
and plot it!
each = 10
p = mp.plot(pts, tris, np.linalg.norm(velocity, axis=1), return_plot=True)
p.add_lines(pts[0:pts.shape[0]:each,:], pts[0:pts.shape[0]:each,:]+velocity[0:pts.shape[0]:each,:]/6)
p
Time dependent simulation¶
Create the mesh using the utility function
pts, faces = create_quad_mesh(50)
create settings, pick linear $Q_1$ elements, and a linear material model.
In this case we want to run a simulation with $t\in[0, 10]$ and have 50 time steps.
settings = pf.Settings(
discr_order=1,
pde=pf.PDEs.LinearElasticity,
tend=10,
time_steps=50
)
Choose Young's modulus and poisson ratio
settings.set_material_params("E", 1)
settings.set_material_params("nu", 0.3)
Next we setup the problem, this doesnt have any parameters. It will...
problem = pf.Gravity()
we set the problem
settings.problem = problem
We create the solver and set the settings
solver = pf.Solver()
solver.settings(settings)
This time we are using pts and faces instead of loading from the disk (For efficienty in the browser we select a coarse vis mesh)
solver.set_mesh(pts, faces, vismesh_rel_area=0.001)
Solve!
solver.solve()
Get the solution and the mises
pts = solver.get_sampled_points_frames()
tris = solver.get_sampled_connectivity_frames()
disp = solver.get_sampled_solution_frames()
mises = solver.get_sampled_mises_avg_frames()
def plot(frame, p=None):
return mp.plot(pts[frame]+disp[frame], tris[frame], mises[frame], return_plot=True, plot=p)
Before the animation, let's plot the solution some frames
plot(0)
plot(5)
plot(11)
Now we are ready to do the animation
p = plot(0)
@mp.interact(frame=(0, len(mises)))
def step(frame=0):
plot(frame, p)