I have tried the HHO code for the Poisson equation using quadrilateral mesh discretization; it is running, but not achieving superconvergence. Here is my code>
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve
import dolfinx
import dolfinx.fem as fem
import dolfinx.mesh as mesh
import basix
import ufl
from mpi4py import MPI
import matplotlib.pyplot as plt
import warnings
# ----------------------------------------------------------------------
# Parameters
# ----------------------------------------------------------------------
k = 1 # degree for cell and face unknowns (Pk)
deg_phi = k + 1 # degree for reconstructed polynomial
n_ref = 4 # number of refinement levels
L2error = np.zeros(n_ref)
DGerror = np.zeros(n_ref)
RGerror = np.zeros(n_ref)
Dx = np.zeros(n_ref)
def u_exact(x):
return np.sin(np.pi * x[0]) * np.sin(np.pi * x[1])
# ----------------------------------------------------------------------
# Helper: map reference interval [0,1] to reference quadrilateral edge
# ----------------------------------------------------------------------
def map_interval_to_quad_edge(t, local_edge):
"""
local_edge: 0 -> bottom (v0->v1), 1 -> right (v1->v2),
2 -> top (v2->v3), 3 -> left (v3->v0)
"""
if local_edge == 0: # bottom: (t, 0)
return t, 0.0
elif local_edge == 1: # right: (1, t)
return 1.0, t
elif local_edge == 2: # top: (1-t, 1)
return 1.0 - t, 1.0
else: # left: (0, 1-t)
return 0.0, 1.0 - t
# ----------------------------------------------------------------------
# Helper: safe inverse of 2x2 matrix (with pseudo‑inverse fallback)
# ----------------------------------------------------------------------
def safe_inv_2x2(J, tol=1e-14):
det = J[0,0]*J[1,1] - J[0,1]*J[1,0]
if abs(det) < tol:
# Use pseudo‑inverse to avoid singular matrix error
return np.linalg.pinv(J)
else:
return np.linalg.inv(J)
# ----------------------------------------------------------------------
# Helper: compute physical normal and length for a facet
# ----------------------------------------------------------------------
def facet_normal_and_length(mesh_, facet_index):
f_verts = mesh_.topology.connectivity(mesh_.topology.dim-1, 0).links(facet_index)
coords = mesh_.geometry.x[f_verts][:, :2]
v = coords[1] - coords[0]
length = np.linalg.norm(v)
f_to_c = mesh_.topology.connectivity(mesh_.topology.dim-1, mesh_.topology.dim)
adj = f_to_c.links(facet_index)
cell_index = adj[0]
cell_verts = mesh_.topology.connectivity(mesh_.topology.dim, 0).links(cell_index)
cell_center = np.mean(mesh_.geometry.x[cell_verts][:, :2], axis=0)
facet_center = np.mean(coords, axis=0)
vec = facet_center - cell_center
n = np.array([v[1], -v[0]])
n = n / np.linalg.norm(n)
if np.dot(n, vec) < 0:
n = -n
return n, length
# ----------------------------------------------------------------------
# Main loop over refinement levels
# ----------------------------------------------------------------------
for n in range(n_ref):
print(f"\n--- Level {n}, cells per edge = {2**(n+4)} ---")
n_cells = 2**(n+4)
mesh_ = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0,0.0],[1.0,1.0]],
[n_cells,n_cells], mesh.CellType.quadrilateral)
h = 1.0 / n_cells
Dx[n] = h
mesh_.topology.create_entities(1)
mesh_.topology.create_connectivity(mesh_.topology.dim-1, 0)
mesh_.topology.create_connectivity(mesh_.topology.dim, mesh_.topology.dim-1)
mesh_.topology.create_connectivity(mesh_.topology.dim-1, mesh_.topology.dim)
# ------------------------------------------------------------------
# 1. Function spaces
# ------------------------------------------------------------------
V_T = fem.functionspace(mesh_, ("DG", k))
V_phi = fem.functionspace(mesh_, ("DG", deg_phi))
num_facets = mesh_.topology.index_map(mesh_.topology.dim-1).size_local
facets = np.arange(num_facets, dtype=np.int32)
facet_mesh, facet_map, _, _ = mesh.create_submesh(mesh_, mesh_.topology.dim-1, facets)
V_F = fem.functionspace(facet_mesh, ("DG", k))
map_T = V_T.dofmap
map_phi = V_phi.dofmap
map_F = V_F.dofmap
N_T = map_T.index_map.size_local
N_phi = map_phi.index_map.size_local
N_F = map_F.index_map.size_local
# ------------------------------------------------------------------
# 2. Basix elements and quadrature
# ------------------------------------------------------------------
elem_T = basix.create_element(basix.ElementFamily.P, basix.CellType.quadrilateral, k,
basix.LagrangeVariant.equispaced)
elem_phi = basix.create_element(basix.ElementFamily.P, basix.CellType.quadrilateral, deg_phi,
basix.LagrangeVariant.equispaced)
elem_F = basix.create_element(basix.ElementFamily.P, basix.CellType.interval, k,
basix.LagrangeVariant.equispaced)
q_pts_cell, q_wts_cell = basix.make_quadrature(basix.CellType.quadrilateral, 2*deg_phi)
q_pts_facet, q_wts_facet = basix.make_quadrature(basix.CellType.interval, 2*deg_phi)
tab_T_cell = elem_T.tabulate(1, q_pts_cell)
tab_phi_cell = elem_phi.tabulate(1, q_pts_cell)
phi_T_vals = tab_T_cell[0]
phi_phi_vals = tab_phi_cell[0]
grad_phi_ref = tab_phi_cell[1]
if grad_phi_ref.ndim == 3:
if grad_phi_ref.shape[2] == 1:
warnings.warn("Gradient tabulation returned only one component; duplicating to make 2D.")
grad_phi_ref = np.concatenate([grad_phi_ref, grad_phi_ref], axis=2)
elif grad_phi_ref.shape[2] != 2:
grad_phi_ref = grad_phi_ref[:, :, :2]
else:
raise ValueError("Unexpected gradient shape from Basix tabulation.")
if phi_phi_vals.ndim == 3:
phi_phi_vals = phi_phi_vals[:, :, 0]
# ------------------------------------------------------------------
# 3. Sparse matrices
# ------------------------------------------------------------------
projT = sp.lil_matrix((N_T, N_phi))
projFd = sp.lil_matrix((N_F, N_phi))
potT = sp.lil_matrix((N_phi, N_T))
potF = sp.lil_matrix((N_phi, N_F))
M_TT_facet = sp.lil_matrix((N_T, N_T))
M_TF_facet = sp.lil_matrix((N_T, N_F))
M_FF_facet = sp.lil_matrix((N_F, N_F))
uPhi = ufl.TrialFunction(V_phi)
vPhi = ufl.TestFunction(V_phi)
S_phi_form = ufl.inner(ufl.grad(uPhi), ufl.grad(vPhi)) * ufl.dx
S_phi = fem.assemble_matrix(fem.form(S_phi_form))
S_phi = S_phi.to_scipy().tocsc()
coords = mesh_.geometry.x
f_to_c = mesh_.topology.connectivity(mesh_.topology.dim-1, mesh_.topology.dim)
cell_facets = mesh_.topology.connectivity(mesh_.topology.dim, mesh_.topology.dim-1)
facet_vertices = mesh_.topology.connectivity(mesh_.topology.dim-1, 0)
vertices_cell = mesh_.topology.connectivity(mesh_.topology.dim, 0)
num_cells = mesh_.topology.index_map(mesh_.topology.dim).size_local
num_facets_global = mesh_.topology.index_map(mesh_.topology.dim-1).size_local
# ------------------------------------------------------------------
# 4. Cell loop: projT, potT, potF
# ------------------------------------------------------------------
for c in range(num_cells):
dofs_T = map_T.cell_dofs(c)
dofs_phi = map_phi.cell_dofs(c)
nT = len(dofs_T)
nphi = len(dofs_phi)
cell_verts = vertices_cell.links(c)
cell_coords = coords[cell_verts][:, :2] # (4,2)
# ---- projT ----
M_T_loc = np.zeros((nT, nT))
M_T_phi_loc = np.zeros((nT, nphi))
for q in range(len(q_pts_cell)):
xi, eta = q_pts_cell[q]
wq = q_wts_cell[q]
N = np.array([(1-xi)*(1-eta), xi*(1-eta), xi*eta, (1-xi)*eta])
dN_dxi = np.array([-(1-eta), (1-eta), eta, -eta])
dN_deta = np.array([-(1-xi), -xi, xi, (1-xi)])
J = np.zeros((2,2))
J[0,0] = np.dot(dN_dxi, cell_coords[:,0])
J[0,1] = np.dot(dN_deta, cell_coords[:,0])
J[1,0] = np.dot(dN_dxi, cell_coords[:,1])
J[1,1] = np.dot(dN_deta, cell_coords[:,1])
detJ = np.linalg.det(J)
w = wq * abs(detJ)
phi_T = phi_T_vals[q, :]
phi_phi = phi_phi_vals[q, :]
for i in range(nT):
for j in range(nT):
M_T_loc[i, j] += w * phi_T[i] * phi_T[j]
for j in range(nphi):
M_T_phi_loc[i, j] += w * phi_T[i] * phi_phi[j]
try:
projT_loc = np.linalg.solve(M_T_loc, M_T_phi_loc)
except np.linalg.LinAlgError:
projT_loc = np.zeros((nT, nphi))
for i, gi in enumerate(dofs_T):
for j, gj in enumerate(dofs_phi):
projT[gi, gj] = projT_loc[i, j]
# ---- potT and potF ----
S_phi_loc = np.zeros((nphi, nphi))
RHS_T_loc = np.zeros((nphi, nT))
facet_contribs = []
local_facets = cell_facets.links(c)
for f in local_facets:
adj = f_to_c.links(f)
side_cell = adj[0]
if side_cell != c:
continue
facet_vert_indices = facet_vertices.links(f)
f_coords = coords[facet_vert_indices][:, :2]
phys_length = np.linalg.norm(f_coords[1] - f_coords[0])
# ---- Determine local edge ----
local_edge = -1
# 1) index-based
for e in range(4):
v1_idx = cell_verts[e]
v2_idx = cell_verts[(e+1)%4]
if (v1_idx == facet_vert_indices[0] and v2_idx == facet_vert_indices[1]) or \
(v1_idx == facet_vert_indices[1] and v2_idx == facet_vert_indices[0]):
local_edge = e
break
# 2) coordinate-based with tolerance
if local_edge == -1:
for e in range(4):
v1 = cell_coords[e]
v2 = cell_coords[(e+1)%4]
if (np.allclose(v1, f_coords[0], rtol=1e-6, atol=1e-8) and
np.allclose(v2, f_coords[1], rtol=1e-6, atol=1e-8)) or \
(np.allclose(v1, f_coords[1], rtol=1e-6, atol=1e-8) and
np.allclose(v2, f_coords[0], rtol=1e-6, atol=1e-8)):
local_edge = e
break
# 3) distance-based fallback (most robust)
if local_edge == -1:
min_dist = np.inf
best_edge = -1
for e in range(4):
v1 = cell_coords[e]
v2 = cell_coords[(e+1)%4]
d1 = np.linalg.norm(f_coords[0] - v1) + np.linalg.norm(f_coords[1] - v2)
d2 = np.linalg.norm(f_coords[0] - v2) + np.linalg.norm(f_coords[1] - v1)
d = min(d1, d2)
if d < min_dist:
min_dist = d
best_edge = e
if best_edge != -1:
local_edge = best_edge
if local_edge == -1:
raise ValueError("Facet not found in cell")
t_pts = q_pts_facet[:, 0]
xi_eta = np.array([map_interval_to_quad_edge(t, local_edge) for t in t_pts])
tab_T_trace = elem_T.tabulate(1, xi_eta)
tab_phi_trace = elem_phi.tabulate(1, xi_eta)
phi_T_trace = tab_T_trace[0]
phi_phi_trace = tab_phi_trace[0]
grad_phi_ref_trace = tab_phi_trace[1]
if grad_phi_ref_trace.ndim == 3:
if grad_phi_ref_trace.shape[2] == 1:
grad_phi_ref_trace = np.concatenate([grad_phi_ref_trace, grad_phi_ref_trace], axis=2)
elif grad_phi_ref_trace.shape[2] != 2:
grad_phi_ref_trace = grad_phi_ref_trace[:, :, :2]
normal, _ = facet_normal_and_length(mesh_, f)
w_f = q_wts_facet * phys_length
f_mesh_cell = facet_map[f]
dofs_F = map_F.cell_dofs(f_mesh_cell)
nF = len(dofs_F)
tab_F = elem_F.tabulate(0, q_pts_facet)
phi_F_vals = tab_F[0]
M_TT_edge = np.zeros((nT, nT))
M_TF_edge = np.zeros((nT, nF))
M_FF_edge = np.zeros((nF, nF))
RHS_T_edge = np.zeros((nphi, nT))
RHS_F_edge = np.zeros((nphi, nF))
for q in range(len(t_pts)):
xi, eta = xi_eta[q]
w = w_f[q]
phi_Tq = phi_T_trace[q, :]
phi_Fq = phi_F_vals[q, :]
N = np.array([(1-xi)*(1-eta), xi*(1-eta), xi*eta, (1-xi)*eta])
dN_dxi = np.array([-(1-eta), (1-eta), eta, -eta])
dN_deta = np.array([-(1-xi), -xi, xi, (1-xi)])
J = np.zeros((2,2))
J[0,0] = np.dot(dN_dxi, cell_coords[:,0])
J[0,1] = np.dot(dN_deta, cell_coords[:,0])
J[1,0] = np.dot(dN_dxi, cell_coords[:,1])
J[1,1] = np.dot(dN_deta, cell_coords[:,1])
# ---- Safe inverse ----
invJ = safe_inv_2x2(J)
grad_phi_q = grad_phi_ref_trace[q, :, :] @ invJ.T
n_dot_grad = np.dot(grad_phi_q, normal)
for i in range(nphi):
for j in range(nT):
RHS_T_edge[i, j] += w * (-phi_Tq[j] * n_dot_grad[i])
for j in range(nF):
RHS_F_edge[i, j] += w * (phi_Fq[j] * n_dot_grad[i])
inv_h = 1.0 / h
for i in range(nT):
for j in range(nT):
M_TT_edge[i, j] += w * inv_h * phi_Tq[i] * phi_Tq[j]
for j in range(nF):
M_TF_edge[i, j] += w * inv_h * phi_Tq[i] * phi_Fq[j]
for i in range(nF):
for j in range(nF):
M_FF_edge[i, j] += w * inv_h * phi_Fq[i] * phi_Fq[j]
facet_contribs.append((dofs_F, M_TT_edge, M_TF_edge, M_FF_edge, RHS_T_edge, RHS_F_edge))
RHS_T_loc += RHS_T_edge
# Build RHS_F_loc preserving order
all_F_dofs = []
for (dofs_F, _, _, _, _, _) in facet_contribs:
for dof in dofs_F:
if dof not in all_F_dofs:
all_F_dofs.append(dof)
total_nF_cell = len(all_F_dofs)
if total_nF_cell > 0:
RHS_F_loc = np.zeros((nphi, total_nF_cell))
for (dofs_F, _, _, _, _, RHS_F_edge) in facet_contribs:
cols = [all_F_dofs.index(dof) for dof in dofs_F]
for i in range(nphi):
for j, col in enumerate(cols):
RHS_F_loc[i, col] += RHS_F_edge[i, j]
RHS = np.hstack([RHS_T_loc, RHS_F_loc])
else:
RHS = RHS_T_loc
# Assemble S_phi_loc
for q in range(len(q_pts_cell)):
xi, eta = q_pts_cell[q]
wq = q_wts_cell[q]
N = np.array([(1-xi)*(1-eta), xi*(1-eta), xi*eta, (1-xi)*eta])
dN_dxi = np.array([-(1-eta), (1-eta), eta, -eta])
dN_deta = np.array([-(1-xi), -xi, xi, (1-xi)])
J = np.zeros((2,2))
J[0,0] = np.dot(dN_dxi, cell_coords[:,0])
J[0,1] = np.dot(dN_deta, cell_coords[:,0])
J[1,0] = np.dot(dN_dxi, cell_coords[:,1])
J[1,1] = np.dot(dN_deta, cell_coords[:,1])
detJ = np.linalg.det(J)
w = wq * abs(detJ)
invJ = safe_inv_2x2(J)
grad_phi_q = grad_phi_ref[q, :, :] @ invJ.T
for i in range(nphi):
for j in range(nphi):
S_phi_loc[i, j] += w * np.dot(grad_phi_q[i, :], grad_phi_q[j, :])
# Average constraint
int_phi = np.zeros(nphi)
for q in range(len(q_pts_cell)):
xi, eta = q_pts_cell[q]
wq = q_wts_cell[q]
N = np.array([(1-xi)*(1-eta), xi*(1-eta), xi*eta, (1-xi)*eta])
dN_dxi = np.array([-(1-eta), (1-eta), eta, -eta])
dN_deta = np.array([-(1-xi), -xi, xi, (1-xi)])
J = np.zeros((2,2))
J[0,0] = np.dot(dN_dxi, cell_coords[:,0])
J[0,1] = np.dot(dN_deta, cell_coords[:,0])
J[1,0] = np.dot(dN_dxi, cell_coords[:,1])
J[1,1] = np.dot(dN_deta, cell_coords[:,1])
detJ = np.linalg.det(J)
int_phi += wq * abs(detJ) * phi_phi_vals[q, :]
n_rhs = RHS.shape[1]
Aug = np.zeros((nphi+1, nphi+1))
Aug[:nphi, :nphi] = S_phi_loc
Aug[:nphi, nphi] = int_phi
Aug[nphi, :nphi] = int_phi
RHS_aug = np.vstack([RHS, np.zeros((1, n_rhs))])
try:
sol_aug = np.linalg.solve(Aug, RHS_aug)
except np.linalg.LinAlgError:
sol_aug = np.zeros((nphi+1, n_rhs))
sol_phi = sol_aug[:nphi, :]
potT_loc = sol_phi[:, :nT]
potF_loc = sol_phi[:, nT:] if total_nF_cell > 0 else np.zeros((nphi, 0))
for i, gi in enumerate(dofs_phi):
for j, gj in enumerate(dofs_T):
potT[gi, gj] = potT_loc[i, j]
for i, gi in enumerate(dofs_phi):
for j, dof_F in enumerate(all_F_dofs):
potF[gi, dof_F] = potF_loc[i, j]
for (dofs_F, M_TT_edge, M_TF_edge, M_FF_edge, _, _) in facet_contribs:
for i, gi in enumerate(dofs_T):
for j, gj in enumerate(dofs_T):
M_TT_facet[gi, gj] += M_TT_edge[i, j]
for j, gj in enumerate(dofs_F):
M_TF_facet[gi, gj] += M_TF_edge[i, j]
for i, gi in enumerate(dofs_F):
for j, gj in enumerate(dofs_F):
M_FF_facet[gi, gj] += M_FF_edge[i, j]
# ------------------------------------------------------------------
# 5. projFd (facet projection)
# ------------------------------------------------------------------
for f in range(num_facets_global):
facet_vert_indices = facet_vertices.links(f)
f_coords = coords[facet_vert_indices][:, :2]
phys_length = np.linalg.norm(f_coords[1] - f_coords[0])
f_mesh_cell = facet_map[f]
dofs_F = map_F.cell_dofs(f_mesh_cell)
nF = len(dofs_F)
adj = f_to_c.links(f)
side_cell = adj[0]
dofs_phi_side = map_phi.cell_dofs(side_cell)
nphi_side = len(dofs_phi_side)
M_FF_loc = np.zeros((nF, nF))
M_F_phi_loc = np.zeros((nF, nphi_side))
tab_F = elem_F.tabulate(0, q_pts_facet)
phi_F_vals = tab_F[0]
cell_vert_indices = vertices_cell.links(side_cell)
cell_coords_side = coords[cell_vert_indices][:, :2]
# ---- Determine local edge ----
local_edge = -1
# index-based
for e in range(4):
v1_idx = cell_vert_indices[e]
v2_idx = cell_vert_indices[(e+1)%4]
if (v1_idx == facet_vert_indices[0] and v2_idx == facet_vert_indices[1]) or \
(v1_idx == facet_vert_indices[1] and v2_idx == facet_vert_indices[0]):
local_edge = e
break
# coordinate-based fallback
if local_edge == -1:
for e in range(4):
v1 = cell_coords_side[e]
v2 = cell_coords_side[(e+1)%4]
if (np.allclose(v1, f_coords[0], rtol=1e-6, atol=1e-8) and
np.allclose(v2, f_coords[1], rtol=1e-6, atol=1e-8)) or \
(np.allclose(v1, f_coords[1], rtol=1e-6, atol=1e-8) and
np.allclose(v2, f_coords[0], rtol=1e-6, atol=1e-8)):
local_edge = e
break
# distance-based fallback
if local_edge == -1:
min_dist = np.inf
best_edge = -1
for e in range(4):
v1 = cell_coords_side[e]
v2 = cell_coords_side[(e+1)%4]
d1 = np.linalg.norm(f_coords[0] - v1) + np.linalg.norm(f_coords[1] - v2)
d2 = np.linalg.norm(f_coords[0] - v2) + np.linalg.norm(f_coords[1] - v1)
d = min(d1, d2)
if d < min_dist:
min_dist = d
best_edge = e
if best_edge != -1:
local_edge = best_edge
if local_edge == -1:
raise ValueError("Facet not found in side cell")
t_pts = q_pts_facet[:, 0]
xi_eta = np.array([map_interval_to_quad_edge(t, local_edge) for t in t_pts])
tab_phi_trace = elem_phi.tabulate(0, xi_eta)
phi_phi_trace = tab_phi_trace[0]
w_f = q_wts_facet * phys_length
for q in range(len(t_pts)):
w = w_f[q]
phi_F = phi_F_vals[q, :]
phi_phi = phi_phi_trace[q, :]
for i in range(nF):
for j in range(nF):
M_FF_loc[i, j] += w * phi_F[i] * phi_F[j]
for j in range(nphi_side):
M_F_phi_loc[i, j] += w * phi_F[i] * phi_phi[j]
try:
projFd_loc = np.linalg.solve(M_FF_loc, M_F_phi_loc)
except np.linalg.LinAlgError:
projFd_loc = np.zeros((nF, nphi_side))
for i, gi in enumerate(dofs_F):
for j, gj in enumerate(dofs_phi_side):
projFd[gi, gj] = projFd_loc[i, j]
# ------------------------------------------------------------------
# 6. Build reduced system
# ------------------------------------------------------------------
projT = projT.tocsc()
projFd = projFd.tocsc()
potT = potT.tocsc()
potF = potF.tocsc()
M_TT_facet = M_TT_facet.tocsc()
M_TF_facet = M_TF_facet.tocsc()
M_FF_facet = M_FF_facet.tocsc()
I_T = sp.eye(N_T, format='csc')
I_F = sp.eye(N_F, format='csc')
CTT = projT * potT - I_T
CTF = projT * potF
CFT = projFd * potT
CFF = projFd * potF - I_F
C = sp.bmat([[CTT, CTF],
[CFT, CFF]], format='csc')
M_facet = sp.bmat([[M_TT_facet, M_TF_facet],
[M_TF_facet.T, M_FF_facet]], format='csc')
K_facet = C.T * M_facet * C
A_TT = potT.T * S_phi * potT
A_TF = potT.T * S_phi * potF
A_FF = potF.T * S_phi * potF
A_top = sp.bmat([[A_TT, A_TF],
[A_TF.T, A_FF]], format='csc')
A = A_top + K_facet
# ------------------------------------------------------------------
# 7. Right-hand side and BCs
# ------------------------------------------------------------------
uT_trial = ufl.TrialFunction(V_T)
vT_test = ufl.TestFunction(V_T)
x = ufl.SpatialCoordinate(mesh_)
f = 2.0 * ufl.pi**2 * ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1])
bT_form = ufl.inner(f, vT_test) * ufl.dx
bT_vec = fem.assemble_vector(fem.form(bT_form))
bT = bT_vec.array[:N_T]
bF = np.zeros(N_F)
b = np.concatenate([bT, bF])
def on_boundary(x):
return np.logical_or.reduce([np.isclose(x[0], 0.0),
np.isclose(x[0], 1.0),
np.isclose(x[1], 0.0),
np.isclose(x[1], 1.0)])
bndry_facets = mesh.locate_entities_boundary(mesh_, mesh_.topology.dim-1, on_boundary)
bndry_F_dofs = []
for f in bndry_facets:
f_mesh_cell = facet_map[f]
dofs = map_F.cell_dofs(f_mesh_cell)
bndry_F_dofs.extend(dofs)
bndry_F_dofs = np.unique(bndry_F_dofs)
A = A.tolil()
for i in bndry_F_dofs:
idx = N_T + i
A[idx, :] = 0
A[:, idx] = 0
A[idx, idx] = 1.0
b[idx] = 0.0
A = A.tocsc()
# ------------------------------------------------------------------
# 8. Solve and recover fields
# ------------------------------------------------------------------
reg = 1e-12
A_reg = A + reg * sp.eye(A.shape[0], format='csc')
sol = spsolve(A_reg, b)
print(f" sol norm: {np.linalg.norm(sol):.2e}, residual norm: {np.linalg.norm(A_reg*sol - b):.2e}")
uT_vec = sol[:N_T]
uF_vec = sol[N_T:]
phiu_vec = potT * uT_vec + potF * uF_vec
uT_func = fem.Function(V_T)
uT_func.x.array[:] = uT_vec
uF_func = fem.Function(V_F)
uF_func.x.array[:] = uF_vec
phiu_func = fem.Function(V_phi)
phiu_func.x.array[:] = phiu_vec
# ------------------------------------------------------------------
# 9. Errors
# ------------------------------------------------------------------
u_exact_func = fem.Function(V_phi)
u_exact_func.interpolate(u_exact)
L2 = fem.assemble_scalar(fem.form((phiu_func - u_exact_func)**2 * ufl.dx))
L2error[n] = np.sqrt(L2)
grad_phi = ufl.grad(phiu_func)
grad_exact = ufl.grad(u_exact_func)
energy_cell = (grad_phi - grad_exact)**2 * ufl.dx
DG_cell = np.sqrt(fem.assemble_scalar(fem.form(energy_cell)))
dS = ufl.Measure("dS", domain=mesh_)
jump = (phiu_func - u_exact_func)('+') - (phiu_func - u_exact_func)('-')
jump_form = (jump**2) / h * dS
DG_jump = np.sqrt(fem.assemble_scalar(fem.form(jump_form)))
DGerror[n] = np.sqrt(DG_cell**2 + DG_jump**2)
grad_uT = ufl.grad(uT_func)
grad_u_exact = ufl.grad(u_exact_func)
RG_cell = (grad_uT - grad_u_exact)**2 * ufl.dx
RG_cell_val = np.sqrt(fem.assemble_scalar(fem.form(RG_cell)))
RGerror[n] = RG_cell_val
print(f" L2error = {L2error[n]:.5e}, DGerror = {DGerror[n]:.5e}, RGerror = {RGerror[n]:.5e}")
# ------------------------------------------------------------------
# 10. Convergence rates and plot
# ------------------------------------------------------------------
print("\nConvergence rates:")
for n in range(1, n_ref):
rate_L2 = np.log(L2error[n-1]/L2error[n]) / np.log(Dx[n-1]/Dx[n])
rate_DG = np.log(DGerror[n-1]/DGerror[n]) / np.log(Dx[n-1]/Dx[n])
rate_RG = np.log(RGerror[n-1]/RGerror[n]) / np.log(Dx[n-1]/Dx[n])
print(f"Level {n}: L2 rate = {rate_L2:.3f}, DG rate = {rate_DG:.3f}, RG rate = {rate_RG:.3f}")
plt.figure()
plt.loglog(Dx, L2error, 'o-', label='L2')
plt.loglog(Dx, DGerror, 's-', label='DG')
plt.loglog(Dx, RGerror, 'd-', label='RG')
plt.xlabel('h')
plt.ylabel('Error')
plt.legend()
plt.grid(True)
plt.show()
Please help!. I would greatly appreciate any help or guidance on how to achieve superconvergence.
Thanks in advance!.