Mixed boundary conditions

[WenjieYao]

I want to solve a problem with periodic boundary conditions on the left and right side, the top and bottom side have nonzero Dirichlet boundary conditions in some segments and Neumann boundary conditions in other segements, how should I do that? I tried some codes from expample 13 and 42, but I got errors Transforming over 1000 vertices to C_CONTIGUOUS. when I execute m = MeshTri1DG.init_tensor( np.linspace(0, 1, 30), np.linspace(0, 1, 30), periodic=[0], ).with_boundaries({ 'ground': lambda x: x[1] == 0., 'positive': lambda x: x[1] == 1, })

[kinnala]

Transforming over 1000 vertices to C_CONTIGUOUS is not an error but a warning. It is of no concern if you do not know what it means. I'll transform this to Discussions.

[kinnala]

I suggest you don't use MeshTri1DG because there might be unexpected results if your segments are touching the periodic boundary.

Try to instead use the following snippet as a starting point:

import numpy as np
from skfem import *
from skfem.models import laplace, unit_load

m = MeshTri().refined(4)
e = ElementTriP1()

basis = Basis(m, e)

left = basis.get_dofs(facets='left').nodal['u']
right = basis.get_dofs(facets='right').nodal['u']
D = basis.get_dofs({'top', 'bottom'}).all()

left = np.setdiff1d(left, D)
right = np.setdiff1d(right, D)

B = np.zeros((len(left), basis.N))
B[range(len(left)), left] = 1
B[range(len(right)), right] = -1

A = asm(laplace, basis)
f = asm(unit_load, basis)

# create a saddle point system with one additional Lagrange multiplier / periodic DOF pair
A = bmat([[A, B.T],
          [B, None]], 'csr')
f = np.concatenate((f, 0 * left))

y = solve(*condense(A, f, D=basis.get_dofs({'top', 'bottom'})))

basis.plot(y[:basis.N], colorbar=True, shading='gouraud', levels=5).show()

It uses Lagrange multipliers to enforce periodicity.

[kinnala]

Here is the file: periodic.msh.txt

[kinnala]

This is how to load it using meshio:

In [15]: import meshio

In [16]: meshio.read("periodic.msh")

Out[16]: 
<meshio mesh object>
  Number of points: 142
  Number of cells:
    line: 10
    line: 10
    line: 10
    line: 10
    triangle: 242
  Cell sets: bottom, right, top, left, all, gmsh:bounding_entities
  Point data: gmsh:dim_tags
  Cell data: gmsh:physical, gmsh:geometrical
  Field data: bottom, right, top, left, all

I cannot find here anything which refers to it being a periodic mesh. Perhaps one of those datasets with gmsh: prefix has the information.

[kinnala]

Anyways, since you have those left and right tags, it should be easy to use those to solve a periodic problem. Just copypaste the above code I gave you:

import numpy as np
from skfem import *
from skfem.models import laplace, unit_load

m = MeshTri.load("periodic.msh")
e = ElementTriP1()

basis = Basis(m, e)

left = basis.get_dofs(facets='left').nodal['u']
right = basis.get_dofs(facets='right').nodal['u']
D = basis.get_dofs({'top', 'bottom'}).all()

left = np.setdiff1d(left, D)
right = np.setdiff1d(right, D)

B = np.zeros((len(left), basis.N))
B[range(len(left)), left] = 1
B[range(len(right)), right] = -1

A = asm(laplace, basis)
f = asm(unit_load, basis)

# create a saddle point system with one additional Lagrange multiplier / periodic DOF pair
A = bmat([[A, B.T],
          [B, None]], 'csr')
f = np.concatenate((f, 0 * left))

y = solve(*condense(A, f, D=basis.get_dofs({'top', 'bottom'})))

basis.plot(y[:basis.N], colorbar=True, shading='gouraud', levels=5).show()

[WenjieYao]

I think mehsio doesn't read the periodic information from a gmsh created file, you need to import gmsh to get the information...

[kinnala]

Ok, then this is a feature request for meshio since we won't import gmsh in scikit-fem due to its license.