Hello, community!
I have a question that may sound very simple.
I want to solve problems with the Helmholtz equation in the acoustic domain and extra equations with extra unknowns only defined on the boundary of the acoustic domain. Those extra unknowns appear in the boundary integral of the main Helmholtz equation in the weak form. And the main unknown, the acoustic pressure, appears on the boundary extra equation.
The idea is to create a sub_mesh which is the boundary of the acoustic domain where the extra equation is defined. Then define a trial function on each mesh and a test function on each mesh. Define the differents forms, and assemble the whole thing in a matrix.
Here is the code. It’s a bit long cause there is all the importation of the mesh from gmsh and the creation of the sub mesh. But the relevant lines are the 12th last ones. I’ve only put the problematic matrix C because the other ones are assembled without error.
import numpy as np
import meshio
import dolfinx.mesh as msh
from dolfinx.io import XDMFFile
from mpi4py import MPI
from dolfinx.fem import Function, FunctionSpace, assemble, form
from ufl import (TestFunction, FiniteElement, TrialFunction, dx, grad, inner, Measure)
# Create a meshio mesh including physical markers for the given type.
def create_mesh(mesh, cell_type):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points
# Create output mesh
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells},
cell_data={"name_to_read": [cell_data]})
return out_mesh
# Read in mesh from Gmsh, create mesh(volume and surface elements),
###################################################################
mesh3D_from_msh_lin = meshio.read("cubic_truncated_domain_small.msh")
mesh_tetra_lin = create_mesh(mesh3D_from_msh_lin, "tetra")
rank = MPI.COMM_WORLD.Get_rank()
if rank == 0:
meshio.write("cubic_truncated_domain_small_msh_tetra_lin.xdmf", mesh_tetra_lin)
# Read mesh from xdmf files: first volume elements, then boundary triangles
with XDMFFile(MPI.COMM_WORLD, "cubic_truncated_domain_small_msh_tetra_lin.xdmf", "r") as xdmf:
mesh_lin = xdmf.read_mesh(name="Grid")
mesh_lin_tags = xdmf.read_meshtags(mesh_lin, name="Grid")
tdim = mesh_lin.topology.dim
fdim = tdim - 1
mesh_lin.topology.create_connectivity(tdim, fdim)
deg = 1
mesh = mesh_lin
mesh_tags = mesh_lin_tags
mesh_bc = mesh_bc_lin
mesh_bc_tags = mesh_bc_lin_tags
#Fonctions to localize the aimed entities
def gamma1_3_facets(x):
return np.isclose(x[0], 0)
#Find the entities which will be in the submesh
gamma1_3_facets_list = msh.locate_entities_boundary(mesh, fdim, gamma1_3_facets)
submesh, entity_map, _, _ = msh.create_submesh(mesh, fdim, gamma1_3_facets_list)
P_el = FiniteElement("Lagrange", mesh.ufl_cell(), deg)
Q_el = FiniteElement("Lagrange", submesh.ufl_cell(), deg)
P = FunctionSpace(mesh, P_el)
Q = FunctionSpace(submesh, Q_el)
p, q = TrialFunction(P), TrialFunction(Q)
v, w = TestFunction(P), TestFunction(Q)
ds = Measure("ds", domain=mesh, subdomain_data=mesh_bc_tags)
k = inner(grad(p), grad(v)) * dx
m = inner(p, v) * dx
c = inner(q, v)*ds(3)
C = petsc.assemble_matrix(form(c))
C.assemble()
I already know that the problem is the trial function q and the test function v don’t have the same size, then the inner product can’t work. But do you have any idea?
Thank you !
And see you in Cagliari!