Boundary constraint on space

[Abelsylowlie]

If define function space as follows
,
how to build the block matrices by using skfem.element.ElementVector?
@kinnala @gdmcbain

[Abelsylowlie]

Just like ex03, but here boundary condition is neither Dirichlet nor Neumann, it's constraint on two spaces.

[kinnala]

I suggest you do one H^1 at a time and assemble two matrices. Then create a block matrix with scipy.sparse.bmat and augment the linear system with Lagrange multipliers that enforce the jump to zero.

[kinnala]

You can also penalize though. I'll try it later.

[kinnala]

What happens to normal gradient?

[Abelsylowlie]

And there is no relation between first two equation, third and fourth build the relation.

[gdmcbain]

I thunk this is an example of a multipoint constraint . #718 (That's a link to the more general discussion, not necessarily the solution.)

What i think i would do here is work out the dimension of the constrained degrees of freedom and then say that if I knew those values, I could compute the rest of the phi and psi fields, as a function of those values treated as Dirichlet data, then I could compute whatever functionals of phi and psi that completed their definition. Finally treating those functionals as functions of the assumed Dirichlet data, I solve to find the Dirichlet data that give the desired solution. Either directly (Lagrange multiplier) or iteratively.

[Abelsylowlie]

Thanks for your rely. I'll try it after I understand this procedure.

[gdmcbain]

Yes it's basically the same as that expressed more clearly by kinnala, which I hadn't read when I wrote.

[kinnala]

Here is example using the penalty method: https://colab.research.google.com/drive/1EzmM-IzMxJLrobE1PE9r8ddgBTU6Y6XG?usp=sharing

I have boundary conditions phi(0, y) = 0 and psi(0, y) = 1 and the jump condition phi - psi penalized for top-right part of the boundary, (x > 0.5) * (y > 0.5), by adding a penalty term to the bilinear form. You could use a Nitsche-type approach instead of penalty if you don't want to fine tune the penalty parameter (Nitsche converges for a larger range of parameters) or use a higher order method.

Requires master version because it uses those unreleased left, right, etc. tags.

[kinnala]

phi:

psi:

[Abelsylowlie]

Thank you very much for these insights.

[kinnala]

I'm adding the source code also here for reference:

from skfem import *
from skfem.visuals.matplotlib import *
from skfem.models import *

m = MeshTri().refined(5)
e = ElementTriP1() * ElementTriP1()
basis = Basis(m, e)
fbasis = FacetBasis(m, e, facets=lambda x: (x[0] > 0.5) * (x[1] > 0.5) * ((x[0] == 1) + (x[1] == 1)))

from skfem.helpers import *

@BilinearForm
def penalty(phi, psi, phit, psit, w):
  return 100. * (phi - psi) * (phit - psit)

@BilinearForm
def double_laplace(phi, psi, phit, psit, w):
  return dot(grad(phi), grad(phit)) + dot(grad(psi), grad(psit))

A = double_laplace.assemble(basis) + penalty.assemble(fbasis)

u = basis.zeros()
u[basis.get_dofs('left').all('u^1')] = 0
u[basis.get_dofs('left').all('u^2')] = 1

x = solve(*condense(A, 0.*u, x=u, D=basis.get_dofs('left')))

(phi, phibasis), (psi, psibasis) = basis.split(x)

import matplotlib.pyplot as plt

plot(phibasis, phi, colorbar=True, vmin=0, vmax=1, cmap=plt.cm.jet)
plot(psibasis, psi, colorbar=True, vmin=0, vmax=1, cmap=plt.cm.jet)

[gdmcbain]

This is very nice. It could be a good addition to docs/examples.