Examples/snippets demonstrating the integrations

[kinnala]

I'd like to keep a track of the libraries that can be integrated to scikit-fem with ease. How about having a single example which simply demonstrates the possible integrations? Or should it be a collection of snippets instead? If the latter, how is it tested in pytest?

[gdmcbain]

I think a collection; the integrations serve diverse ends.

For pytest, I've seen this idiom:

netCDF4 = pytest.importorskip("netCDF4")

https://github.com/nschloe/meshio/blob/3a380fc1100827ed38b8f38960087a066c4aab26/test/test_exodus.py#L6

[kinnala]

So basically in practice we could have have something like docs/integrations/pygmsh.py, docs/integrations/dmsh.py, and so on, make sure that these are collected by pytest, and then maybe have docs/snippets.rst which includes all these in the Sphinx documentation?

[gdmcbain]

Yes, I think so. That directory structure should make it easy to find the case of interest.

Is there a catch with trying to import pygmsh from inside a file.called pygmsh.py or from another file in a directory with that file?

[kinnala]

Yeah, let's call them mesh_pygmsh.py, mesh_dmsh.py, solver_pacopy.py, and so on, i.e. try to have a 'category' which describes the type of integration?

[kinnala]

Here are some integrations that come to mind:

• pygmsh
• pacopy
• dmsh
• pyamgcl
• sksparse
• optimesh
• pyfvm (moving mesh to and from, e.g., for iterative coupling)
• vedo (actually I should continue on vedo, formerly known as vtkplotter #361 instead)

I'm sure there are others. Feel free to add to the list.

Edit: I'm sure there are also those that might be hard to test in our test suite such as reading a mesh from FEniCS or ngsolve. Not sure if we should include those or only what we are willing to test?

[kinnala]

Some additional urls for consideration:

• https://github.com/andrenarchy/krypy
• https://pymesh.readthedocs.io/en/latest/api_sparse_solver.html
• https://github.com/nschloe/pygalmesh
• https://github.com/nschloe/meshzoo
• https://github.com/pyvista/pyvista
• https://pypi.org/project/petsc4py/

[gdmcbain]

There's a start on PETSc in #236 .

After that recent question about pressure boundary conditions for the Stokes equation #466 , which referred to a FEniCS example using MUMPS, I started looking at PyMumps.

nschloe has a new package picking up from krypy: Krylov.

I've been looking at pyMOR; besides it being useful for scikit-fem applications, I think scikit-fem.could be useful for it. It has bindings to FEniCS, deal.ii, and NGSolve, and its own little pythonic FEM implementation but that's rather limited.

I like the idea of coupling to a finite volume solver. For finite differences, there's scikit-fdiff. A boundary element code could be really useful too.

[gdmcbain]

There's been some progress on the pyMOR side in pymor/pymor#1507 ‘Experimental scikit-fem support’; I'm very keen to see how this goes and try it out myself too.

[gdmcbain]

I had been that thinking optimesh might be useful with the new MeshTri.init_circle #475 in examples #476 but actually it doesn't seem to be required #477 .

[kinnala]

I think that in near future I might be reimplementing one specific algorithm from optimesh that I'm very fond of and also using in https://github.com/kinnala/adaptmesh.

[gdmcbain]

For external reasons but possibly of interest here, I've posted an example of integration with OpenFOAM.

• https://github.com/gdmcbain/OpenFOAM-9/tree/d4afc0eb4e000d82ecd81dd72dd5da3c9c808e03/tutorials/basic/potentialFoam/cylinderSkFem

It generates a quasi-two-dimensional hexahedral mesh using OpenFOAM's blockMesh as in the OpenFOAM tutorial basic/potentialFoam/cylinder, one hexahedron thick in the z-direction, which is how two dimensions is done in OpenFOAM. This is then parsed in a small general Python submodule to create a skfem.MeshQuad for the domain with .boundaries for the inlet. The inlet boundary is also saved as a trace mesh (#530 #687) and the discrete finite volume inlet velocity boundary condition computed in scikit-fem and written for potentialFoam. The same problem is then solved in both codes.

The creation of the skfem.Mesh for domain and boundary and their writing and rereading was greatly simplified following the recent improvements to skfem.io.meshio #680.

[adtzlr]

Hi, as a quick follow-up from #839 here is a minimal-modified version of example 43 with a hyperelastic Neo-Hookean material formulation defined by matADi (pip install matadi==0.0.28) instead of the hyperelastic form coded by hand.

scikit-fem with matadi

import numpy as np

from skfem import *
from skfem.helpers import grad, identity

from matadi import MaterialHyperelastic
from matadi.models import neo_hooke

umat = MaterialHyperelastic(fun=neo_hooke, C10=0.5, bulk=2.0)
mesh = MeshHex().refined(2).with_boundaries({
    'left': lambda x: x[0] == 0.,
    'right': lambda x: x[0] == 1.,
})
element = ElementVector(ElementHex1())
basis = Basis(mesh, element, intorder=1)
u = basis.zeros()

dofs = basis.get_dofs('right')
dirichlet = basis.get_dofs({'right', 'left'})
active = basis.complement_dofs(dirichlet)

Eye = identity(grad(basis.interpolate(u)))
F = grad(basis.interpolate(u)) + Eye

@LinearForm
def L(v, w):
    return umat.gradient_vector_product([F], [grad(v)])[0]

@BilinearForm
def a(u, v, w):
    return umat.hessian_vector_product([F], [grad(v)], [grad(u)])[0]

b = None
right = -0.5

for iteration in range(8):

    K = a.assemble(basis)

    du_D = u.copy()
    du_D[dofs.nodal['u^1']] -= right

    u += solve(*condense(K, b, x=-du_D, D=dirichlet))

    F = grad(basis.interpolate(u)) + Eye
    b = -L.assemble(basis)
    norm_f = np.linalg.norm(-b[active])

    print(1 + iteration, norm_f)

    if norm_f < 1e-10:
        break

output:

1 0.2313927273813761
2 0.018997771204881598
3 4.554610746120609e-05
4 7.916934983340329e-09
5 7.069702684553263e-16