Define a boundary from DG elements

1

Hello,

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 = []
	coo.append(coc.x())
	coo.append(coc.y())
	centroid.append(coo)

# Store the centroids of the inlet

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


# 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)
sub_domains.set_all(0)
inflow = Inflow()
inflow.mark(sub_domains, 1)
inflowBC = DirichletBC(VP.sub(1), pIn , sub_boundaries , 1 )

Thank you

2

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
import dolfinx.io
from mpi4py import MPI
import ufl
import numpy as np

mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 5, 5)
mesh.topology.create_connectivity_all()
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.append(vertex_to_cell.links(vertex))
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 dolfinx.io.XDMFFile(MPI.COMM_WORLD, "p_bc.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(p_test.sub(1))
1 Like
3

Dear @dokken, I tried to run the code you provided. I’m running FENICSX Ubuntu’s PPA.

First, I got the error

File “subdomains.py”, line 8, in
mesh.topology.create_connectivity_all()
AttributeError: ‘dolfinx.cpp.mesh.Topology’ object has no attribute ‘create_connectivity_all’

Changing mesh.topology.create_connectivity_all() to mesh.topology.create_connectivity(1,2) gives a new error:

Traceback (most recent call last):
File “subdomains.py”, line 19, in
boundary_cells.append(vertex_to_cell.links(vertex))
AttributeError: ‘NoneType’ object has no attribute ‘links’

Am I missing something? . Thanks!

4

You need to change it to:

mesh.topology.create_connectivity(0, mesh.topology.dim)

as this is the connectivity that you request in the code.

1 Like
5

This work perfectly. On the other hand, I’m having a strange issue when trying to replicate this in an imported mesh. I will open a different topic for this.