Skip to content

Commit

Permalink
modernize ex09 and fix vedo to work with the latest release
Browse files Browse the repository at this point in the history
  • Loading branch information
kinnala committed Jun 17, 2024
1 parent bd16209 commit e5b0fac
Show file tree
Hide file tree
Showing 2 changed files with 28 additions and 33 deletions.
54 changes: 24 additions & 30 deletions docs/examples/ex09.py
Original file line number Diff line number Diff line change
Expand Up @@ -35,35 +35,41 @@
* Demidov, D. (2019). AMGCL: an efficient, flexible, and extensible algebraic multigrid implementation. `arXiv:1811.05704 <https://arxiv.org/abs/1811.05704>`_
"""

from skfem import *
from skfem.models.poisson import *
from skfem.helpers import *
import numpy as np
from scipy.sparse import spmatrix
from scipy.sparse.linalg import LinearOperator


p = np.linspace(0, 1, 16)
m = MeshTet.init_tensor(*(p,) * 3)
basis = Basis(m, ElementTetP1())

e = ElementTetP1()
basis = Basis(m, e)

A = asm(laplace, basis)
b = asm(unit_load, basis)
@BilinearForm
def laplace(u, v, _):
return dot(grad(u), grad(v))

I = m.interior_nodes()

@LinearForm
def unit_load(v, _):
return 1. * v


A = laplace.assemble(basis)
b = unit_load.assemble(basis)

I = m.interior_nodes()
x = 0. * b

if __name__ == "__main__":
verbose = True
else:
verbose = False

Aint, bint = condense(A, b, I=I, expand=False)

preconditioners = [None, build_pc_ilu(Aint, drop_tol=1e-3)]


# try importing algebraic multigrid from external package
try:
from pyamg import smoothed_aggregation_solver

Expand All @@ -76,29 +82,17 @@ def build_pc_amgsa(A: spmatrix, **kwargs) -> LinearOperator:
except ImportError:
print('Skipping PyAMG')

try:
import pyamgcl

def build_pc_amgcl(A: spmatrix, **kwargs) -> LinearOperator:
"""AMG preconditioner"""

if hasattr(pyamgcl, 'amgcl'): # v. 1.3.99+
return pyamgcl.amgcl(A, **kwargs)
else:
return pyamgcl.amg(A, **kwargs)

preconditioners.append(build_pc_amgcl(Aint))

except ImportError:
print('Skipping pyamgcl')

# solve for each preconditioner
for pc in preconditioners:
x[I] = solve(Aint, bint, solver=solver_iter_pcg(verbose=verbose, M=pc))
x[I] = solve(Aint, bint, solver=solver_iter_pcg(verbose=True, M=pc))


if verbose:
if __name__ == "__main__":
from os.path import splitext
from sys import argv

m.draw('vedo', point_data={'potential': x}).show()
#m.save(splitext(argv[0])[0] + ".vtk", {'potential': x})
# use vedo: press 5 to visualize potential, X for cutter tool
basis.plot(x, 'vedo').show()
# preferred: save to vtk for visualization in Paraview
#m.save(splitext(argv[0])[0] + ".vtk", point_data{'potential': x})
7 changes: 4 additions & 3 deletions skfem/visuals/vedo.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,17 +3,18 @@

def draw(m, backend=False, **kwargs):
"""Visualize meshes."""
from vedo import Plotter, UGrid
from vedo import Plotter, UnstructuredGrid, show
vp = Plotter()
grid = None
with tempfile.NamedTemporaryFile() as tmp:
m.save(tmp.name + '.vtk',
encode_cell_data=False,
encode_point_data=True,
**kwargs)
grid = UGrid(tmp.name + '.vtk')
grid = UnstructuredGrid(tmp.name + '.vtk')
vp += grid.tomesh()
# save these for further use
grid.show = lambda: vp.show([grid.tomesh()]).close()
grid.show = lambda: vp.show()
grid.plotter = vp
return grid

Expand Down

0 comments on commit e5b0fac

Please sign in to comment.