How to find the degrees of freedom on the boundary

General
1

Hi everyone,

I want to find the degrees of freedom on the boundary, and I have tried the following method.

from dolfin import * 
# dolfin.__version__ : 2019.2.0.dev0
import numpy as np

mesh = RectangleMesh(Point(0, 0), Point(1, 1), 2, 2)
P = FiniteElement("DG", triangle, 1)
V = FunctionSpace(mesh, P)
bc = DirichletBC(V, Constant(0.0) , "on_boundary")
u = Function(V)
u.vector()[:] = 1.0
bc.apply(u.vector())
boundary_dofs = np.where(u.vector()[:] == 0.0)
print(boundary_dofs)

But this method doesnā€™t seem correct

output:
boundary_dofs =  (array([], dtype=int64),)

How should I modify my code, or is there another method I can use?

Thanks!

2

A DG space has all DOFs in the interior, and Dirichlet BC are typically imposed weakly.

3

I see, thank you!

Sorry, may I ask another question? So, if I want to solve the linear system Ax=b, for example A=assemble(inner(grad(u),grad(v))*dx), b=assemble(f*v*dxāˆ’g*v*ds), is it unnecessary to apply bc.apply(A) and bc.apply(b)?

4

Because I want to add them respectively to a PETSc matrix/vector that I converted from NumPy, I need to solve the DG system using a linear system.

5

If you do not add any boundary conditions, you will have a singular system.

If you want to solve Poisson with pure Neumann, this is possible, see for instance
https://fenicsproject.org/olddocs/dolfin/latest/python/demos/singular-poisson/demo_singular-poisson.py.html
and
https://fenicsproject.org/olddocs/dolfin/latest/python/demos/neumann-poisson/demo_neumann-poisson.py.html

But please do note that you do not have a proper variational formulation for a DG system, as you are missing all the coupling terms between the cells.

1 Like
6

Thank you for your reply.

Sorry, I gave a not-so-good example. My question is, in the case of only Dirichlet boundary uDā€‹, for example

A=assemble(inner(grad(u),grad(v))*dx+uD*u*v*ds), b=assemble(f*v*dx+inner(uD*n,v*n)*ds)

This is my simplified question.

If itā€™s just the above example, I can directly use bc.apply() and successfully solve Ax = b with solve(A, sol.vector(), b). However, since I want to add them respectively to a PETSc matrix/vector that I converted from NumPy, using bc.apply() causes an error (perhaps there is a structural issue with the NumPy array I constructed). Thatā€™s why Iā€™m asking whether bc.apply() is needed in this situation.

7

I donā€™t know why you want the term:

in your problem.

Secondly, you need to apply Dirichlet bcs if you want to solve the correct problem.

Without a reproducible example of this behavior we can unfortunately not help you with this. There shouldnā€™t be anything particular about applying bcs.

The process of applying DirichletBCs in a symmetric/unsymmetric fashion is explained in
http://jsdokken.com/FEniCS23-tutorial/src/lifting.html

8

I am currently trying to solve a minimization problem, and this is the paper I am referring to. Based on equation (3.14) from the paper, I attempted to store the first term on left-hand side and the right-hand side using numpy.

Below is my code. However, after summing the first term with the other terms, I am not sure how to apply the boundary conditions, which is the issue I am facing.

Or maybe the way I am storing things in numpy is incorrect.
Although the program can run and it seems fine in the case of p=2 (the linear case).

from fenics import *
from dolfin import *
from petsc4py import PETSc
import numpy as np
import matplotlib.pyplot as plt
import math

p = 4 # p-Laplace
pdeg = 1
size = 4

# p=2
mesh = RectangleMesh(Point(0, 0), Point(1, 1), size, size)
f = Expression(" 2*pow(pi,2)*sin(pi*x[0])*sin(pi*x[1])", degree=5,pi=math.pi)
uex = Expression(" sin(pi*x[0])*sin(pi*x[1])",degree=5,pi=math.pi)
qex = Expression((" pi*sin(pi*x[1])*cos(pi*x[0]) "," pi*sin(pi*x[0])*cos(pi*x[1]) "),degree=5,pi=math.pi)
sigmaex = Expression((" pi*sin(pi*x[1])*cos(pi*x[0]) "," pi*sin(pi*x[0])*cos(pi*x[1]) "),degree=5,pi=math.pi)
Coe = Constant(1.0)

