Skip to content

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()
# set necessary settings
# e.g. settings.discr_order = 2

settings.set_pde(polyfempy.PDEs.LinearElasticity) # or any other PDE

problem = polyfempy.Problem() # or any other problem
# set problem related data
# e.g. problem.set_displacement(1, [0, 0], [True, False])

settings.set_problem(problem)

#now we can create a solver and solve
solver = polyfempy.Solver()

solver.settings(settings)
solver.load_mesh_from_path(mesh_path)

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 setting.normalize_mesh = False 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 settings.vismesh_rel_area.

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

settings = pf.Settings()

pick linear P_1 elements. If the mesh would be a quad it would be Q_1

settings.discr_order = 1

normalize the mesh to be in [0,1]^2

settings.normalize_mesh = True

and choose Young’s modulus and poisson ratio

settings.set_material_params("E", 210000)
settings.set_material_params("nu", 0.3)

We are use a linear material model

settings.set_pde(pf.PDEs.LinearElasticity)

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.set_problem(problem)

Solve!

solver = pf.Solver()

solver.settings(settings)
solver.load_mesh_from_path(mesh_path)

solver.solve()
[2019-07-16 15:29:45.119] [polyfem] [info] Loading mesh...
[2019-07-16 15:29:45.119] [geogram] [info] Loading file plane_hole.obj...
[2019-07-16 15:29:45.131] [geogram] [info] (FP64) nb_v:8549 nb_e:0 nb_f:16797 nb_b:299 tri:1 dim:3
[2019-07-16 15:29:45.131] [geogram] [info] Attributes on vertices: point[3]
[2019-07-16 15:29:45.140] [polyfem] [info] mesh bb min [0.500183, 0.5], max [100.5, 50.5]
[2019-07-16 15:29:45.141] [polyfem] [info]  took 0.0210737s
[2019-07-16 15:29:45.142] [polyfem] [info] simplex_count:   16797
[2019-07-16 15:29:45.142] [polyfem] [info] regular_count:   0
[2019-07-16 15:29:45.142] [polyfem] [info] regular_boundary_count:  0
[2019-07-16 15:29:45.142] [polyfem] [info] simple_singular_count:   0
[2019-07-16 15:29:45.142] [polyfem] [info] multi_singular_count:    0
[2019-07-16 15:29:45.142] [polyfem] [info] boundary_count:  0
[2019-07-16 15:29:45.142] [polyfem] [info] multi_singular_boundary_count:   0
[2019-07-16 15:29:45.142] [polyfem] [info] non_regular_count:   0
[2019-07-16 15:29:45.142] [polyfem] [info] non_regular_boundary_count:  0
[2019-07-16 15:29:45.142] [polyfem] [info] undefined_count:     0
[2019-07-16 15:29:45.142] [polyfem] [info] total count:  16797
[2019-07-16 15:29:45.142] [polyfem] [info] Building isoparametric basis...
[2019-07-16 15:29:45.163] [polyfem] [info] Computing polygonal basis...
[2019-07-16 15:29:45.163] [polyfem] [info]  took 1.4877e-05s
[2019-07-16 15:29:45.164] [polyfem] [info] hmin: 0.420699
[2019-07-16 15:29:45.164] [polyfem] [info] hmax: 1.875
[2019-07-16 15:29:45.164] [polyfem] [info] havg: 0.831949
[2019-07-16 15:29:45.164] [polyfem] [info]  took 0.020509s
[2019-07-16 15:29:45.164] [polyfem] [info] flipped elements 0
[2019-07-16 15:29:45.164] [polyfem] [info] h: 1.875
[2019-07-16 15:29:45.164] [polyfem] [info] n bases: 8549
[2019-07-16 15:29:45.164] [polyfem] [info] n pressure bases: 0
[2019-07-16 15:29:45.164] [polyfem] [info] Assigning rhs...
[2019-07-16 15:29:45.164] [polyfem] [info]  took 0.00070312s
[2019-07-16 15:29:45.164] [polyfem] [info] Assembling stiffness mat...
[2019-07-16 15:29:45.197] [polyfem] [info]  took 0.0324238s
[2019-07-16 15:29:45.197] [polyfem] [info] sparsity: 236956/292341604
[2019-07-16 15:29:45.197] [polyfem] [info] Solving LinearElasticity with
[2019-07-16 15:29:45.197] [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, shading={"wireframe": False}, return_plot=True)
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(50.530025…

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, shading={"wireframe": False}, return_plot=True)
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(50.530025…

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, mute_log=True)
1/1 n sub paths: 5

