search

firedrakeproject / firedrake

Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM)
https://firedrakeproject.org
Other
498 stars 157 forks source link

BUG: DirichletBC is not correctly applied for facet element of polynomial degree greater than 1 over quadrilateral meshes #3077

Open a-brugnoli opened 1 year ago

a-brugnoli commented 1 year ago

The imposition of inhomogeneous essential boundary condition on facet spaces over quadrilateral does not work when degree higher than 1. This is probably related to the fact that projection is not supported for facet elements. The following code solves the Poisson equation using the hybrid formulation introduced in "Locally conservative fluxes for the Continuous Galerkin method". The result clearly shows a problem with the imposition of the inhomogeneous boundary condition:

import firedrake as fdrk
import matplotlib.pyplot as plt

n_elem = 4
deg = 2
quad = True
mesh = fdrk.UnitSquareMesh(n_elem, n_elem, quadrilateral=quad)
normal_versor = fdrk.FacetNormal(mesh)
x, y = fdrk.SpatialCoordinate(mesh)

cell = mesh.ufl_cell()
CG_element = fdrk.FiniteElement("CG", cell, deg) 
broken_CG_element = fdrk.BrokenElement(CG_element)
facet_CG_element = CG_element[fdrk.facet]
brokenfacet_CG_element = fdrk.BrokenElement(facet_CG_element)

broken_CG_space = fdrk.FunctionSpace(mesh, broken_CG_element)
brokenfacet_CG_space = fdrk.FunctionSpace(mesh, brokenfacet_CG_element)
facet_CG_space = fdrk.FunctionSpace(mesh, facet_CG_element)

mixed_space = broken_CG_space * brokenfacet_CG_space * facet_CG_space

test_pressure, test_normaltrace, test_tangentialtrace = fdrk.TestFunctions(mixed_space)
trial_pressure, trial_normaltrace, trial_tangentialtrace = fdrk.TrialFunctions(mixed_space)

control_local = fdrk.inner(test_pressure, trial_normaltrace)
control_local_adj = fdrk.inner(test_normaltrace, trial_pressure)

control_global = fdrk.inner(test_normaltrace, trial_tangentialtrace)
control_global_adj = fdrk.inner(test_tangentialtrace, trial_normaltrace)

constr_local = + (control_local('+') + control_local('-')) * fdrk.dS + control_local * fdrk.ds \
                - ((control_local_adj('+') + control_local_adj('-')) * fdrk.dS + control_local_adj * fdrk.ds)

constr_global = + (control_global('+') + control_global('-')) * fdrk.dS + control_global * fdrk.ds \
                - ((control_global_adj('+') + control_global_adj('-')) * fdrk.dS + control_global_adj * fdrk.ds)

A_operator = fdrk.inner(fdrk.grad(test_pressure), fdrk.grad(trial_pressure)) * fdrk.dx \
                    - constr_local - constr_global

exact_solution = fdrk.sin(x)*fdrk.sin(y)
forcing = -fdrk.div(fdrk.grad(exact_solution))

b_functional = test_pressure*forcing*fdrk.dx 

solution = fdrk.Function(mixed_space)
bc_dirichlet = fdrk.DirichletBC(mixed_space.sub(2), exact_solution, "on_boundary")
problem = fdrk.LinearVariationalProblem(A_operator, b_functional, solution, bcs=bc_dirichlet)
solver =  fdrk.LinearVariationalSolver(problem)
solver.solve()

fdrk.trisurf(solution.sub(0))

plt.show()
dham commented 1 year ago

It appears you are indeed correct. I think I have successfully worked around this by explicitly projecting the BC into the right space:

import firedrake as fdrk
import matplotlib.pyplot as plt

n_elem = 4
deg = 2
quad = True
mesh = fdrk.UnitSquareMesh(n_elem, n_elem, quadrilateral=quad)
normal_versor = fdrk.FacetNormal(mesh)
x, y = fdrk.SpatialCoordinate(mesh)

cell = mesh.ufl_cell()
CG_element = fdrk.FiniteElement("CG", cell, deg) 
broken_CG_element = fdrk.BrokenElement(CG_element)
facet_CG_element = CG_element[fdrk.facet]
brokenfacet_CG_element = fdrk.BrokenElement(facet_CG_element)