# p=4
mesh = RectangleMesh(Point(0, 0), Point(1, 1), size, size)
uex = Expression(" sin(pi*x[0])*sin(pi*x[1])",degree=5,pi=math.pi)
qex = Expression((" pi*sin(pi*x[1])*cos(pi*x[0]) "," pi*sin(pi*x[0])*cos(pi*x[1]) "),degree=5,pi=math.pi)
Coe = Expression("pow((pow(pi,2)*pow(sin(pi*x[0]),2)*pow(cos(pi*x[1]),2) + pow(pi,2)*pow(sin(pi*x[1]),2)*pow(cos(pi*x[0]),2)),1.0)", degree=5,pi=math.pi)
sigmaex = Expression((" Coe*pi*sin(pi*x[1])*cos(pi*x[0]) "," Coe*pi*sin(pi*x[0])*cos(pi*x[1]) "),degree=5,pi=math.pi,Coe=Coe)

f = Expression(" (0.5*pow(pi,4.0)*(-cos(pi*(2*x[0] - 2*x[1])) - cos(pi*(2*x[0] + 2*x[1])) + 2) - 2.0*pow(pi,4)*pow(cos(pi*x[0]),2)*cos(2*pi*x[1]) - 2.0*pow(pi,4)*cos(2*pi*x[0])*pow(cos(pi*x[1]),2))*sin(pi*x[0])*sin(pi*x[1])", degree=5,pi=math.pi)

P1 = VectorElement("DG", triangle, pdeg)
P2 = VectorElement("DG", triangle, pdeg)
P3 = FiniteElement("DG", triangle, pdeg)
P1P2P3 = MixedElement([P1, P2, P3])
W = FunctionSpace(mesh, P1P2P3)

V1 = FunctionSpace(mesh, P1)
V2 = FunctionSpace(mesh, P2)
V3 = FunctionSpace(mesh, P3)

(q,sigma,u) = TrialFunctions(W)
(zeta,tau,v) = TestFunctions(W)

n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2

# norm
def vecnorm(u):
  absu = pow(inner(u,u),0.5)
  return absu

def scalarnorm(u):
  absu = pow(pow(u,2),0.5)
  return absu


# DG discrete weak gradient operators
def Discrete_weak_gradient(v,g):
  (q,sigma,uu) = TrialFunctions(W)
  (zeta,tau,vv) = TestFunctions(W)
  v = project(v,V3) # v function
  g = project(g,V3) # g function

  # # choose C12
  C12 = Constant((1.0,0))
  dot_product = inner(C12,n('+'))
  sign = conditional(dot_product > 0, 1.0, -1.0)

  bc = DirichletBC(W.sub(2), g, "on_boundary")

  L = inner(q,zeta)*dx

  R = -v*div(zeta)*dx + avg(v)*inner(jump(zeta),n('+'))*dS\
    + inner(0.5*sign,jump(v))*inner(jump(zeta),n('+'))*dS + g*inner(zeta,n)*ds

  F = L - R
  L = lhs(F)
  R = rhs(F)

  w = Function(W)
  solve(L==R,w,bc)
  (q,sigma,u) = w.split()

  return q

# energy functional
def J(u,g,eta):
  eta = Constant(eta)
  Dwg_u_g = Discrete_weak_gradient(u,g)
  Ju =  1/p*pow(vecnorm(Dwg_u_g),p)*dx\
     + 1/p*eta*pow(h_avg,1-p)*pow(vecnorm(jump(u,n)),p)*dS\
     + 1/p*eta*pow(h,1-p)*pow(scalarnorm(u-g),p)*ds - f*u*dx

  return Ju

# nonlinear coefficient
def nonlinear(vec):
  nonlinear = pow((vecnorm(vec)),p-2)*vec
  return nonlinear


