Dirichlet boundary of DG/RT/N1curl element?

General
1

Hi everyone!
I attempted to assemble Dirichlet boundary conditions for DG finite element analysis, but in my testing, dofs was empty. It seems that it cannot locate the freedom of the node. Why is this?

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
from petsc4py import PETSc as _PETSc
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, 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()

(unv,Anv,Pn)=ufl.TrialFunctions(W)   #
(v,fai,q)=ufl.TestFunctions(W)       #

un=fem.Function(W0)                   

def initial_u(x):
    values=np.zeros((3,x.shape[1]))
    values[0,:]=2    #y
    values[1,:]=2    #z
    values[2,:]=2    #x
    return values

un.interpolate(initial_u)  

An=fem.Function(W1)
def initial_A(x):
    values=np.zeros((3,x.shape[1]))
    values[0,:]=3       #z
    values[1,:]=3       #0
    values[2,:]=3       #y
    return values
An.interpolate(initial_A)

P=fem.Function(W2)
def initial_P(x):
    values=np.zeros((x.shape[1]))
    return values
P.interpolate(initial_P)

#---------------------------------------------------------------------------------
def local_u(x):
    return np.isclose(x[1],0)

def u_D(x):
    values=np.zeros((3,x.shape[1]))
    values[0,:]=0           
    values[1,:]=0           
    values[2,:]=0           
    return values

Vh=fem.FunctionSpace(msh,u_element)
u_boundary = Function(W0)
u_boundary.interpolate(u_D)
facets = locate_entities_boundary(msh,2, local_u)
dofs = locate_dofs_topological((W.sub(0), W0), 2, facets)
bc0 = dirichletbc(u_boundary, dofs, W.sub(0))

print('facets is',facets)
bcs=[bc0]
print('un initial is',un.x.array)
set_bc(un.vector,bcs)
print('un last is',un.x.array)
print('dofs is',dofs)

Do I have any other mandatory measures besides adding a penalty item in the weak form?
Many Thanks!
If there is relevant information, it would be even better!

1 Like
2

A DG element does not have any dofs associated with its, facets, thus, locate_dofs_topological will always be empty.
A dof is associated with a facet would mean that the degree of freedom is shared with all other cells that has this facet. If you instead use locate_dofs_geometrical:

import numpy as np
import ufl
from dolfinx import fem
from dolfinx.fem import (Function, dirichletbc,
                         locate_dofs_geometrical)
from dolfinx.mesh import (CellType, GhostMode, create_unit_cube)
from mpi4py import MPI
from dolfinx.fem.petsc import (set_bc)
import ufl
from petsc4py import PETSc as _PETSc

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 local_u(x):
    return np.isclose(x[1], 0)


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


u_boundary = Function(W0)
u_boundary.interpolate(u_D)
dofs = locate_dofs_geometrical((W.sub(0), W0), local_u)
bc0 = dirichletbc(u_boundary, dofs, W.sub(0))

bcs = [bc0]
un = Function(W)
print(f'un initial is {un.x.array}')
set_bc(un.vector, bcs)
print(f'un last is {un.x.array}')
print(f'{dofs=}')
2 Likes