Now we proceed as before

#create settings
settings = pf.Settings()

#pick linear P_1 elements, even if the geometry is P_3
settings.discr_order = 1

#Normalize
settings.normalize_mesh = True

#Material parameters
settings.set_material_params("E", 210000)
settings.set_material_params("nu", 0.3)

#Linear elasticity
settings.set_pde(pf.PDEs.LinearElasticity)

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.set_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)
[2019-07-16 15:29:47.764] [polyfem] [info] Loading mesh...
[2019-07-16 15:29:47.764] [polyfem] [info] mesh bb min [0, 0], max [1, 0.5]
[2019-07-16 15:29:47.764] [polyfem] [info]  took 3.5333e-05s

And finally, solve!

solver.solve()
[2019-07-16 15:29:47.776] [polyfem] [info] simplex_count:   314
[2019-07-16 15:29:47.776] [polyfem] [info] regular_count:   0
[2019-07-16 15:29:47.776] [polyfem] [info] regular_boundary_count:  0
[2019-07-16 15:29:47.776] [polyfem] [info] simple_singular_count:   0
[2019-07-16 15:29:47.776] [polyfem] [info] multi_singular_count:    0
[2019-07-16 15:29:47.776] [polyfem] [info] boundary_count:  0
[2019-07-16 15:29:47.776] [polyfem] [info] multi_singular_boundary_count:   0
[2019-07-16 15:29:47.776] [polyfem] [info] non_regular_count:   0
[2019-07-16 15:29:47.776] [polyfem] [info] non_regular_boundary_count:  0
[2019-07-16 15:29:47.776] [polyfem] [info] undefined_count:     0
[2019-07-16 15:29:47.776] [polyfem] [info] total count:  314
[2019-07-16 15:29:47.776] [polyfem] [info] Building not isoparametric basis...
[2019-07-16 15:29:47.777] [polyfem] [info] Computing polygonal basis...
[2019-07-16 15:29:47.777] [polyfem] [info]  took 6.363e-06s
[2019-07-16 15:29:47.777] [polyfem] [info] hmin: 0.0341246
[2019-07-16 15:29:47.777] [polyfem] [info] hmax: 0.0904418
[2019-07-16 15:29:47.777] [polyfem] [info] havg: 0.0597291
[2019-07-16 15:29:47.777] [polyfem] [info]  took 0.000643824s
[2019-07-16 15:29:47.777] [polyfem] [info] flipped elements 0
[2019-07-16 15:29:47.777] [polyfem] [info] h: 0.0904418
[2019-07-16 15:29:47.777] [polyfem] [info] n bases: 183
[2019-07-16 15:29:47.777] [polyfem] [info] n pressure bases: 0
[2019-07-16 15:29:47.777] [polyfem] [info] Assigning rhs...
[2019-07-16 15:29:47.777] [polyfem] [info]  took 0.000175593s
[2019-07-16 15:29:47.777] [polyfem] [info] Assembling stiffness mat...
[2019-07-16 15:29:47.778] [polyfem] [info]  took 0.00112067s
[2019-07-16 15:29:47.778] [polyfem] [info] sparsity: 4700/133956
[2019-07-16 15:29:47.778] [polyfem] [info] Solving LinearElasticity with
[2019-07-16 15:29:47.778] [polyfem] [info] Hypre...

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, shading={"wireframe": False}, return_plot=True)
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(0.5002871…

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)
[2019-07-16 15:29:47.929] [polyfem] [info] Loading mesh...
[2019-07-16 15:29:47.953] [polyfem] [info] mesh bb min [-10, -10, 0], max [10, 10, 100]
[2019-07-16 15:29:47.954] [polyfem] [info]  took 0.0248576s

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, shading={"wireframe": False}, return_plot=True)
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(0.0, 0.0,…

Now we create the settings, as before

settings = pf.Settings()

It is an hex-mesh so we are using Q_1

settings.discr_order = 1

Choose Young’s modulus and Poisson’s ratio, as before

settings.set_material_params("E", 200)
settings.set_material_params("nu", 0.35)

Differently from before we want a non-linear material model: NeoHookean

settings.set_pde(pf.PDEs.NonLinearElasticity)

and we want to do 5 steps of incremental loading to avoid ambiguities in the rotation

settings.nl_solver_rhs_steps = 5

Now we setup problem with fixed sideset, rotating an number of tours

problem = pf.Torsion()

sideset 5 is fixed

problem.fixed_boundary = 5

sideset 6 rotates

problem.turning_boundary = 6

around the z-axis (2)

problem.axis_coordiante = 2

by half a tour

problem.n_turns = 0.5

and set the problem and solve

settings.set_problem(problem)

#solving!
solver.settings(settings)

solver.solve()
[2019-07-16 15:29:48.101] [polyfem] [info] simplex_count:   0
[2019-07-16 15:29:48.101] [polyfem] [info] regular_count:   108
[2019-07-16 15:29:48.101] [polyfem] [info] regular_boundary_count:  356
[2019-07-16 15:29:48.101] [polyfem] [info] simple_singular_count:   0
[2019-07-16 15:29:48.101] [polyfem] [info] multi_singular_count:    0
[2019-07-16 15:29:48.101] [polyfem] [info] boundary_count:  0
[2019-07-16 15:29:48.101] [polyfem] [info] multi_singular_boundary_count:   0
[2019-07-16 15:29:48.101] [polyfem] [info] non_regular_count:   0
[2019-07-16 15:29:48.101] [polyfem] [info] non_regular_boundary_count:  0
[2019-07-16 15:29:48.101] [polyfem] [info] undefined_count:     0
[2019-07-16 15:29:48.101] [polyfem] [info] total count:  464
[2019-07-16 15:29:48.101] [polyfem] [info] Building isoparametric basis...
[2019-07-16 15:29:48.107] [polyfem] [info] Computing polygonal basis...
[2019-07-16 15:29:48.107] [polyfem] [info]  took 1.3212e-05s
[2019-07-16 15:29:48.107] [polyfem] [info] hmin: 3.4482
[2019-07-16 15:29:48.107] [polyfem] [info] hmax: 5
[2019-07-16 15:29:48.107] [polyfem] [info] havg: 4.41558
[2019-07-16 15:29:48.107] [polyfem] [info]  took 0.00593987s
[2019-07-16 15:29:48.107] [polyfem] [info] flipped elements 0
[2019-07-16 15:29:48.107] [polyfem] [info] h: 5
[2019-07-16 15:29:48.107] [polyfem] [info] n bases: 750
[2019-07-16 15:29:48.107] [polyfem] [info] n pressure bases: 0
[2019-07-16 15:29:48.107] [polyfem] [info] Assigning rhs...
[2019-07-16 15:29:48.108] [polyfem] [info]  took 0.000998475s
[2019-07-16 15:29:48.108] [polyfem] [info] Assembling stiffness mat...
[2019-07-16 15:29:48.108] [polyfem] [info]  took 3.524e-06s
[2019-07-16 15:29:48.108] [polyfem] [info] sparsity: 0/0
[2019-07-16 15:29:48.108] [polyfem] [info] Solving NeoHookean with
[2019-07-16 15:29:48.108] [polyfem] [info] t: 0.2 prev: 0 step: 0.2
[2019-07-16 15:29:55.892] [polyfem] [info] t: 0.4 prev: 0.2 step: 0.2
[2019-07-16 15:30:04.704] [polyfem] [info] t: 0.6 prev: 0.4 step: 0.2
[2019-07-16 15:30:12.333] [polyfem] [info] t: 0.8 prev: 0.6 step: 0.2
[2019-07-16 15:30:20.148] [polyfem] [info] t: 1 prev: 0.8 step: 0.2

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={"wireframe": False, "flat":True}, return_plot=True)
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(0.0, 0.0,…

Fluid simulation

Create the mesh using the utility function

pts, faces = create_quad_mesh(50)

create settings

settings = pf.Settings()

settings.vismesh_rel_area = 0.001

pick linear Q_2 elements for velocity and Q_1 for pressure

settings.discr_order = 2
settings.pressure_discr_order = 1

Set the viscosity of the fluid

settings.set_material_params("viscosity", 1)

We select stokes as material model

settings.set_pde(pf.PDEs.Stokes)

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.set_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)
[2019-07-16 15:30:28.385] [polyfem] [info] Loading mesh...
[2019-07-16 15:30:28.385] [polyfem] [info] mesh bb min [0, 0], max [1, 1]
[2019-07-16 15:30:28.385] [polyfem] [info]  took 0.000239624s

