# Polyfem Tutorial¶

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

Polyfem relies on 3 main objects:

`Settings`

that contains the main settings such discretization order (e.g., P_1 or P_2), material parameters, formulation, etc.`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`

.`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.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:
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 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()

[2019-10-01 10:18:11.046] [polyfem] [info] Loading mesh... [2019-10-01 10:18:11.047] [geogram] [info] Loading file plane_hole.obj... [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)