multipoint constraints

[gdmcbain]

@bhaveshshrimali, from #436:

I am interested in applying scikit-fem to solve problems with general periodic boundary conditions, which can be implemented using multipoint constraints. I have done this at a higher level in Abaqus, where the kinematics of nodes on the boundary is described using some reference nodes. However I am still looking at how to implement these within a FE code written in pure python (and also, sooner than later in scikit-fem).

I wanted to know if there are plans to support these through a higher level interface such as in dolfinx_mpc, i.e. assembling FE matrices for problems with constrained kinematics.

[gdmcbain]

These were also asked about in #426:

The library do provide a way to handle Dirichlet boundary conditions by elimination. This shall be at least mentioned. I imagine that the treatment of more complex conditions, such as multiple points constraints in mechanics, have to be treated by hand after the assembly, isn't it ?

There was some discussion of ‘how DOFs are handled’ in #429; I sketched some ideas there but was without a clear example of the problem so it's a bit abstract.

What's a ‘general periodic boundary condition’? Is that the value at some point on the boundary is equal to that on another plus a jump? (Like maybe in a viscous incompressible fluid problem with the pressure in a length of duct, with velocity being periodic from outlet to inlet but the pressure being more or less periodic but with an overall longitudinal gradient?) Or is it more general than that, with some functional relation between the values (and normal gradients?) at the two points? Or are more than two points coupled?

I'll look at dolfinx_mpc. Maybe something like https://github.com/jorgensd/dolfinx_mpc/blob/master/python/demos/demo_periodic.py ?

[bhaveshshrimali]

What's a ‘general periodic boundary condition’? Is that the value at some point on the boundary is equal to that on another plus a jump? (Like maybe in a viscous incompressible fluid problem with the pressure in a length of duct, with velocity being periodic from outlet to inlet but the pressure being more or less periodic but with an overall longitudinal gradient?) Or is it more general than that, with some functional relation between the values (and normal gradients?) at the two points? Or are more than two points coupled?

Thanks for such a prompt reply. Yes that's exactly what I was referring to by periodic. By general periodic, I meant something like here.
I used the term probably a bit more loosely than should be. By general I wanted to denote a set of linear constraints that tie the dofs on opposite facets of a domain. For instance in elasticity, on a simple unit square, it would amount to saying

# u_R: displacement on the right edge (x=1)
# u_L: displacement on the left edge (x=0)
# u_T: displacement on the right edge (y=1)
# u_B: displacement on the left edge (y=0)

# u_BR: displacement of the bottom-right corner (1,0)
# u_TL: displacement of the top-left corner (0,1)

# u_o: displacement of the origin (0,0)

u_R - u_L = u_BR - u_o
u_T - u_B = u_TL - u_o

More often than not, I specify Dirichlet BC on u_o (can very well be another point in the domain), and treat u_BR/u_TL as control points, either (or both) of them could have prescribed (Dirichlet) values.

Implementation-wise this would be slightly more general than simply periodic boundary conditions (u_L = u_R and u_T = u_B), as it would even allow us to impose a combination of dirichlet and neumann bcs. If, say, the material is incompressible and I apply u_BR = lmbda (Dirichlet) and nothing else (which is equivalent to Neumann BC on the top edge), then u_TL = 1./(1+lmbda)**(0.5) and accordingly the displacements at the top and right edges.

Like maybe in a viscous incompressible fluid problem with the pressure in a length of duct, with velocity being periodic from outlet to inlet but the pressure being more or less periodic but with an overall longitudinal gradient?

Yes. Infact if I am thinking properly, this would be analogous to the incompressible elasticity problem (with mixed Dirichlet/Neumann BCs) that I sketched above.

Thanks again for your reply. Sorry if this sounds too niche. However once I pick up some practice, I can maybe assist in some implementation too.

I will start this week by implementing a simple compressible hyperelasticity (like FEniCS demo), then hopefully incompressible hyperelasticity and eventually revisit this issue.

[gdmcbain]

No, I don't think this is too niche. I haven't needed it myself yet, but doubtless will sooner or later.

