3D boundary conditions in mixed formulation

dolfinx
1

Hi,

I’m setting a displacement component (u_x=0.1) on the front surface(x=1) for a unit cubic domain to test boundary conditions

V = fem.FunctionSpace(domain, element_mix)
V_disp = V.sub(0)
V_disp_x = V_disp.sub(0)

def front(x):
    return np.isclose(x[0], 1, atol=1.0e-12)

and the u_x on front surface is

u_x =  fem.Constant(domain, 0.1)
def expression(x):
    return np.full(x.shape[1], u_x)

fcn_x = fem.Function(V_disp_x_collapse)
fcn_x.interpolate(expression)
dofs_front_x = fem.locate_dofs_geometrical((V_disp_x, V_disp_x_collapse), front)
bc_front_x = fem.dirichletbc(fcn_x, dofs_front_x, V_disp_x)

However, when I plot the boundary condition, it looks like:

image
image1284×1086 49 KB

Thanks

2

There’s not enough context for anyone to be able to help:

  • the code is just a small snippet
  • the image does not have the coordinate axis, how are we supposed to know which side of the square is aligned with x, y and z?
  • the image does not have a color bar, hence we are left to imagine by our own what red and blue mean (0 and 0.1, one can guess)

It would helpful if you provide a SMALL (but complete) code showing the issue. There are several posts in this forum about BCs with mixed elements, but (to my memory) this one is slightly different because it is a component of a mixed element, and having a complete code could be helpful for people finding this post in future with the search button.

3

Here’s the complete code:

import numpy as np
import matplotlib.pyplot as plt

from dolfinx import mesh, fem, io, plot, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
import ufl  

import gmsh
import os
from tqdm import tqdm
import csv
import sys

### (1a) define the 3D domain
coords_lower_left = [0.0, 0.0, 0.0]
Lx, Ly, Lz = 1.0, 1.0, 1.0
coords_upper_right = [Lx, Ly, Lz]
num_elements = [1, 1, 1]
domain = mesh.create_box(MPI.COMM_WORLD, [coords_lower_left, coords_upper_right], num_elements, mesh.CellType.hexahedron)

### (1b) set element types 
element_u = ufl.VectorElement("CG", domain.ufl_cell(), degree=2) ## vector: displacement
element_p= ufl.FiniteElement("CG", domain.ufl_cell(), degree=1) ## scalar: pressure
element_mix = ufl.MixedElement(element_u, element_p)

V = fem.FunctionSpace(domain, element_mix)
V_disp = V.sub(0)
V_press = V.sub(1)

### displacement component in x direction
V_disp_x = V_disp.sub(0)

def front(x):
    return np.isclose(x[0], 1, atol=1.0e-12)

fdim = domain.topology.dim - 1
value_zero=  fem.Constant(domain, 0.1)
facet_front = mesh.locate_entities_boundary(domain, fdim, front)
dofs_front_x = fem.locate_dofs_topological(V_disp_x, fdim, facet_front)
bc_front_x = fem.dirichletbc(value_zero, dofs_front_x, V_disp_x)

bcs = [bc_front_x]

I tried “fem.locate_dofs_topological” in the above to get dofs.

For plot the figure in Paraview, the code is

##### plot boundary conditions to Paraview for mixed formulation
Mixed_Space = V
u_boundary = fem.Function(Mixed_Space, name='boundary')
fem.petsc.set_bc(u_boundary.vector, bcs)
u_boundary.x.scatter_forward()

with io.XDMFFile(domain.comm, "boundary/bc_geometrical.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(u_boundary.sub(0))
5

That’s not the complete code: at the very least the are missing imports. It needs to be something that anyone can copy and paste in order to run on their own.

6

See also Dirichletbc not working as expected for slip condition - #2 by francesco-ballarin

7

Hi Francesco, I complement the imports in the above code

8

The mistake comes from writing out the bc to xdmf file, it should be

state = fem.Function(V)
fem.petsc.set_bc(state.vector, bcs)

with io.XDMFFile(domain.comm, "boundary/test_mixed_bc_topo.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(state.sub(0).collapse())

rather than

xdmf.write_function(state.sub(0))