broken_CG_space = fdrk.FunctionSpace(mesh, broken_CG_element)
brokenfacet_CG_space = fdrk.FunctionSpace(mesh, brokenfacet_CG_element)
facet_CG_space = fdrk.FunctionSpace(mesh, facet_CG_element)

mixed_space = broken_CG_space * brokenfacet_CG_space * facet_CG_space

test_pressure, test_normaltrace, test_tangentialtrace = fdrk.TestFunctions(mixed_space)
trial_pressure, trial_normaltrace, trial_tangentialtrace = fdrk.TrialFunctions(mixed_space)

control_local = fdrk.inner(test_pressure, trial_normaltrace)
control_local_adj = fdrk.inner(test_normaltrace, trial_pressure)

control_global = fdrk.inner(test_normaltrace, trial_tangentialtrace)
control_global_adj = fdrk.inner(test_tangentialtrace, trial_normaltrace)

constr_local = + (control_local('+') + control_local('-')) * fdrk.dS + control_local * fdrk.ds \
                - ((control_local_adj('+') + control_local_adj('-')) * fdrk.dS + control_local_adj * fdrk.ds)

constr_global = + (control_global('+') + control_global('-')) * fdrk.dS + control_global * fdrk.ds \
                - ((control_global_adj('+') + control_global_adj('-')) * fdrk.dS + control_global_adj * fdrk.ds)

A_operator = fdrk.inner(fdrk.grad(test_pressure), fdrk.grad(trial_pressure)) * fdrk.dx \
                    - constr_local - constr_global

exact_solution = fdrk.sin(x)*fdrk.sin(y)
forcing = -fdrk.div(fdrk.grad(exact_solution))

b_functional = test_pressure*forcing*fdrk.dx 

bc_value = fdrk.Function(facet_CG_space)
bc_test = fdrk.TestFunction(facet_CG_space)
bc_trial = fdrk.TrialFunction(facet_CG_space)
fdrk.solve(bc_test*bc_trial*fdrk.ds + bc_test("+")*bc_trial("+")*fdrk.dS
           == bc_test*exact_solution*fdrk.ds
           + bc_test("+")*exact_solution("+")*fdrk.dS, bc_value)

solution = fdrk.Function(mixed_space)
bc_dirichlet = fdrk.DirichletBC(mixed_space.sub(2), bc_value, "on_boundary")
problem = fdrk.LinearVariationalProblem(A_operator, b_functional, solution, bcs=bc_dirichlet)
solver = fdrk.LinearVariationalSolver(problem)
solver.solve()

fdrk.trisurf(solution.sub(0))

plt.show()
a-brugnoli commented 1 year ago

Thank you so much for the help. The workaround works. A little discrepancy between the projected field and the corresponding facet projected field still remains:

import firedrake as fdrk

n_elem = 4
deg = 2
quad = True

mesh = fdrk.UnitSquareMesh(n_elem, n_elem, quadrilateral=quad)
cell_diameter = fdrk.CellDiameter(mesh)
normal_versor = fdrk.FacetNormal(mesh)
x, y = fdrk.SpatialCoordinate(mesh)
cell = mesh.ufl_cell()

CG_element = fdrk.FiniteElement("CG", cell, deg) 
facet_CG_element = CG_element[fdrk.facet]

CG_space = fdrk.FunctionSpace(mesh, CG_element)
facet_CG_space = fdrk.FunctionSpace(mesh, facet_CG_element)

exact_solution = fdrk.sin(x)*fdrk.sin(y)

projected_facet_exact = fdrk.Function(facet_CG_space)
facet_test = fdrk.TestFunction(facet_CG_space)
facet_trial = fdrk.TrialFunction(facet_CG_space)
fdrk.solve(facet_test*facet_trial*fdrk.ds + facet_test("+")*facet_trial("+")*fdrk.dS
           == facet_test*exact_solution*fdrk.ds
           + facet_test("+")*exact_solution("+")*fdrk.dS, projected_facet_exact)

projected_exact = fdrk.project(exact_solution, CG_space)

boundary_integrand = cell_diameter * (projected_exact - projected_facet_exact) ** 2

square_norm = boundary_integrand('+') * fdrk.dS + boundary_integrand * fdrk.ds

print(f"Error: {fdrk.sqrt(fdrk.assemble(square_norm))}")