For identifying the periodic boundaries, my first thought would be to make use (in two dimensions) of Periodic Curve in Gmsh. I think that part of a Gmsh MSH is read by meshio.read, but not yet by skfem.io.meshio.from_meshio.

[bhaveshshrimali]

Thanks! I was starting to take a look at implementing the hyperelasticity demo. The most natural examples to follow seem to be linear elasticity, nonlinear poisson and laplace with inhomogeneous bcs.

What would be the best place to look for functions like log, inv, det (like in UFL) when writing the weak forms? For instance the stress would have an expression that looks like the following in UFL:

def firstPKStress(u):
    F = Identity(len(u)) + grad(u) 
    J = det(F)
    return mu * F - mu * inv(F).T + J * (J-1) * inv(F).T

I can see the helper functions eye and grad, which should help in defining F as eye(1, u.shape[0]) + grad(u), but how about the determinant (det) and inverse (inv) ?

For identifying the periodic boundaries, my first thought would be to make use (in two dimensions) of Periodic Curve in Gmsh. I think that part of a Gmsh MSH is read by meshio.read, but not yet by skfem.io.meshio.from_meshio.

I see; will take a look. This should save any post-processing like looking up matching vertices. Thanks!

For a simple 2D domain, generating a periodic mesh itself shouldn't be a problem as Christophe mentioned in a discussion sometime back. I have tested this with FEniCS/Abaqus. For generating the linear constraints however, I always did a simple lookup on the generated mesh (using numpy.where) to get the coordinates of matching points on opposite facets. For 3D meshes I have used Netgen which gives the corresponding node-pairs and moreover the list of linear constraints

[gdmcbain]

I'll attempt to answer the question on writing forms in #439.

[kinnala]

I once enforced a periodic solution like this:

    # periodic boundary condition
    Nleft=m.nodes_satisfying(lambda x,y:x==0.0)
    meps = 1e-12
    Nright=m.nodes_satisfying(lambda x,y:(x<2.0*np.pi+meps)*(x>2.0*np.pi-meps))
    Nbottom=m.nodes_satisfying(lambda x,y:y==0.0)
    Nleft=np.setdiff1d(Nleft,N1)
    Nright=np.setdiff1d(Nright,N1)
    Nelse=np.setdiff1d(I,Nleft)
    Nelse=np.setdiff1d(Nelse,Nright)

    I1=np.arange(Nelse.shape[0])
    I2=np.arange(Nleft.shape[0])+I1[-1]+1

    def periodify_matrix(Z):
        Zxx=Z[Nelse].T[Nelse].T
        Zx1=Z[Nelse].T[Nleft].T
        Z1x=Zx1.T
        Zx2=Z[Nelse].T[Nright].T
        Z2x=Zx2.T
        Z11=Z[Nleft].T[Nleft].T
        Z12=Z[Nleft].T[Nright].T
        Z21=Z12.T
        Z22=Z[Nright].T[Nright].T

        Zp1=spsp.hstack((Zxx,Zx1+Zx2))
        Zp2=spsp.hstack((Z1x+Z2x,Z11+Z12+Z21+Z22))
        return spsp.vstack((Zp1,Zp2)).tocsr()

    def periodify_vector(y):
        return np.hstack((y[Nelse],y[Nleft]+y[Nright]))

    def unpack_periodic(y):
        Y=np.zeros(f.shape[0])
        Y[Nelse]=y[I1]
        Y[Nleft]=y[I2]
        Y[Nright]=y[I2]

    Ap=periodify_matrix(A)
    x0=f*0.0
    x0[Nbottom]=other_BC
    fp=periodify_vector(f-A.dot(x0))

    up[0]=fsol.direct(Ap,fp,use_umfpack=True)
    # unpack periodic
    U=unpack_periodic(up[0])

It's from the numerical example of this paper: https://arxiv.org/abs/1711.04274. I think it works only for linear elements though. Idea is to eliminate the matching DOF's.

There are several other options:

• Lagrange multipliers
• weak enforcement (penalty or Nitsche's method)
• modify mesh.t so that matching vertices have same index (haven't tested this so there might be caveats)