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