The Problem of Dirichlet Boundary Conditions for DG

Hi everyone! Dear community!
I specify the Dirichlet boundary condition in the mixed element.

msh = create_unit_cube(MPI.COMM_WORLD,
                      3,3,3,
                       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)      

define Dirichlet boundary:

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)
def u_D(x):
    values = np.zeros_like(x, dtype=_PETSc.ScalarType)
    values[0, :] = 1
    values[1, :] = 0
    values[2, :] = 0
    return values

set the bcs,have two method for DG.
Method 1

facets = locate_entities_boundary(msh,2,boundarymark)
dofs1 = locate_dofs_topological((W.sub(0), W0), 2, facets)
bc1 = dirichletbc(un_D, dofs1, W.sub(0))
bcs=[bc1]

Method 2

dofs0 = locate_dofs_geometrical((W0, W0), boundarymark)
bc0 = dirichletbc(un_D, dofs0, W0)
bcs=[bc0]
print('dofs0 is',dofs0)

The method 1 is from many official demo :
(Divergence conforming discontinuous Galerkin method for the Navier–Stokes equations — DOLFINx 0.8.0.0 documentation)
For mixed elements, this is how they are implemented, but the difference is that no one uses DG.
In my testing,Method 2 is correct. It will make the boundary of the output mandatory.
So can someone explain to me why this is?
Method 1and 2,What kind of step is it.
The main code is :

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,
                       5, 5, 5,
                       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, :] = 0
    values[2, :] = 0
    return values

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

u_boundary = Function(W0)
un_D=fem.Function(W0)   
un_D.interpolate(u_D)
An_D=fem.Function(W1)
An_D.interpolate(A_D)

(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)
print('dofs0 is',dofs0)

facets = locate_entities_boundary(msh,2,boundarymark)
dofs1 = locate_dofs_topological((W.sub(0), W0), 2, facets)
bc1 = dirichletbc(un_D, dofs1, W.sub(0))
print('dofs1 is',dofs1)      #dofs1 is empty    It's means bc1 is empty

As I explained in: Dirichlet boundary of DG/RT/N1curl element? - #2 by dokken
locate_dofs_topological works on dofs associated with the a set of entities (say facets in your case).
A DG element has all its dofs associated with the cell, not the facets or vertices, as they are not shared with neighbouring elements.

locate_dofs_geometrical tabulates the coordinates of the degrees of freedom (their location in physical space if they are defined as point evaluations) and check if they satisfies the conditional sent in.

2 Likes

Thanks dokken ! I got it! :smile:

By the way, is only DG not shared ? Or rather, only DG is different from other Finitre elements?

Any Finite element space that is represented as discontinuous (called broken in mathematical textbooks) consist of degrees of freedom that is only represented on a single element.
This means that no degrees of freedom are shared with neighbouring cells.

The DG (discontinuous Galerkin) space is an example of the «broken» Lagrange/P space.