Set domain length as a variable

1

Hi all.

I want to make a parametrisation study of a channel length (L). I am defining length as L=Constant(value) where value is a number at the beginning of my code.

Then, when I need to define the boundaries I have `outflow = ‘near(x[0], L)’.

This gives me the following error:
------------------- Start compiler output ------------------------
/tmp/tmp8dma58xo/dolfin_subdomain_351ee9b816aebffd744e4e69091dcb79.cpp: In member function ‘virtual bool dolfin::dolfin_subdomain_351ee9b816aebffd744e4e69091dcb79::inside(Eigen::Ref<const Eigen::Matrix<double, -1, 1> >, bool) const’:
/tmp/tmp8dma58xo/dolfin_subdomain_351ee9b816aebffd744e4e69091dcb79.cpp:63:28: error: ‘L’ was not declared in this scope
63 | return near(x[0], L);
| ^
------------------- End compiler output ------------------------

I think it has to do with C++ strings and variables but I haven’t found a way to write it.

Any suggestions? Thanks in advance :slight_smile:

2

Please supply a minimal working code example,following the guidelines of: Read before posting: How do I get my question answered?

3

Sorry dokken, I didn’t think it was necessary for this question. Please find an example below:

import dolfin
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter
import matplotlib
from matplotlib import rc # use LaTeX text

# Plotting parameters
matplotlib.rcParams['contour.negative_linestyle'] = 'solid'
matplotlib.rcParams['text.usetex'] = True
rc('text', usetex=True)
rc('font', family = 'serif')


# Problem parameters
T = 10.0            # final time
num_steps = 500     # number of time steps
dt = T / num_steps  # time step size
mu = 0.01           # kinematic viscosity [N s/m2]
rho = 1             # density [kg/m3]
L = Constant(2)     # channel length
H = Constant(1)     # channel height
r = 0.25            # cell radius 
x0,y0 = -0.5,0.5    # cell centre

# Define geometries
domain = Rectangle(Point(-1,0),Point(L,H))
circles = [[r,x0, y0]]; circles = np.array(circles)

for circle in circles:     
circleD = Circle(Point(circle[1:]), circle[0], 60)

# Set subdomains
domain.set_subdomain(1,circleD)

# Create mesh
mesh=generate_mesh(domain,30)
coor = mesh.coordinates()
xcoor, ycoor = coor[:,0], coor[:,1]

# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

# Define boundaries
inflow  = 'near(x[0], -1)'
outflow = 'near(x[0], L)'
walls   = 'near(x[1], 0) || near(x[1], 1)'

# Define boundary conditions
bcu_noslip  = DirichletBC(V, Constant((0, 0)), walls)
bcp_inflow  = DirichletBC(Q, Constant(8), inflow)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_noslip]
bcp = [bcp_inflow, bcp_outflow]

Apologies again. Thanks!

4

Consider using the following:

class Outflow(SubDomain):
    def inside(self, x, on_bounary):
        return np.isclose(x[0], float(L))

outflow = Outflow()
2 Likes
5

That is great, thank you very much! I used near instead of np.isclose just consistency with my other commands.

I am new to Fenics, so could you please break-up for me the components of this class? In particular, the use of the keyword SudDomain and the command(?) on_boundary.

6

See for instance page 89 of the FEniCS tutorial

7

Thank you.

Don’t I need to explicitly say return on_boundary too?

8

In this case you do not need it, as there are no other possible dofs close to this point and are not on the boundary. In general you should add an and on_boundary

1 Like