Struggling with vector values bcdirichlet

I am trying to formulate a dirichlet boundary condition for a problem with two unknown variables, and the boundary condition should only be applied to the second unknown variable, but not the first. In the sample code below, it appears that my code functions if I use a boundary condition with
value=default_scalar_type(0)
However, if I set
value= u_D
where u_D is an interpolated function, the code crashes.
How must I formulate u_D to get the boundary condition to work?

Here is the code (I’m sorry, I don’t know how to paste code into a separate window, as everyone else here seems to do…). If you change iflag from 0 to 1 at the bottom of the code, you will have the option of using value=default_scalar_type(0) or value= u_D.
You can see that the code runs with iflag=0 but crashes with iflag=1.


‘’’ python
from dolfinx.fem import (Constant, dirichletbc, Function, functionspace, locate_dofs_geometrical,
locate_dofs_topological)
from mpi4py import MPI
from dolfinx import fem, mesh, default_scalar_type
import numpy as np
import ufl

Define mesh

nx, ny = 50, 50
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([-2, -2]), np.array([2, 2])],
[nx, ny], mesh.CellType.triangle)

VP=fem.functionspace(domain, (“Lagrange”, 1, (2 , )))

Create facet to cell connectivity required to determine boundary facets

tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs_x = locate_dofs_topological(VP.sub(1), domain.topology.dim - 1, boundary_facets)

V=fem.functionspace(domain, (“Lagrange”, 1))
x = ufl.SpatialCoordinate(domain)
u_D=Function(V)
u_D.interpolate(lambda x:x[0])
iflag=1
if iflag==0:
value=default_scalar_type(0)
else:
value=u_D

bcx = dirichletbc(value, boundary_dofs_x, VP.sub(1))
‘’’

Please format the code appropriately, i.e using

```python
# Add code here

```

I’m not sure I understand what you requested. I went back and used ‘’’ python code ‘’’ Is that what you were requesting?

