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("mesh.obj") mp.plot(V, F, shading={"wireframe": True})
Tetrahedralization with Wildmeshing¶
wm.tetrahedralize("mesh.obj", "out.mesh", mute_log=True)
Loading the Mesh in Polyfem¶
solver = pf.Solver() solver.load_mesh_from_path("out.mesh")
[2019-11-29 15:44:36.036] [polyfem] [info] Loading mesh... [2019-11-29 15:44:36.038] [geogram] [info] Loading file out.mesh... [2019-11-29 15:44:36.087] [geogram] [info] (FP64) nb_v:3956 nb_e:0 nb_f:4858 nb_b:0 tri:1 dim:3 [2019-11-29 15:44:36.087] [geogram] [info] nb_tets:15402 [2019-11-29 15:44:36.087] [geogram] [info] Attributes on vertices: point[3] [2019-11-29 15:44:36.331] [polyfem] [info] mesh bb min [-92.075, -98.5272, -42.037], max [92.075, 28.6836, 42.037] [2019-11-29 15:44:36.331] [polyfem] [info] took 0.294995s
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 a 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()
[2019-11-29 15:44:36.511] [polyfem] [info] simplex_count: 15402 [2019-11-29 15:44:36.511] [polyfem] [info] regular_count: 0 [2019-11-29 15:44:36.511] [polyfem] [info] regular_boundary_count: 0 [2019-11-29 15:44:36.511] [polyfem] [info] simple_singular_count: 0 [2019-11-29 15:44:36.511] [polyfem] [info] multi_singular_count: 0 [2019-11-29 15:44:36.511] [polyfem] [info] boundary_count: 0 [2019-11-29 15:44:36.511] [polyfem] [info] multi_singular_boundary_count: 0 [2019-11-29 15:44:36.511] [polyfem] [info] non_regular_count: 0 [2019-11-29 15:44:36.511] [polyfem] [info] non_regular_boundary_count: 0 [2019-11-29 15:44:36.511] [polyfem] [info] undefined_count: 0 [2019-11-29 15:44:36.511] [polyfem] [info] total count: 15402 [2019-11-29 15:44:36.513] [polyfem] [info] Building isoparametric basis... [2019-11-29 15:44:36.582] [polyfem] [info] Computing polygonal basis... [2019-11-29 15:44:36.582] [polyfem] [info] took 0.000498612s [2019-11-29 15:44:36.584] [polyfem] [info] hmin: 0.915916 [2019-11-29 15:44:36.584] [polyfem] [info] hmax: 17.9633 [2019-11-29 15:44:36.584] [polyfem] [info] havg: 5.92034 [2019-11-29 15:44:36.584] [polyfem] [info] took 0.0691386s [2019-11-29 15:44:36.584] [polyfem] [info] flipped elements 0 [2019-11-29 15:44:36.584] [polyfem] [info] h: 17.9633 [2019-11-29 15:44:36.584] [polyfem] [info] n bases: 3956 [2019-11-29 15:44:36.584] [polyfem] [info] n pressure bases: 0 [2019-11-29 15:44:36.584] [polyfem] [info] Assigning rhs... [2019-11-29 15:44:36.593] [polyfem] [info] took 0.00884672s [2019-11-29 15:44:36.593] [polyfem] [info] Assembling stiffness mat... [2019-11-29 15:44:36.769] [polyfem] [info] took 0.17643s [2019-11-29 15:44:36.769] [polyfem] [info] sparsity: 427704/140849424 [2019-11-29 15:44:36.769] [polyfem] [info] Solving LinearElasticity with [2019-11-29 15:44:36.770] [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)