# Polyfem Tutorial¶

This is a jupyter notebook. The “real” notebook can be found here.

Polyfem relies on 3 main objects:

1. Settings that contains the main settings such discretization order (e.g., $P_1$ or $P_2$), material parameters, formulation, etc.
2. Problem that describe the problem you want to solve, that is the boundary conditions and right-hand side. There are some predefined problems, such as DrivenCavity, or generic problems, such as GenericTensor.
3. Solver that 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.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: 1. it solves the problem of nodes reordering and additional nodes in the same way 2. 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 finallypolyfempy

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.solve()

[2019-10-01 10:18:11.046] [polyfem] [info] Loading mesh...
[2019-10-01 10:18:11.056] [geogram] [info] (FP64) nb_v:8549 nb_e:0 nb_f:16797 nb_b:299 tri:1 dim:3
[2019-10-01 10:18:11.056] [geogram] [info] Attributes on vertices: point[3]
[2019-10-01 10:18:11.067] [polyfem] [info] mesh bb min [0, 0], max [1, 0.500001]
[2019-10-01 10:18:11.067] [polyfem] [info]  took 0.0201773s
[2019-10-01 10:18:11.068] [polyfem] [info] simplex_count:   16797
[2019-10-01 10:18:11.068] [polyfem] [info] regular_count:   0
[2019-10-01 10:18:11.068] [polyfem] [info] regular_boundary_count:  0
[2019-10-01 10:18:11.068] [polyfem] [info] simple_singular_count:   0
[2019-10-01 10:18:11.068] [polyfem] [info] multi_singular_count:    0
[2019-10-01 10:18:11.068] [polyfem] [info] boundary_count:  0
[2019-10-01 10:18:11.068] [polyfem] [info] multi_singular_boundary_count:   0
[2019-10-01 10:18:11.068] [polyfem] [info] non_regular_count:   0
[2019-10-01 10:18:11.068] [polyfem] [info] non_regular_boundary_count:  0
[2019-10-01 10:18:11.068] [polyfem] [info] undefined_count:     0
[2019-10-01 10:18:11.068] [polyfem] [info] total count:  16797
[2019-10-01 10:18:11.068] [polyfem] [info] Building isoparametric basis...
[2019-10-01 10:18:11.089] [polyfem] [info] Computing polygonal basis...
[2019-10-01 10:18:11.089] [polyfem] [info]  took 1.6827e-05s
[2019-10-01 10:18:11.089] [polyfem] [info] hmin: 0.004207
[2019-10-01 10:18:11.089] [polyfem] [info] hmax: 0.01875
[2019-10-01 10:18:11.089] [polyfem] [info] havg: 0.00831951
[2019-10-01 10:18:11.089] [polyfem] [info]  took 0.0209313s
[2019-10-01 10:18:11.089] [polyfem] [info] flipped elements 0
[2019-10-01 10:18:11.089] [polyfem] [info] h: 0.01875
[2019-10-01 10:18:11.089] [polyfem] [info] n bases: 8549
[2019-10-01 10:18:11.089] [polyfem] [info] n pressure bases: 0
[2019-10-01 10:18:11.089] [polyfem] [info] Assigning rhs...
[2019-10-01 10:18:11.091] [polyfem] [info]  took 0.00127016s
[2019-10-01 10:18:11.091] [polyfem] [info] Assembling stiffness mat...
[2019-10-01 10:18:11.124] [polyfem] [info]  took 0.0336872s
[2019-10-01 10:18:11.124] [polyfem] [info] sparsity: 236956/292341604
[2019-10-01 10:18:11.124] [polyfem] [info] Solving LinearElasticity with
[2019-10-01 10:18:11.124] [polyfem] [info] Hypre...


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)