Skip to content

Ultimate Example

Imports

import numpy as np
import math
import igl
import meshplot as mp
import wildmeshing as wm
import polyfempy as pf

Read the Mesh

V, F = igl.read_triangle_mesh("data/mesh.obj")
mp.plot(V, F, shading={"wireframe": True})

Tetrahedralization with Wildmeshing

wm.tetrahedralize("data/mesh.obj", "out.mesh", mute_log=True)
True

Loading the Mesh in Polyfem

solver = pf.Solver()
solver.load_mesh_from_path("out.mesh")
[2020-05-25 21:26:14.381] [polyfem] [info] Loading mesh...
[2020-05-25 21:26:14.381] [geogram] [info] Loading file out.mesh...
[2020-05-25 21:26:14.413] [geogram] [info] (FP64) nb_v:3598 nb_e:0 nb_f:4696 nb_b:0 tri:1 dim:3
[2020-05-25 21:26:14.413] [geogram] [info]  nb_tets:13462
[2020-05-25 21:26:14.413] [geogram] [info] Attributes on vertices: point[3]
[2020-05-25 21:26:14.488] [polyfem] [info] mesh bb min [-92.075, -98.5274, -42.037], max [92.075, 28.6772, 42.037]
[2020-05-25 21:26:14.488] [polyfem] [info]  took 0.107315s

Advanced Sidesets

minn = np.min(V, axis=0)
maxx = np.max(V, axis=0)

center = (minn+maxx)/2

Function to select sidesets:

  • Sideset 2: z-component close to the bottom boundary minn[2] and larger than y of center
  • Sideset 3: z-component close to the bottom boundary minn[2] and smaller than y of center
def sideset(p):
    if abs(p[2] - minn[2]) < 0.5:
        if p[1] > center[1]:
            return 2
        else:
            return 3
    return 1

Loading and Visualizing Sidesets

solver.set_boundary_side_set_from_bary(sideset)
ps, ts, s = solver.get_boundary_sidesets()

col = np.zeros_like(s); col[s==2] = 2; col[s==3] = 3
mp.plot(ps, ts, col)

Setup Problem

settings = pf.Settings()
problem = pf.Problem()

settings.set_pde(pf.PDEs.LinearElasticity)

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


problem.set_displacement(2, [0, 0, 0])
problem.set_force(3, [0, 0.5, 0])

settings.set_problem(problem)

Solve!

solver.settings(settings)
solver.solve()
[2020-05-25 21:26:14.614] [polyfem] [info] simplex_count:  13462
[2020-05-25 21:26:14.615] [polyfem] [info] regular_count:   0
[2020-05-25 21:26:14.615] [polyfem] [info] regular_boundary_count:  0
[2020-05-25 21:26:14.615] [polyfem] [info] simple_singular_count:   0
[2020-05-25 21:26:14.615] [polyfem] [info] multi_singular_count:    0
[2020-05-25 21:26:14.615] [polyfem] [info] boundary_count:  0
[2020-05-25 21:26:14.615] [polyfem] [info] multi_singular_boundary_count:   0
[2020-05-25 21:26:14.615] [polyfem] [info] non_regular_count:   0
[2020-05-25 21:26:14.615] [polyfem] [info] non_regular_boundary_count:  0
[2020-05-25 21:26:14.615] [polyfem] [info] undefined_count:     0
[2020-05-25 21:26:14.615] [polyfem] [info] total count:  13462
[2020-05-25 21:26:14.615] [polyfem] [info] Building isoparametric basis...
[2020-05-25 21:26:14.650] [polyfem] [info] Computing polygonal basis...
[2020-05-25 21:26:14.651] [polyfem] [info]  took 1.6e-05s
[2020-05-25 21:26:14.652] [polyfem] [info] hmin: 0.989271
[2020-05-25 21:26:14.652] [polyfem] [info] hmax: 18.93
[2020-05-25 21:26:14.652] [polyfem] [info] havg: 6.19137
[2020-05-25 21:26:14.652] [polyfem] [info]  took 0.035892s
[2020-05-25 21:26:14.652] [polyfem] [info] flipped elements 0
[2020-05-25 21:26:14.652] [polyfem] [info] h: 18.93
[2020-05-25 21:26:14.652] [polyfem] [info] n bases: 3598
[2020-05-25 21:26:14.652] [polyfem] [info] n pressure bases: 0
[2020-05-25 21:26:14.652] [polyfem] [info] Assigning rhs...
[2020-05-25 21:26:14.655] [polyfem] [info]  took 0.003053s
[2020-05-25 21:26:14.655] [polyfem] [info] Assembling stiffness mat...
[2020-05-25 21:26:14.756] [polyfem] [info]  took 0.100833s
[2020-05-25 21:26:14.756] [polyfem] [info] sparsity: 381532/116510436
[2020-05-25 21:26:14.756] [polyfem] [info] Solving LinearElasticity with
[2020-05-25 21:26:14.756] [polyfem] [info] Hypre...

Visualization of Results

p, t, d = solver.get_sampled_solution()
m = np.linalg.norm(d, axis=1)

mp.plot(p+d, t, m)

Isolines

p_uni, indices, inverse = np.unique(p, return_index=True, return_inverse=True, axis=0)

t_uni = np.array([inverse[t[:, 0]], inverse[t[:, 1]], inverse[t[:, 2]], inverse[t[:, 3]]]).transpose()
d_uni = d[indices, :]
m_uni = m[indices]
adj, _ = igl.tet_tet_adjacency(t_uni)

igl_faces = [[0,1,2], [0,1,3], [1,2,3], [2,0,3]]
surf = []
for t in range(adj.shape[0]):
    for f in range(adj.shape[1]):
        face = igl_faces[f]
        if adj[t, f] == -1:
            surf += [[t_uni[t, face[0]], t_uni[t, face[1]], t_uni[t, face[2]]]]
surf = np.array(surf)
iso_p, iso_l = igl.isolines(p_uni+d_uni, surf, m_uni, 20)

Plot!

p = mp.plot(p_uni+d_uni, surf, m_uni, return_plot=True)
p.add_edges(iso_p, iso_l, shading={"line_color": "white"})