# calculate derivative of energy functional
def dJudv_value(u,g,v,eta,Dwg1,Dwg2): # v function
  eta = Constant(eta)
  nonlinear1 = nonlinear(Dwg1)
  nonlinear2 = nonlinear((1/h_avg)*jump(u,n))
  nonlinear3 = nonlinear((1/h)*(u-g)*n)

  dJudv_value = inner(nonlinear1,Dwg2)*dx\
            + eta*inner(nonlinear2,jump(v,n))*dS\
            + eta*inner(nonlinear3,v*n)*ds - f*v*dx

  return assemble(dJudv_value)


# derivative of energy functional
def dJudv_fn(u,g,v,eta): # v testfunction
  eta = Constant(eta)
  nonlinear2 = nonlinear((1/h_avg)*jump(u,n))
  nonlinear3 = nonlinear((1/h)*(u-g)*n)

  dJudv_fn =  eta*inner(nonlinear2,jump(v,n))*dS\
            + eta*inner(nonlinear3,v*n)*ds - f*v*dx

  return dJudv_fn

Amat = np.zeros((V3.dim(), V3.dim()))
bvec = np.zeros(V3.dim())
dofmap = V3.dofmap()
global_dofs = dofmap.dofs()
Dwg_phi_j_0_list = []
for j in range(V3.dim()):
    phi_j = Function(V3)
    phi_j.vector()[j] = 1
    Dwg_phi_j_0_list.append(Discrete_weak_gradient(phi_j, Constant(0.0)))

print("Dimension V3: ",V3.dim())

# set boundary condition and collect boundary dofs
uD = project(uD,V3)
bc = DirichletBC(V3, uD , "on_boundary")
# boundary_dofs = list(set(bc.get_boundary_values().keys()))


eps = 1.0
tol = 1.0e-6
iter = 0
maxiter = 10
epsi = 1.0e-14 # avoid degeneracy

while eps > tol and iter < maxiter:

  w = TrialFunction(V3)
  v = TestFunction(V3)
  alpha = 1.0

# Find the steepest direction
# The first term of LHS and RHS for p=4 (equation (3.14) )

  Dwg_u0_g = Discrete_weak_gradient(u0,uD)
  Dwg_u0_g_coe = nonlinear(Dwg_u0_g)
  Dwg_u0_g_norm = pow(vecnorm(Dwg_u0_g),p-2)

  for i in range(V3.dim()):
    Dwg_phi_i_0 = Dwg_phi_j_0_list[i]
    phivalue = assemble(inner(Dwg_u0_g_coe, Dwg_phi_i_0) * dx)
    bvec[global_dofs[i]] = phivalue

    for j in range(V3.dim()):
      Dwg_phi_j_0 = Dwg_phi_j_0_list[j]
      phi_ij_value = assemble((epsi + Dwg_u0_g_norm) * inner(Dwg_phi_i_0, Dwg_phi_j_0) * dx)
      Amat[global_dofs[i], global_dofs[j]] = phi_ij_value


  weighted = assemble(Constant(eta)*(1/h_avg) * (epsi +  pow(vecnorm((1/h_avg)*jump(u0,n)),p-2) ) *inner(jump(w,n),jump(v,n))*dS\
            + Constant(eta)*(1/h)* (epsi + pow(vecnorm((1/h)*(u0-uD)),p-2) ) * w * v * ds)
  Apetsc = as_backend_type(weighted).mat()
  Apetsc = Apetsc.getValues(range(0, Apetsc.getSize()[0]), range(0,  Apetsc.getSize()[1]))
  Amat = Amat + Apetsc

  dJudv = assemble(dJudv_fn(u0,uD,v,eta))
  dJudv_vec = dJudv.get_local()
  bvec = bvec + dJudv_vec

 # for dof in boundary_dofs:
   # Amat[dof] = 0
   # Amat[dof, dof] = 1 
   # bvec[dof] = 0


# convert numpy matrix to dolfin.cpp.la.PETScMatrix
  AMat = PETSc.Mat().createAIJ(Amat.shape)
  AMat.setUp()
  AMat.setValues(range(0, Amat.shape[0]), range(0, Amat.shape[1]), Amat)
  AMat.assemble()
  AMat = PETScMatrix(AMat)

