Hcurl and DG in mixed element assemble the boundary conditions

General
1

Hi everyone! Dear community！
I am trying to learn the boundary condition processing of DG and curl elements.
Because I haven’t found it in any tutorials.Progress is very slow.
First I created a mesh and defined a mixed finite element space.

msh = create_unit_cube(MPI.COMM_WORLD,
                       1, 1, 1,
                       CellType.tetrahedron, ghost_mode=GhostMode.none)

u_element = ufl.VectorElement("DG", ufl.tetrahedron, 1)
A_element = ufl.FiniteElement("N1curl", ufl.tetrahedron, 1)
P_element = ufl.FiniteElement("DG", ufl.tetrahedron, 0)

TH = ufl.MixedElement([u_element, A_element, P_element])
W = fem.FunctionSpace(msh, TH)

W0, map0 = W.sub(0).collapse()
W1, map1 = W.sub(1).collapse()
W2, map2 = W.sub(2).collapse()

def u_D(x):
    values = np.zeros_like(x, dtype=_PETSc.ScalarType)
    values[0, :] = 1
    values[1, :] = 2
    values[2, :] = 3
    return values

def A_D(x):
    values = np.zeros_like(x, dtype=_PETSc.ScalarType)
    values[0, :] = 3
    values[1, :] = 2
    values[2, :] = 1
    return values

u_D and A_D is Dirichlet boundary condition that I want.
Second I mark the boundary and assemble the boundary conditions

un_D=fem.Function(W0)   
un_D.interpolate(u_D)
un=fem.Function(W0)  
An_D=fem.Function(W1)
An_D.interpolate(A_D)
An=fem.Function(W1)
(unv,Anv,Pn)=ufl.TrialFunctions(W)    #
(v,fai,q)=ufl.TestFunctions(W)        #

def boundarymark(x):
    return np.isclose(x[0],0)|np.isclose(x[0],1)|np.isclose(x[1],0)|np.isclose(x[1],1)|np.isclose(x[2],0)|np.isclose(x[2],1)

dofs0 = locate_dofs_geometrical((W0, W0), boundarymark)    
bc0 = dirichletbc(un_D, dofs0, W0)    
bcsu=[bc0]

dofs1 = locate_dofs_geometrical((W0, W0), boundarymark)   
bc1 = dirichletbc(An_D, dofs1, W.sub(1)) 
bcsa=[bc1]

set_bc(un.vector,bcsu)
set_bc(An.vector,bcsa)
The output is :Un

image
image761×264 3.51 KB

An:
image
This code seems different from the Dolfinx tutorial.
Is this correct?
**The most important thing is that bcsa and bcsu the last thing I need to use when using Dirichlet boundary conditions?**i.e. apply_lifting(b,[a],[bcs]) ?

2

The main code:

import numpy as np
import ufl
from dolfinx import cpp as _cpp
from dolfinx import fem
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc,
                         extract_function_spaces, form,
                         locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_rectangle,
                          locate_entities_boundary,create_unit_cube)
from ufl import div, dx, grad, inner,dot,outer,ds,dS,avg,jump,rot,cross
from mpi4py import MPI
from petsc4py import PETSc as _PETSc
from dolfinx.io import XDMFFile
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, 
                               create_vector, set_bc)

msh = create_unit_cube(MPI.COMM_WORLD,
                       1, 1, 1,
                       CellType.tetrahedron, ghost_mode=GhostMode.none)

u_element = ufl.VectorElement("DG", ufl.tetrahedron, 1)
A_element = ufl.FiniteElement("N1curl", ufl.tetrahedron, 1)
P_element = ufl.FiniteElement("DG", ufl.tetrahedron, 0)

TH = ufl.MixedElement([u_element, A_element, P_element])
W = fem.FunctionSpace(msh, TH)

W0, map0 = W.sub(0).collapse()
W1, map1 = W.sub(1).collapse()
W2, map2 = W.sub(2).collapse()

def u_D(x):
    values = np.zeros_like(x, dtype=_PETSc.ScalarType)
    values[0, :] = 1
    values[1, :] = 2
    values[2, :] = 3
    return values

def A_D(x):
    values = np.zeros_like(x, dtype=_PETSc.ScalarType)
    values[0, :] = 3
    values[1, :] = 2
    values[2, :] = 1
    return values


un_D=fem.Function(W0)   
un_D.interpolate(u_D)
un=fem.Function(W0)  
An_D=fem.Function(W1)
An_D.interpolate(A_D)
An=fem.Function(W1)
(unv,Anv,Pn)=ufl.TrialFunctions(W)    #
(v,fai,q)=ufl.TestFunctions(W)        #

def boundarymark(x):
    return np.isclose(x[0],0)|np.isclose(x[0],1)|np.isclose(x[1],0)|np.isclose(x[1],1)|np.isclose(x[2],0)|np.isclose(x[2],1)

dofs0 = locate_dofs_geometrical((W0, W0), boundarymark)    
bc0 = dirichletbc(un_D, dofs0, W0)    

bcsu=[bc0]

dofs1 = locate_dofs_geometrical((W0, W0), boundarymark)   
bc1 = dirichletbc(An_D, dofs1, W.sub(1)) 

bcsa=[bc1]
set_bc(un.vector,bcsu)
set_bc(An.vector,bcsa)

print('un is',un.x.array)
print('dofs1(A) is',dofs1)
print('An is',An.x.array)

with XDMFFile(MPI.COMM_WORLD, "boundary/un.xdmf", "w") as ufile_xdmf:
    un.x.scatter_forward()
    ufile_xdmf.write_mesh(msh)
    ufile_xdmf.write_function(un)

with XDMFFile(MPI.COMM_WORLD, "boundary/An.xdmf", "w") as pfile_xdmf:
    An.x.scatter_forward()
    pfile_xdmf.write_mesh(msh)
    pfile_xdmf.write_function(An)
3

I would suggest you start by visualizing the un and An vector (for instance in Paraview) to verify if the boundary conditions are set correctly.