Define a boundary from DG elements


I have a 2D rectangular geometry. The boundaries are the inlet, the outlet and the wall. I want to impose a Dirichlet BC for an unknown using DG elements. By using the classical DirichletBC() function I can only impose a “pointwise” BC on the centroid of an element. To this end I would like to define a boundary that contains the centroids of the elements, lets say at the inlet, in order to use the DirichletBC() function for a boundary than for a point. Initially, I exported the centroids of the mesh and then I made a list with those that belong to the inlet boundary, the inletCentroids list. Finally, I do not get any error but the inflow sub_boundary does not contain any centroid. The centroids are calculated correctly, so I guess that the problem is inside class Inflow(). Here is the sample of my code

# Find centroids of the mesh

centroid = []

for cell_no in range(mesh.num_cells()):
	coc = Cell(mesh, cell_no).midpoint()
	coo = []

# Store the centroids of the inlet

inletCentroids  = []
for i in range(len(centroid)):
	if (abs(centroid[i][0]-0.0041)<1.e-4):

# Sub domain for inlet 
class Inflow(SubDomain):
	def inside(self, x, on_boundary):
		tol = 1.e-4
            for i in inletCentroids:
		      p = near( x[0], centroid[i][0] , tol) and near(x[1],  centroid[i][1]  , tol) 
		      if ( p is True):
			return (p)

sub_domains = MeshFunction("size_t", mesh,0)
inflow = Inflow()
inflow.mark(sub_domains, 1)
inflowBC = DirichletBC(VP.sub(1), pIn , sub_boundaries , 1 ) 

Thank you

It is not common to use strong Dirichlet conditions with DG elements (and I do not see any simple way of modifying the DirichletBC interface to accommodate for your use-case).

With dolfin I would suggest to weakly impose the boundary condition using Nitsche’s method.

You can also try to use dolfinx, where setting such a condition is quite straigthforward as shown below:

import dolfinx
from mpi4py import MPI
import ufl
import numpy as np

mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 5, 5)
vertex_to_cell = mesh.topology.connectivity(0, mesh.topology.dim)

def left_boundary(x):
    return np.isclose(x[0], 0)

boundary_vertices = dolfinx.mesh.locate_entities(mesh, 0, left_boundary)
boundary_cells = []
for vertex in boundary_vertices:
boundary_cells = np.hstack(boundary_cells)

el = ufl.MixedElement([ufl.VectorElement(
    "CG", mesh.ufl_cell(), 2), ufl.FiniteElement("DG", mesh.ufl_cell(), 0)])
VP = dolfinx.FunctionSpace(mesh, el)
P = VP.sub(1).collapse()
p = dolfinx.Function(P)
p.interpolate(lambda x: np.cos(np.pi*x[1]))
boundary_dofs = dolfinx.fem.locate_dofs_topological(
    (VP.sub(1), P), mesh.topology.dim, boundary_cells)

inflowBC = dolfinx.DirichletBC(p, boundary_dofs, VP.sub(1))

p_test = dolfinx.Function(VP)
dolfinx.fem.set_bc(p_test.vector, [inflowBC])

with, "p_bc.xdmf", "w") as xdmf:
1 Like