# convert numpy vector to dolfin.cpp.la.PETScVector
  bVec = PETSc.Vec().createSeq(V3.dim())
  bVec.setValues(range(0,V3.dim()), bvec)
  bVec = PETScVector(bVec) 

  w = Function(V3)
  solve(AMat, w.vector(), -bVec)


# doing Backtrack line search

  tol_l = 1e-4
  c2 = 1e-4
  c1 = 0.6
  stepsize = alpha
  starting = u0
  direction = w

  Jw_value = assemble(J(w,uD,eta))
  Jz0_value = assemble(J(u0,uD,eta))
  Dwg_u0_g = Discrete_weak_gradient(direction,uD)
  Dwg_dire_0 = Discrete_weak_gradient(direction,Constant(0.0))
  dJudz0 = dJudv_value(starting,uD,direction,eta,Dwg_u0_g,Dwg_dire_0)

  zn = Function(V3)
  iter_l = 0

  while stepsize>tol_l:
    print("****line search iter=%d: stepsize=%.16f ****"%(iter_l,stepsize))
    iter_l  = iter_l + 1
    zn = starting+stepsize*direction
    Jzn_value = assemble(J(zn,uD,eta))
    if Jzn_value < Jz0_value+c2*stepsize*dJudz0:
      break
    else:
      stepsize = c1*stepsize

  print("After doing line search: ")
  print("iter_stepsize=%d:" %iter_l)
  print ("stepsize=%.16f Jz0=%.16f Jw=%.16f Jzn=%.16f Juex=%.16f" %(stepsize, Jz0_value, Jw_value, Jzn_value, Juex_value))

  alpha = stepsize
  u = u0 + alpha * w

  u_er = project(w, V3)
  u_exer = project(u-uex, V3)
  u_L2iter = sqrt(assemble(inner(u_er,u_er)*dx))
  uex_L2iter = sqrt(assemble(inner(u_exer,u_exer)*dx))
  eps = u_L2iter
  eps_ex =  uex_L2iter
  print ("iter=%d: error_iter=%.16f error_ex=%.16f" %(iter, eps, eps_ex))
  print (" ")

  u0 = u
  iter = iter + 1


# error computation

plt.figure()
c = plot(u)
plt.colorbar(c)
plt.title('u DG solution')
plt.legend()
plt.show()

q = project(grad(u), V1)
sigma = pow(sqrt(inner(grad(u), grad(u))),p-2)*grad(u)
sigma = project(sigma, V2)
u = project(u, V3)

V1 = FunctionSpace(mesh, P1); V2 = FunctionSpace(mesh, P2); V3 = FunctionSpace(mesh, P3)
qer = project(q-qex, V1)
sigmaer = project(sigma-sigmaex, V2)
uer = project(u-uex, V3)

q_L2er = sqrt( assemble(inner(qer,qer)*dx) )
sigma_L2er = sqrt( assemble(inner(sigmaer,sigmaer)*dx) )
u_L2er = sqrt( assemble(uer*uer*dx) )

print ("\n")
print("p = %d " %p)
print("u : degree %d , q : degree %d, sigma : degree %d" %(pdeg,pdeg,pdeg))
print("mesh = %d x %d x 2" %(size,size))
print ("|u-uh|_{L2} \t |q-qh|_{L2} \t |sigma-sigmah|_{L2}")
print ("%.16f \t  %.16f \t  %.16f" % (u_L2er, q_L2er, sigma_L2er))
print ("\n")
9

This example is too lengthy for me to have time to look at its entirety.

Have you first tried to solve a simpler problem, like the equivalent Poisson equation, and made sure that you understand how to apply boundary conditions to that problem?

10

Yes, I have tried solving the Poisson equation using CG, and it seems that this method has no issues.

from dolfin import *
import numpy as np
from petsc4py import PETSc
from scipy.sparse import csr_matrix

mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)

f = Expression("2*pow(pi, 2)*sin(pi*x[0])*sin(pi*x[1])", degree=5, pi=np.pi)
uex = Expression("sin(pi*x[0])*sin(pi*x[1])", degree=5, pi=np.pi)
uD = project(uex, V)
bc = DirichletBC(V, uD, "on_boundary")