Solve!

solver.solve()
[2019-07-16 15:30:28.393] [polyfem] [info] simplex_count:   0
[2019-07-16 15:30:28.393] [polyfem] [info] regular_count:   2209
[2019-07-16 15:30:28.393] [polyfem] [info] regular_boundary_count:  192
[2019-07-16 15:30:28.393] [polyfem] [info] simple_singular_count:   0
[2019-07-16 15:30:28.393] [polyfem] [info] multi_singular_count:    0
[2019-07-16 15:30:28.393] [polyfem] [info] boundary_count:  0
[2019-07-16 15:30:28.393] [polyfem] [info] multi_singular_boundary_count:   0
[2019-07-16 15:30:28.393] [polyfem] [info] non_regular_count:   0
[2019-07-16 15:30:28.393] [polyfem] [info] non_regular_boundary_count:  0
[2019-07-16 15:30:28.393] [polyfem] [info] undefined_count:     0
[2019-07-16 15:30:28.393] [polyfem] [info] total count:  2401
[2019-07-16 15:30:28.393] [polyfem] [info] Building not isoparametric basis...
[2019-07-16 15:30:28.408] [polyfem] [info] Computing polygonal basis...
[2019-07-16 15:30:28.408] [polyfem] [info]  took 1.4753e-05s
[2019-07-16 15:30:28.409] [polyfem] [info] hmin: 0.0204082
[2019-07-16 15:30:28.409] [polyfem] [info] hmax: 0.0204082
[2019-07-16 15:30:28.409] [polyfem] [info] havg: 0.0204082
[2019-07-16 15:30:28.409] [polyfem] [info]  took 0.0151766s
[2019-07-16 15:30:28.409] [polyfem] [info] flipped elements 0
[2019-07-16 15:30:28.409] [polyfem] [info] h: 0.0204082
[2019-07-16 15:30:28.409] [polyfem] [info] n bases: 9801
[2019-07-16 15:30:28.409] [polyfem] [info] n pressure bases: 2500
[2019-07-16 15:30:28.409] [polyfem] [info] Assigning rhs...
[2019-07-16 15:30:28.410] [polyfem] [info]  took 0.00184905s
[2019-07-16 15:30:28.410] [polyfem] [info] Assembling stiffness mat...
[2019-07-16 15:30:28.478] [polyfem] [info]  took 0.0677301s
[2019-07-16 15:30:28.478] [polyfem] [info] sparsity: 549568/488498404
[2019-07-16 15:30:28.478] [polyfem] [info] Solving Stokes with
[2019-07-16 15:30:28.478] [polyfem] [info] N5Eigen9UmfPackLUINS_12SparseMatrixIdLi0EiEEEE...

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), shading={"wireframe": False}, return_plot=True)
p.add_lines(pts[0:pts.shape[0]:each,:], pts[0:pts.shape[0]:each,:]+velocity[0:pts.shape[0]:each,:])
p
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(0.5, 0.5,…

Time dependent simulation

Create the mesh using the utility function

pts, faces = create_quad_mesh(50)

create settings

settings = pf.Settings()

pick linear Q_1 elements.

settings.discr_order = 1

and choose Young’s modulus and poisson ratio

settings.set_material_params("E", 1)
settings.set_material_params("nu", 0.3)

We are use a linear material model

settings.set_pde(pf.PDEs.LinearElasticity)

For efficienty in the browser we select a coarse vis mesh

settings.vismesh_rel_area = 0.001

We simulate from 0 to 10s and 50 steps.

settings.tend = 10
settings.time_steps = 50

Next we setup the problem, this doesnt have any parameters. It will…

problem = pf.Gravity()

we set the problem

settings.set_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)
[2019-07-16 15:30:28.997] [polyfem] [info] Loading mesh...
[2019-07-16 15:30:28.997] [polyfem] [info] mesh bb min [0, 0], max [1, 1]

