Locating the inflow boundary

mesh
1

Hello, I want to integrate on the inflow boundary \partial \Omega_{-} of domain Ω = [-1,1]*[-1, 1] , which is defined by \partial \Omega_{-} = \{x\in \Omega : \beta\cdot \mathbf{n}(\mathbf{x}) < 0\}.
\mathbf{n} is the normal vector of the boundary, and \beta(x) is \beta(x) = (2, 1+0.8\sin(8\pi x))^T .
I have difficulty in defining the function used by dolfinx.mesh.locate_entities_boundary. Its last argument is a function like this

def some_facets(x):
    return (x[1] <= 0.51) & np.isclose(x[0], 0)

The function returns a numpy array, but the ufl.FacetNormal gives a FacetNormal object. How to solve this problem ?

2

Hi Smith_Jack,

An option is to use ufl.conditional().

First create dolfinx.mesh.meshtags with your locator function some_facets(). Then, define your boundary integral measure as
ds = ufl.Measure('ds', domain=your_mesh, subdomain_data=your_meshtags)

Let’s say you marked the boundary facets located by some_facets() with facets_tag. If you want to integrate the flux \beta\cdot\mathbf{n} over the inflow part of the boundary, the following code is a possible way to add the flux integral to the linear form L (assuming necessary stuff like L, \beta and \mathbf{n} are already defined previously in the code):

flux = ufl.dot(beta, n) # The flux you want to integrate
cond = ufl.lt(flux, 0.0) # "lower than" condition, returns True where flux < 0 and False elsewhere
conditional_flux = ufl.conditional(cond, flux, 0.0) # Conditional which is equal to flux where "cond" holds True and equal to zero elsewhere
L += conditional_flux * v * ds(facets_tag)

Here v is a test function. In this code, the integral will only take values over the inflow boundary where \beta\cdot\mathbf{n} < 0. On facets of the boundary where \beta\cdot\mathbf{n} \geq 0, conditional_flux will return zero and resultingly the integral over these facets will be zero.

Cheers,
Halvor

3 Likes
3

Thanks for your answer, it works. : D

1 Like
4

This is an implementation of deal.II step-9 by fenics （without mesh refinement）

from mpi4py import MPI
import ufl
import dolfinx
import numpy as np 
from dolfinx.fem.petsc import LinearProblem

nx = 32
nrefine = 6

domain = dolfinx.mesh.create_rectangle(
    MPI.COMM_WORLD, 
    points=((-1.0, -1.0), (1.0, 1.0)),
    n=(nx, nx),
    diagonal=dolfinx.cpp.mesh.DiagonalType.crossed)

V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

x = ufl.SpatialCoordinate(domain)
delta = 0.1
s = 0.1 
x0 = [-0.75, -0.75]

# |x - x0|^2 < s^2
cond_f = ufl.lt((x[0] - x0[0])**2 + (x[1] - x0[1])**2 , s*s)
f = ufl.conditional(cond_f, 1.0/(10*s*s), 0)

beta = ufl.as_vector((2.0, 1 + 0.8*ufl.sin(8.0*ufl.pi*x[0])))
n = ufl.FacetNormal(domain)
flux = ufl.dot(beta, n)
cond = ufl.lt(flux, 0.0)
conditional_flux = ufl.conditional(cond, flux, 0.0)

xn2 = x[0]**2 + x[1]**2
g = ufl.exp(5*(1-xn2))*ufl.sin(16*ufl.pi*xn2)
#def func(x):
#    flag1 = np.logical_or(np.isclose(x[0], 1), np.isclose(x[0], -1))
#    flag2 = np.logical_or(np.isclose(x[1], 1), np.isclose(x[1], -1))
#    return np.logical_or(flag1, flag2)
#fdim = domain.topology.dim - 1
#boundary = dolfinx.mesh.locate_entities_boundary(domain, fdim, func)
#values = np.full_like(boundary, 1, dtype=np.int32)
#facet_tags = dolfinx.mesh.meshtags(domain, fdim, boundary, values)

a = (v + delta*ufl.dot(beta, ufl.grad(v))) * ufl.dot(beta, ufl.grad(u)) * ufl.dx - conditional_flux*v*u*ufl.ds
L = (v + delta*ufl.dot(beta, ufl.grad(v))) * f * ufl.dx - conditional_flux*v*g*ufl.ds
problem = LinearProblem(a, L, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
5

May I ask you another question ?
In the tutorial I am reading, it imposes the Dirichlet condition by wrting it in the weak form (see the second term on both side of the equation):

(\beta \cdot \nabla u_h, v_h + \delta \beta\cdot\nabla v_h)_\Omega-(u_h, \beta\cdot {\mathbf n} v_h)_{\partial\Omega_-}=(f, v_h + \delta \beta\cdot\nabla v_h)_\Omega -(g, \beta\cdot {\mathbf n} v_h)_{\partial\Omega_-}

Why Dirichlet condition is written in this way ? Can we just impose it by setting the value of the solution on the boudary (Just as what we do in the Poisson equation tutorial.) ?

6

Sure, no problem.

Essentially, what has been done is that the equation

(u_h, \beta\cdot\mathbf{n} v_h)_{\partial\Omega_\_} = (g, \beta\cdot\mathbf{n} v_h)_{\partial\Omega_\_}

has been subtracted from the weak form, such that this condition holds when solving for u_h. Since on the inflow boundary \beta\cdot\mathbf{n} < 0, it is important that the term is subtracted, and not added, in order to maintain coercivity of the bilinear form.

Usually though, weak imposition of Dirichlet boundary conditions is done by including these terms scaled with a penalty parameter and a mesh-dependent metric (often the mesh cell diameter), to ensure that the terms scale properly with mesh refinement, ensuring stability of the scheme with respect to the discretization. See e.g. Weak imposition of Dirichlet conditions for the Poisson problem — FEniCSx tutorial and the reference therein to Nitsche’s method, which is a method for imposing Dirichlet boundary conditions weakly.

7

OK, Thanks again ! (:slight_smile:)