L = f * v * dx
b = assemble(L)
b = as_backend_type(b)
bc.apply(b)

Amat = np.zeros((V.dim(), V.dim()))
dofmap = V.dofmap() 
global_dofs = dofmap.dofs()


# The term inner(grad(u), grad(v))
for i in range(V.dim()):
    phi_i = Function(V)
    phi_i.vector()[i] = 1
    for j in range(V.dim()):
        phi_j = Function(V)
        phi_j.vector()[j] = 1
        phijvec = assemble(inner(grad(phi_i), grad(phi_j)) * dx)
        Amat[global_dofs[i], global_dofs[j]] = phijvec


boundary_dofs = list(bc.get_boundary_values().keys())
print(boundary_dofs)
for dof in boundary_dofs:
  Amat[dof] = 0
  Amat[dof, dof] = 1


AMat = PETSc.Mat().createAIJ(Amat.shape)
AMat.setUp()
AMat.setValues(range(0, Amat.shape[0]), range(0, Amat.shape[1]), Amat)
AMat.assemble()
AMat = PETScMatrix(AMat)


sol = Function(V)
solve(AMat, sol.vector(), b)

uer = project(uex-sol,V)
u_L2er = sqrt( assemble(uer*uer*dx) )

print("u_L2 error:",u_L2er)
plot(sol)

But I encountered a problem applying the boundary conditions when using DG.

from fenics import *
from dolfin import *
from petsc4py import PETSc
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt


mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, "DG", 1)
u = TrialFunction(V)
v = TestFunction(V)
alpha = Constant(10.0)
n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2

f = Expression("2*pow(pi, 2)*sin(pi*x[0])*sin(pi*x[1])", degree=5, pi=np.pi)
uex = Expression("sin(pi*x[0])*sin(pi*x[1])", degree=5, pi=np.pi)
uD = project(uex, V)
bc = DirichletBC(V, uD ,"on_boundary")


Amat = np.zeros((V.dim(), V.dim()))
bvec = np.zeros(V.dim())
dofmap = V.dofmap() 
global_dofs = dofmap.dofs()


# The first term of a and L    
# dot( grad (u), grad (v) )*dx , f*v*dx
for i in range(V.dim()):
    phi_i = Function(V)
    phi_i.vector()[i] = 1
    phivalue = assemble(f * phi_i * dx)
    bvec[global_dofs[i]] = phivalue
    for j in range(V.dim()):
        phi_j = Function(V)
        phi_j.vector()[j] = 1
        phi_ij = assemble(inner(grad(phi_i), grad(phi_j)) * dx)
        Amat[global_dofs[i], global_dofs[j]] = phi_ij

a = assemble(-dot(avg(grad(u)),jump(v,n))*dS-dot(avg(grad(v)),jump(u,n))*dS\
    +alpha / h_avg *dot( jump (v, n), jump (u, n))*dS\
    -dot(grad(u),v*n)*ds-dot(grad(v),u*n)*ds\
    +alpha /h*u*v*ds )

L = assemble(-dot(grad(v),uD*n)*ds\
    +alpha /h*uD*v*ds )


Apetsc = as_backend_type(a).mat()
Apetsc = Apetsc.getValues(range(0, Apetsc.getSize()[0]), range(0,  Apetsc.getSize()[1]))
Amat = Amat - Apetsc


L_vec = L.get_local()
bvec = bvec - L_vec

# boundary_dofs = list(bc.get_boundary_values().keys())
# for dof in boundary_dofs:
#   Amat[dof] = 0
#   Amat[dof, dof] = 1

AMat = PETSc.Mat().createAIJ(Amat.shape)
AMat.setUp()
AMat.setValues(range(0, Amat.shape[0]), range(0, Amat.shape[1]), Amat)
AMat.assemble()
AMat = PETScMatrix(AMat)


# convert numpy matrix to dolfin.cpp.la.PETScMatrix
AMat = PETSc.Mat().createAIJ(Amat.shape)
AMat.setUp()
AMat.setValues(range(0, Amat.shape[0]), range(0, Amat.shape[1]), Amat)
AMat.assemble()
AMat = PETScMatrix(AMat)

