Him I am currently working on a project involving the variational formulation of Maxwell’s equations using Discontinuous Galerkin (DG) finite elements with a particular emphasis on incorporating the cross operator. I am reaching out to you as I have encountered some challenges in my implementation.
here:
the {.} is the avarange operator, [.]_N is the normal jump and [.]_T is the tangant jump operator. since in Error for more or less standard variational formulation using two cross products the cross operator is only defined for vectors of length three! i can’t define this operator in the variational form ! and i got the this error:
(A, b) = assemble_system(a, L, bc)
................................................
return Product(a[i], b[j]) - Product(a[j], b[i])
File "/usr/lib/python3/dist-packages/ufl_legacy/exproperators.py", line 438, in _getitem
all_indices, slice_indices, repeated_indices = create_slice_indices(component, shape, self.ufl_free_indices)
File "/usr/lib/python3/dist-packages/ufl_legacy/index_combination_utils.py", line 152, in create_slice_indices
error("Index out of bounds.")
File "/usr/lib/python3/dist-packages/ufl_legacy/log.py", line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl_legacy.log.UFLException: Index out of bounds.
Here is a minimal code:
import matplotlib.pyplot as plt
from dolfin import *
import numpy as np
import time
def solve_pde(n_mesh, bc, FE_CHOICE, Plot):
t0 = time.time()
nx = 1
ny = 1
deg = 2
mesh = UnitSquareMesh(n_mesh, n_mesh) # ,"crossed")
# Define function space
CHOICE_FEM = FE_CHOICE
P2 = VectorElement(CHOICE_FEM, triangle, deg)
P1 = FiniteElement(CHOICE_FEM, triangle, deg - 1)
TH = MixedElement([P2, P1])
Mix_h = FunctionSpace(mesh, TH)
f = Expression(("-4*x[1]*exp(x[1])-(pow(x[1],2)-x[1])*exp(x[0]*x[1])*(x[0]*x[1]*(x[0]-1)+2*x[0]-1)",
"-4*x[0]*exp(x[0])-(pow(x[0],2)-x[0])*exp(x[0]*x[1])*(x[0]*x[1]*(x[1]-1)+2*x[1]-1)"),
element=VectorElement("Lagrange", triangle, 6))
f = interpolate(f, VectorFunctionSpace(mesh, "Lagrange", degree=3))
u_exact = Expression(("x[1]*(x[1]-1)*exp(x[1])", "x[0]*(x[0]-1)*exp(x[0])"),
element=VectorElement("Lagrange", triangle, 6))
p_exact = Expression("x[0]*x[1]*(x[0]-1)*(x[1]-1)*exp(x[0]*x[1])",
element=FiniteElement("Lagrange", triangle, 5))
w0 = Constant("0.0")
def project_boundary(x, on_boundary):
tol = 1.0e-14
return on_boundary
def u0_boundary(x, on_boundary):
tol = 1.0e-14
return on_boundary and (abs(x[1] - 1) < tol or abs(x[1]) < tol or (abs(x[1]) < tol and x[0] > tol))
def u1_boundary(x, on_boundary):
tol = 1.0e-14
return on_boundary and (abs(x[0] - 1) < tol or abs(x[0]) < tol or (abs(x[0]) < tol and x[1] < tol))
w0 = Constant((0.0))
bc0 = DirichletBC(Mix_h.sub(0).sub(0), w0, u0_boundary)
bc1 = DirichletBC(Mix_h.sub(0).sub(1), w0, u1_boundary)
bc_uxn = [bc0, bc1]
bc_proj_u = DirichletBC(Mix_h.sub(0), Constant((0, 0)), project_boundary)
bc_proj_p = DirichletBC(Mix_h.sub(1), p_exact, project_boundary)
bc_projection = [bc_proj_u, bc_proj_p]
t1 = time.time()
def boundary(x, on_boundary):
return on_boundary
bcu_dex = DirichletBC(Mix_h.sub(0), u_exact, boundary)
bc_pex = DirichletBC(Mix_h.sub(1), p_exact, boundary)
bc_ex = [bcu_dex, bc_pex]
if bc == "uxn":
bc = bc_uxn
elif bc == "projection":
bc = bc_projection
elif bc == "exact":
bc = bc_ex
U = TrialFunction(Mix_h)
V = TestFunction(Mix_h)
(u_h, p) = split(U)
(v_h, q) = split(V)
h = mesh.hmax()
n = FacetNormal(mesh)
n = FacetNormal(mesh)
omega = Constant('1.0')
lamda = 1
eta_a = 1
r = 1
sigma = 1
mu = 1
# Define bilinear form
a_form = (omega ** 2 * inner(u_h, v_h) * dx + inner(curl(u_h), curl(v_h)) * dx + inner(div(u_h), div(v_h)) * dx
+ eta_a * inner(jump(u_h, n), jump(v_h, n)) * dS
- inner(jump(u_h, n), avg(curl(v_h))) * dS # +-
+ inner(jump(v_h, n), avg(curl(u_h))) * dS # +-
+inner(cross(u_h, n), cross(v_h, n))*dS
)
b_form = (inner(p, div(v_h)) * dx - eta_a * inner(avg(p), jump(v_h, n)) * dS # integrale frme B
- inner(div(u_h), q) * dx + eta_a * inner(avg(q), jump(u_h, n)) * dS # integrale frme B
)
L = inner(f, v_h) * dx
a = a_form + b_form
(A, b) = assemble_system(a, L, bc)
solver = KrylovSolver("cg", "ilu")
solver.parameters["absolute_tolerance"] = 1E-14
solver.parameters["relative_tolerance"] = 1E-14
solver.parameters["maximum_iterations"] = 1000
U_0 = Function(Mix_h)
solve(A, U_0.vector(), b)
(u0, p0) = U_0.split(deepcopy=True)
return u0, p0
bc = "exact" # or "projection" or "exact"
FE_CHOICE = "DG" # or "DG" or other choices
n_mesh = 32
solve_pde(n_mesh, bc, FE_CHOICE, True)
Warm regards,