Boundary Condition for Mixed Elements

dolfinx
1

Hi there,

I am switching from Fenics to FenicsX, which provides some issues regarding the boundary conditions (bcs) for MixedElement. I am a little confused about how to apply, for example, a constraint to u_x at x=0. And how would it look like if I have a constraint for both components, i.e., u_x and u_y?

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
from basix.ufl import element, mixed_element
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import LinearProblem
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions,
                 div, exp, inner)

# Mesh parameters
Lx = 20.0
Ly = 2.0
Nx = 40
Ny = 4

# Define a structured quadrilateral mesh
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([Lx, Ly])], [Nx, Ny], mesh.CellType.quadrilateral)

Q_el = element("Lagrange", domain.basix_cell(), 2)
P_el = element("Lagrange", domain.basix_cell(), 1)
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(domain, V_el)

(u, Theta) = TrialFunctions(V)
(u_, Theta_) = TestFunctions(V)

fdim = domain.topology.dim - 1
facets_left = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[0], 0))

Q, _ = V.sub(0).collapse()
dofs_left = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_left)
2

You would use a similar procedure to

i.e.

P, _ = V.sub(1).collapse()
dofs_left_1 = fem.locate_dofs_topological((V.sub(1), P), fdim, facets_left)
3

Thanks for the response. I thought V.sub(0) would only access the space of u and V.sub(1) therefore would regard the variable scalar \theta here. But u=u_x, u_y is a 2D vector and therefore i should be able to somehow constraint each dimension separately, no?.

4

In your example both spaces are scalar valued.

To make Q_el a vector element you would have to call


[quote="toDumb, post:3, topic:15215, full:true"]
Thanks for the response. I thought V.sub(0) would only access the space of u and V.sub(1) therefore would regard the variable scalar \thetaθ\theta here. But u=u_x, u_yu=ux,uyu=u_x, u_y is a 2D vector and therefore i should be able to somehow constraint each dimension separately, no?.
[/quote]

In your example both spaces are scalar valued.

[quote="toDumb, post:1, topic:15215"]

Q_el = element(“Lagrange”, domain.basix_cell(), 2, shape=(domain.geometry.dim, ))

which would make a second order vector function space.

With this definition you would use 
`V.sub(0).collapse()` to get both components,
`V.sub(0).sub(0).collapse()` to get the x component only and `V.sub(0).sub(1)` to het the y component only
5

Thanks again. I corrected the mistake immediately, and it was very helpful. Now, I’ve set up the example. The resulting plot does not satisfy my left boundary condition. Where did I go wrong?

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
from basix.ufl import element, mixed_element
from dolfinx import fem, io, mesh, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions,
                 div, exp, inner)

L = 20
H = 2

domain = mesh.create_rectangle(MPI.COMM_WORLD, np.array([[0, 0], [L, H]]), [40, 4], cell_type=mesh.CellType.quadrilateral)

Q_el = element("Lagrange", domain.basix_cell(), 2, shape=(domain.geometry.dim, ))
P_el = element("Lagrange", domain.basix_cell(), 1)
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(domain, V_el)

fdim = domain.topology.dim - 1
facets_left = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[0], 0.0))

Vx, _ = V.sub(0).sub(0).collapse()
dofs_left = fem.locate_dofs_topological((V.sub(0).sub(0), Vx), fdim, facets_left)

def u_ones(x):
    values = np.ones((1, x.shape[1]))
    return values

ux_bc = fem.Function(Vx)
ux_bc.interpolate(u_ones)
bc = fem.dirichletbc(ux_bc, dofs_left, V.sub(0).sub(0))

# Test if BC are fullfilled
ux = fem.Function(Vx)
fem.petsc.set_bc(ux.vector, [bc])

import pyvista
from dolfinx import plot 
pyvista.start_xvfb()

u_linear = fem.Function(fem.functionspace(domain, element("Lagrange", domain.basix_cell(), 1)))
u_linear.interpolate(ux)
topology, cell_types, geometry = plot.vtk_mesh(domain, domain.topology.dim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
grid.point_data["u"] = u_linear.x.array.real

plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True, scalars="u")
plotter.view_xy()
plotter.save_graphic("filename.pdf")

image
image1334×334 58.3 KB