# convert numpy vector to dolfin.cpp.la.PETScVector
bVec = PETSc.Vec().createSeq(V.dim())
bVec.setValues(range(0,V.dim()), bvec)
bVec = PETScVector(bVec) 

sol = Function(V)
solve(AMat, sol.vector(), bVec)

uer = project(uex-sol,V)
u_L2er = sqrt( assemble(uer*uer*dx) )

print("u_L2 error:",u_L2er)
plot(sol)

Thatā€™s also why I wanted to ask how to find the degrees of freedom on the boundary when using DG.

11

Havenā€™t we been through this already?
Strong Dirichlet conditions does not make sense for DG.
You can see this with:

from dolfin import *
mesh = UnitSquareMesh(5, 5)
V = FunctionSpace(mesh, "CG", 1)

uh = Function(V)

u_ex = Expression("x[0]*x[0] + x[1]*x[1]", degree=2)
bc = DirichletBC(V, u_ex ,"on_boundary")
bc.apply(uh.vector())
with XDMFFile("uh.xdmf") as file:
    file.write_checkpoint(uh, "uh", 0.0, append=False)

which generates
image

while

from dolfin import *
mesh = UnitSquareMesh(5, 5)
V = FunctionSpace(mesh, "DG", 1)

uh = Function(V)

u_ex = Expression("x[0]*x[0] + x[1]*x[1]", degree=2)
bc = DirichletBC(V, u_ex ,"on_boundary")
bc.apply(uh.vector())
with XDMFFile("uh.xdmf") as file:
    file.write_checkpoint(uh, "uh", 0.0, append=False)

generates
image
as ā€œon_boundaryā€ only works for dofs associated with facets/vertices/edges, which no DG dof is.
You could do:

from dolfin import *
mesh = UnitSquareMesh(5, 5)
V = FunctionSpace(mesh, "DG", 1)

uh = Function(V)

u_ex = Expression("x[0]*x[0] + x[1]*x[1]", degree=2)
bc = DirichletBC(V, u_ex ,"near(x[0], 0)| near(x[0], 1) | near(x[1], 0) | near(x[1],1)", method="pointwise")
bc.apply(uh.vector())
with XDMFFile("uh.xdmf") as file:
    file.write_checkpoint(uh, "uh", 0.0, append=False)

resulting in a correct strong enforcement.
However, note as already mentioned, that one usually treats DG boundary conditions weakly with a Nitsche type formulation.

1 Like
12

Thank you for your reply!

Regarding using DG to solve the Poisson problem, I understand that when I treat the problem as a variational problem, the Dirichlet BCs are typically imposed weakly.

However, when I solve it as a linear system, do I need to handle the assembled matrix/vector differently?

If not, I think the issue might lie in how numpy assembles the components.

13

You would need to add the extra terms to your variational form.

As shown above, you donā€™t get Dirichlet enforcement strongly with your current code, leading to it being interpreted as a pure Neumann condition

1 Like
14

This is the result I got when I directly assembled the bilinear form
š‘Ž and the linear form š‘ from the variational problem and solved the linear system, which gave me the desired result.

from fenics import *
from dolfin import *
from petsc4py import PETSc
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt


mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, "DG", 1)
u = TrialFunction(V)
v = TestFunction(V)
alpha = Constant(10.0)
n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2

f = Expression("2*pow(pi, 2)*sin(pi*x[0])*sin(pi*x[1])", degree=5, pi=np.pi)
uex = Expression("sin(pi*x[0])*sin(pi*x[1])", degree=5, pi=np.pi)
uD = project(uex, V)
# bc = DirichletBC(V, uD ,"on_boundary")
# bc = DirichletBC(V, uD ,"near(x[0], 0)| near(x[0], 1) | near(x[1], 0) | near(x[1],1)", method="pointwise")


a = assemble(dot( grad (u), grad (v))*dx\
    -dot(avg(grad(u)),jump(v,n))*dS-dot(avg(grad(v)),jump(u,n))*dS\
    +alpha / h_avg *dot( jump (v, n), jump (u, n))*dS\
    -dot(grad(u),v*n)*ds-dot(grad(v),u*n)*ds\
    +alpha /h*u*v*ds)