Solve!

solver.solve()
[2019-07-16 15:30:29.006] [polyfem] [info] simplex_count:   0
[2019-07-16 15:30:29.006] [polyfem] [info] regular_count:   2209
[2019-07-16 15:30:29.006] [polyfem] [info] regular_boundary_count:  192
[2019-07-16 15:30:29.006] [polyfem] [info] simple_singular_count:   0
[2019-07-16 15:30:29.006] [polyfem] [info] multi_singular_count:    0
[2019-07-16 15:30:29.006] [polyfem] [info] boundary_count:  0
[2019-07-16 15:30:29.006] [polyfem] [info] multi_singular_boundary_count:   0
[2019-07-16 15:30:29.006] [polyfem] [info] non_regular_count:   0
[2019-07-16 15:30:29.006] [polyfem] [info] non_regular_boundary_count:  0
[2019-07-16 15:30:29.006] [polyfem] [info] undefined_count:     0
[2019-07-16 15:30:29.006] [polyfem] [info] total count:  2401
[2019-07-16 15:30:29.006] [polyfem] [info] Building isoparametric basis...
[2019-07-16 15:30:29.009] [polyfem] [info] Computing polygonal basis...
[2019-07-16 15:30:29.009] [polyfem] [info]  took 1.093e-05s
[2019-07-16 15:30:29.009] [polyfem] [info] hmin: 0.0204082
[2019-07-16 15:30:29.009] [polyfem] [info] hmax: 0.0204082
[2019-07-16 15:30:29.009] [polyfem] [info] havg: 0.0204082
[2019-07-16 15:30:29.009] [polyfem] [info]  took 0.00349494s
[2019-07-16 15:30:29.009] [polyfem] [info] flipped elements 0
[2019-07-16 15:30:29.009] [polyfem] [info] h: 0.0204082
[2019-07-16 15:30:29.009] [polyfem] [info] n bases: 2500
[2019-07-16 15:30:29.009] [polyfem] [info] n pressure bases: 0
[2019-07-16 15:30:29.009] [polyfem] [info] Assigning rhs...
[2019-07-16 15:30:29.018] [polyfem] [info]  took 0.00823314s
[2019-07-16 15:30:29.018] [polyfem] [info] Assembling stiffness mat...
[2019-07-16 15:30:29.036] [polyfem] [info]  took 0.0186862s
[2019-07-16 15:30:29.036] [polyfem] [info] sparsity: 87616/25000000
[2019-07-16 15:30:29.036] [polyfem] [info] Solving LinearElasticity with
[2019-07-16 15:30:29.044] [polyfem] [info] Hypre...
[2019-07-16 15:30:29.206] [polyfem] [info] 1/50
[2019-07-16 15:30:29.298] [polyfem] [info] 2/50
[2019-07-16 15:30:29.386] [polyfem] [info] 3/50
[2019-07-16 15:30:29.477] [polyfem] [info] 4/50
[2019-07-16 15:30:29.573] [polyfem] [info] 5/50
[2019-07-16 15:30:29.671] [polyfem] [info] 6/50
[2019-07-16 15:30:29.767] [polyfem] [info] 7/50
[2019-07-16 15:30:29.857] [polyfem] [info] 8/50
[2019-07-16 15:30:29.948] [polyfem] [info] 9/50
[2019-07-16 15:30:30.038] [polyfem] [info] 10/50
[2019-07-16 15:30:30.129] [polyfem] [info] 11/50
[2019-07-16 15:30:30.227] [polyfem] [info] 12/50
[2019-07-16 15:30:30.324] [polyfem] [info] 13/50
[2019-07-16 15:30:30.409] [polyfem] [info] 14/50
[2019-07-16 15:30:30.498] [polyfem] [info] 15/50
[2019-07-16 15:30:30.596] [polyfem] [info] 16/50
[2019-07-16 15:30:30.687] [polyfem] [info] 17/50
[2019-07-16 15:30:30.778] [polyfem] [info] 18/50
[2019-07-16 15:30:30.868] [polyfem] [info] 19/50
[2019-07-16 15:30:30.959] [polyfem] [info] 20/50
[2019-07-16 15:30:31.054] [polyfem] [info] 21/50
[2019-07-16 15:30:31.147] [polyfem] [info] 22/50
[2019-07-16 15:30:31.237] [polyfem] [info] 23/50
[2019-07-16 15:30:31.322] [polyfem] [info] 24/50
[2019-07-16 15:30:31.419] [polyfem] [info] 25/50
[2019-07-16 15:30:31.513] [polyfem] [info] 26/50
[2019-07-16 15:30:31.603] [polyfem] [info] 27/50
[2019-07-16 15:30:31.692] [polyfem] [info] 28/50
[2019-07-16 15:30:31.825] [polyfem] [info] 29/50
[2019-07-16 15:30:31.925] [polyfem] [info] 30/50
[2019-07-16 15:30:32.014] [polyfem] [info] 31/50
[2019-07-16 15:30:32.113] [polyfem] [info] 32/50
[2019-07-16 15:30:32.213] [polyfem] [info] 33/50
[2019-07-16 15:30:32.301] [polyfem] [info] 34/50
[2019-07-16 15:30:32.386] [polyfem] [info] 35/50
[2019-07-16 15:30:32.473] [polyfem] [info] 36/50
[2019-07-16 15:30:32.562] [polyfem] [info] 37/50
[2019-07-16 15:30:32.651] [polyfem] [info] 38/50
[2019-07-16 15:30:32.756] [polyfem] [info] 39/50
[2019-07-16 15:30:32.873] [polyfem] [info] 40/50
[2019-07-16 15:30:32.980] [polyfem] [info] 41/50
[2019-07-16 15:30:33.094] [polyfem] [info] 42/50
[2019-07-16 15:30:33.208] [polyfem] [info] 43/50
[2019-07-16 15:30:33.320] [polyfem] [info] 44/50
[2019-07-16 15:30:33.419] [polyfem] [info] 45/50
[2019-07-16 15:30:33.511] [polyfem] [info] 46/50
[2019-07-16 15:30:33.599] [polyfem] [info] 47/50
[2019-07-16 15:30:33.688] [polyfem] [info] 48/50
[2019-07-16 15:30:33.775] [polyfem] [info] 49/50
[2019-07-16 15:30:33.863] [polyfem] [info] 50/50

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

Before the animation, let’s plot the solution at frame 25

frame = 25

mp.plot(pts[frame]+disp[frame], tris[frame], mises[frame], shading={"wireframe": False}, return_plot=True)
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(0.5000000…

Now we are ready to do the animation

p = mp.plot(pts[frame]+disp[frame], tris[frame], mises[frame], shading={"wireframe": False}, return_plot=True)

@mp.interact(frame=(0, len(mises)))
def step(frame=0):
#     p.update_object(vertices=pts[frame]+disp[frame]) 
    mp.plot(pts[frame]+disp[frame], tris[frame], mises[frame], shading={"wireframe": False}, plot=p)
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(0.5000000…



interactive(children=(IntSlider(value=0, description='frame', max=51), Output()), _dom_classes=('widget-intera…