No, please use the backticks i had in the previous post `, as it is standard markdown for rendering code, making indentation and various other things format correctly .

Here is a second try to reproduce the code

from dolfinx.fem import (Constant, dirichletbc, Function, functionspace, locate_dofs_geometrical,
                         locate_dofs_topological)
from mpi4py import MPI
from dolfinx import fem, mesh, default_scalar_type
import numpy as np
import ufl


# Define mesh
nx, ny = 50, 50
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([-2, -2]), np.array([2, 2])],
                               [nx, ny], mesh.CellType.triangle)

VP=fem.functionspace(domain, ("Lagrange", 1, (2 , )))

# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs_x = locate_dofs_topological(VP.sub(1), domain.topology.dim - 1, boundary_facets)
V=fem.functionspace(domain, ("Lagrange", 1))
x = ufl.SpatialCoordinate(domain)
u_D=Function(V)
u_D.interpolate(lambda x:x[0])
iflag=1
if iflag==0:
    value=default_scalar_type(0)
else:
    value=u_D

bcx = dirichletbc(value, boundary_dofs_x, VP.sub(1))

Thanks, I learned something new today. My keyboard doesn’t have a backtick mark, but I can type one with Alt 96

The first usage is also not correct.
Here is a code explaining how each bc function works:

from mpi4py import MPI
from dolfinx import fem, mesh, default_scalar_type, io as _io
import numpy as np


# Define mesh
nx, ny = 50, 50
domain = mesh.create_rectangle(
    MPI.COMM_WORLD,
    [np.array([-2, -2]), np.array([2, 2])],
    [nx, ny],
    mesh.CellType.triangle,
)

VP = fem.functionspace(domain, ("Lagrange", 1, (2,)))

# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)

iflag = 1
if iflag == 0:
    # Create boundary condition based on a sub-space

    # Get subspace and create a function in those space to put BC in
    V_1, _ = VP.sub(1).collapse()
    u_D = fem.Function(V_1)
    u_D.interpolate(lambda x: x[0])

    # Locate degrees of freedom in the collapsed sub space, and the VP sub space
    boundary_dofs_x = fem.locate_dofs_topological(
        (VP.sub(1), V_1), fdim, boundary_facets
    )
    # We observe that boundary_dof_x consist of two arrays
    print(boundary_dofs_x)


else:
    # Create a boundary condition based on a constant value
    # We do not need to collapse to a sub space, as a constant is not in a subspace
    boundary_dofs_x = fem.locate_dofs_topological(VP.sub(1), fdim, boundary_facets)
    u_D = default_scalar_type(10)

bcx = fem.dirichletbc(u_D, boundary_dofs_x, VP.sub(1))


# Check validity of bc
u_tmp = fem.Function(VP)
fem.set_bc(u_tmp.x.array, [bcx])

with _io.VTXWriter(domain.comm, f"u_{iflag}.bp", [u_tmp], engine="BP4") as bp:
    bp.write(0.0)

@dokken, you have the patience of a saint! Thank you so much. There is a lot about the grammar and rules you have used that I don’t understand, but can imitate. Is there some better way for me to learn about these structures?

Now that the boundary conditions don’t crash, I am finding that the NonlinearProblem is still not converging when I use them. I will spend some time seeing if I can fix this problem myself first. If I fail, I will post again to ask more questions.

Thanks again for your help

I have made various tutorials at:
https://jsdokken.com/dolfinx-tutorial/

I am a little confused about the following

[quote=“dokken, post:6, topic:14881”]

   # Create boundary condition based on a sub-space

    # Get subspace and create a function in those space to put BC in
    V_1, _ = VP.sub(1).collapse()
    u_D = fem.Function(V_1)
    u_D.interpolate(lambda x: x[0])

    # Locate degrees of freedom in the collapsed sub space, and the VP sub space
    boundary_dofs_x = fem.locate_dofs_topological(
        (VP.sub(1), V_1), fdim, boundary_facets
    )

In your tutorial
https://jsdokken.com/dolfinx-tutorial/chapter3/component_bc.html
you have what seems to be a similar boundary condition, but the formulation is

[Quote]
Next, we locate the degrees of freedom on the top boundary. However, as the boundary condition is in a sub space of our solution, we need to supply both the parent space 𝑉 and the sub space 𝑉0 to dolfinx.locate_dofs_topological.

def right(x): return np.logical_and(np.isclose(x[0], L), x[1] < H)

boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, right) boundary_dofs_x = locate_dofs_topological(V.sub(0), mesh.topology.dim - 1, boundary_facets)

We can now create our Dirichlet condition

bcx = dirichletbc(default_scalar_type(0), boundary_dofs_x, V.sub(0))

[Quote]
Why is first argument of fem.locate_dofs_topological simply V.sub(0) in your example whereas above the corresponding argument is (VP.sub(1), V_1)?

As mentioned above, one is for using dolfinx.fem.Function from a collapsed sub-space as a boundary condition. In that case you supply the mapping between the dofs of the collapsed space and the parent sub space, as well as what sub space one is applying it to.

The one from the tutorial applies a constant value, which does not belong to a function space, thus you need no map, but you need to supply the sub space.

Suppose that in my variational formulation I want to give a function f, where
del^2 phi = f
If this function f is also an analytic function, would I also use

V_1, _ = VP.sub(1).collapse()
 f = fem.Function(V_1)
 f.interpolate(lambda x: x[0])

?

Since f is going to be part of your variational form, you can either define it as a function, or by using ufl.SpatialCoordinate, as shown in
https://jsdokken.com/dolfinx-tutorial/chapter1/membrane_code.html#defining-a-spatially-varying-load

Thanks for your help!

@dokken I have two question
1.What is the error in imposing the input condition in this way?
2.How can I draw the boundary conditions to verify them?
This is the code

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista

from dolfinx.fem import Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx.io import VTXWriter
from dolfinx import mesh
from dolfinx.plot import vtk_mesh
from basix.ufl import element
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
                 div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym)

t = 0
T = 10
num_steps = 500
dt = T / num_steps

lenght = 10
width = 10

#create the rectangular mesh
mesh = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([lenght, width])],
                               [20, 20], cell_type=mesh.CellType.triangle)
t = 0
T = 10
num_steps = 500
dt = T / num_steps

v_cg2 = element("Lagrange", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim, ))
s_cg1 = element("Lagrange", mesh.topology.cell_name(), 1)
V = functionspace(mesh, v_cg2)
Q = functionspace(mesh, s_cg1)

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

fdim = mesh.topology.dim - 1


def walls(x):
    return np.isclose(x[1], 0)


wall_dofs = locate_dofs_geometrical(V, walls)
u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bcu_noslip = dirichletbc(u_noslip, wall_dofs, V)


def inflow(x):
    return np.isclose(x[0], 0)


u_ref = 5.5
z_ref = 2.0


class InletVelocity():
    def __call__(self, x):
        values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
        values[0] = u_ref * (x[1] / z_ref) ** 0.16
        return values


inflow_dofs = locate_dofs_geometrical(V, inflow)
# Inlet
u_inlet = Function(V)
inlet_velocity = InletVelocity()
u_inlet.interpolate(inlet_velocity)
bcu_inflow = dirichletbc(u_inlet, inflow_dofs, V)

Please note that your code has way too many imports and unused code. Here is a cleaned up version of your code:

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np

from dolfinx.fem import (
    Function,
    functionspace,
    dirichletbc,
    locate_dofs_geometrical,
)
from dolfinx import mesh
from basix.ufl import element
from ufl import (
    TestFunction,
    TrialFunction,
)

t = 0
T = 10
num_steps = 500
dt = T / num_steps

lenght = 10
width = 10

# create the rectangular mesh
mesh = mesh.create_rectangle(
    MPI.COMM_WORLD,
    [np.array([0, 0]), np.array([lenght, width])],
    [20, 20],
    cell_type=mesh.CellType.triangle,
)
t = 0
T = 10
num_steps = 500
dt = T / num_steps

v_cg2 = element("Lagrange", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim,))
s_cg1 = element("Lagrange", mesh.topology.cell_name(), 1)
V = functionspace(mesh, v_cg2)
Q = functionspace(mesh, s_cg1)

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

fdim = mesh.topology.dim - 1


def walls(x):
    return np.isclose(x[1], 0)


wall_dofs = locate_dofs_geometrical(V, walls)
u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bcu_noslip = dirichletbc(u_noslip, wall_dofs, V)


def inflow(x):
    return np.isclose(x[0], 0)


u_ref = 5.5
z_ref = 2.0


class InletVelocity:
    def __call__(self, x):
        values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
        values[0] = u_ref * (x[1] / z_ref) ** 0.16
        return values


inflow_dofs = locate_dofs_geometrical(V, inflow)
# Inlet
u_inlet = Function(V)
inlet_velocity = InletVelocity()
u_inlet.interpolate(inlet_velocity)
bcu_inflow = dirichletbc(u_inlet, inflow_dofs,V)

You get an error (which you should post next time)

Traceback (most recent call last):
  File "/root/shared/mwe.py", line 80, in <module>
    bcu_inflow = dirichletbc(u_inlet, inflow_dofs, V)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/bcs.py", line 189, in dirichletbc
    bc = bctype(_value, dofs, V._cpp_object)
TypeError: __init__(): incompatible function arguments. The following argument types are supported:
    1. __init__(self, g: ndarray[dtype=float64, writable=False, order='C'], dofs: ndarray[dtype=int32, writable=False, shape=(*), order='C'], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None
    2. __init__(self, g: dolfinx.cpp.fem.Constant_float64, dofs: ndarray[dtype=int32, writable=False, shape=(*), order='C'], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None
    3. __init__(self, g: dolfinx.cpp.fem.Function_float64, dofs: ndarray[dtype=int32, writable=False, shape=(*), order='C']) -> None
    4. __init__(self, g: dolfinx.cpp.fem.Function_float64, dofs: list[ndarray[dtype=int32, writable=False, shape=(*), order='C']], V: dolfinx.cpp.fem.FunctionSpace_float64) -> None

Invoked with types: dolfinx.cpp.fem.DirichletBC_float64, dolfinx.cpp.fem.Function_float64, ndarray, dolfinx.cpp.fem.FunctionSpace_float64

because you are sending in V as the last argument to dirichletbc.
The reason for this throwing an error is because you are sending in a function from the space V, and then you only need two inputs to dirichletbc:

bcu_inflow = dirichletbc(u_inlet, inflow_dofs)

Consider the following MWE

from dolfinx.fem import set_bc
from dolfinx.io import VTXWriter

u_bc = Function(V)
u_bc.x.array[:] = 3.14  # Some unique value
set_bc(u_bc.x.array, [bcu_inflow, bcu_noslip])
with VTXWriter(mesh.comm, "bc_test.bp", [u_bc], engine="BP4") as bp:
    bp.write(0.0)

which yields

1 Like

I have tried to use the code, but I have problems with paraview. This is my code. It is the same as the previous one.

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np

from dolfinx.fem import (
    Function,
    functionspace,
    dirichletbc,
    locate_dofs_geometrical,
)
from dolfinx import mesh
from basix.ufl import element
from ufl import (
    TestFunction,
    TrialFunction,
)

t = 0
T = 10
num_steps = 500
dt = T / num_steps

lenght = 10
width = 10

# create the rectangular mesh
mesh = mesh.create_rectangle(
    MPI.COMM_WORLD,
    [np.array([0, 0]), np.array([lenght, width])],
    [20, 20],
    cell_type=mesh.CellType.triangle,
)
t = 0
T = 10
num_steps = 500
dt = T / num_steps

v_cg2 = element("Lagrange", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim,))
s_cg1 = element("Lagrange", mesh.topology.cell_name(), 1)
V = functionspace(mesh, v_cg2)
Q = functionspace(mesh, s_cg1)

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

fdim = mesh.topology.dim - 1


def walls(x):
    return np.isclose(x[1], 0)


wall_dofs = locate_dofs_geometrical(V, walls)
u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bcu_noslip = dirichletbc(u_noslip, wall_dofs, V)


def inflow(x):
    return np.isclose(x[0], 0)


u_ref = 5.5
z_ref = 2.0


class InletVelocity:
    def __call__(self, x):
        values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
        values[0] = u_ref * (x[1] / z_ref) ** 0.16
        return values


inflow_dofs = locate_dofs_geometrical(V, inflow)
# Inlet
u_inlet = Function(V)
inlet_velocity = InletVelocity()
u_inlet.interpolate(inlet_velocity)
bcu_inflow = dirichletbc(u_inlet, inflow_dofs)

from dolfinx.fem import set_bc
from dolfinx.io import VTXWriter

u_bc = Function(V)
u_bc.x.array[:] = 3.14  # Some unique value
set_bc(u_bc.x.array, [bcu_inflow, bcu_noslip])
with VTXWriter(mesh.comm, "bcu_test.bp", [u_bc], engine="BP4") as bp:
    bp.write(0.0)

You need to open the whole .bp folder as a Paraview file, not the files within the folder

Thank you for your reply, its necessary add this line?

No, it was just to illustrate it when you inspect the boundary visually

I have a question. When I solve in a square domain of 10 x 10 and 30 elements per side, I have this result.

But when I change the domain to a larger one, for example a renctangle of 100 x 100 and 200 elements per side, I have some glyphs.

Is the reason for this difference because paraview plotter?
I have values in this border even if it doesn’t look like it?