L = assemble(f*v*dx\
    -dot(grad(v),uD*n)*ds\
    +alpha /h*uD*v*ds)


sol = Function(V)
solve(a, sol.vector() ,L)

uer = project(uex-sol,V)
u_L2er = sqrt( assemble(uer*uer*dx) )

print("u_L2 error:",u_L2er)

plt.figure()
c = plot(sol)
plt.colorbar(c)
plt.title('u DG solution')
plt.legend()
plt.show()

However, the following is the result when I used the code I provided above, but it produced an issue. I also tried modifying it based on your suggestion, but I still donā€™t know where the problem is.

from fenics import *
from dolfin import *
from petsc4py import PETSc
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt


mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, "DG", 1)
u = TrialFunction(V)
v = TestFunction(V)
alpha = Constant(10.0)
n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2

f = Expression("2*pow(pi, 2)*sin(pi*x[0])*sin(pi*x[1])", degree=5, pi=np.pi)
uex = Expression("sin(pi*x[0])*sin(pi*x[1])", degree=5, pi=np.pi)
uD = project(uex, V)
# bc = DirichletBC(V, uD ,"on_boundary")
# bc = DirichletBC(V, uD ,"near(x[0], 0)| near(x[0], 1) | near(x[1], 0) | near(x[1],1)", method="pointwise")

Amat = np.zeros((V.dim(), V.dim()))
bvec = np.zeros(V.dim())
dofmap = V.dofmap() 
global_dofs = dofmap.dofs()

# The first term of a and L    
# dot( grad (u), grad (v) )*dx , f*v*dx
for i in range(V.dim()):
    phi_i = Function(V)
    phi_i.vector()[i] = 1
    phivalue = assemble(f * phi_i * dx)
    bvec[global_dofs[i]] = phivalue
    for j in range(V.dim()):
        phi_j = Function(V)
        phi_j.vector()[j] = 1
        phi_ij = assemble(inner(grad(phi_i), grad(phi_j)) * dx)
        Amat[global_dofs[i], global_dofs[j]] = phi_ij

a = assemble(-dot(avg(grad(u)),jump(v,n))*dS-dot(avg(grad(v)),jump(u,n))*dS\
    +alpha / h_avg *dot( jump (v, n), jump (u, n))*dS\
    -dot(grad(u),v*n)*ds-dot(grad(v),u*n)*ds\
    +alpha /h*u*v*ds )

L = assemble(-dot(grad(v),uD*n)*ds\
    +alpha /h*uD*v*ds )


Apetsc = as_backend_type(a).mat()
Apetsc = Apetsc.getValues(range(0, Apetsc.getSize()[0]), range(0,  Apetsc.getSize()[1]))
Amat = Amat - Apetsc

L_vec = L.get_local()
bvec = bvec - L_vec

# boundary_dofs = list(bc.get_boundary_values().keys())
# for dof in boundary_dofs:
#   Amat[dof] = 0
#   Amat[dof, dof] = 1.0
#   bvec[dof] = 0

# convert numpy matrix to dolfin.cpp.la.PETScMatrix
AMat = PETSc.Mat().createAIJ(Amat.shape)
AMat.setUp()
AMat.setValues(range(0, Amat.shape[0]), range(0, Amat.shape[1]), Amat)
AMat.assemble()
AMat = PETScMatrix(AMat)

# convert numpy vector to dolfin.cpp.la.PETScVector
bVec = PETSc.Vec().createSeq(V.dim())
bVec.setValues(range(0,V.dim()), bvec)
bVec = PETScVector(bVec) 

sol = Function(V)
solve(AMat, sol.vector(), bVec)

uer = project(uex-sol,V)
u_L2er = sqrt( assemble(uer*uer*dx) )

print("u_L2 error:",u_L2er)

plt.figure()
c = plot(sol)
plt.colorbar(c)
plt.title('u DG solution')
plt.legend()
plt.show()
15

Sorry, I found that the issue was that I got the signs wrong.
Amat = Amat + Apetsc
bvec = bvec + L_vec

Thank you for your help!