### How to restrict a Function to boundary dofs

370
views
0
12 months ago by

EDIT

In the MWE below I compute and export to a file the pressure coefficient for a pressure field around a cylinder, which I initialised with an arbitrary field for simplicity. How can I modify this MWE such that it computes cp only along the boundary of the cylinder instead of at all nodes in the entire mesh?

from fenics import *
from mshr import *
import numpy as np

def pressure_coefficient(p, pInf, qInf):
Q = p.function_space()
cp = Function(Q)
cp.assign(project((p - pInf)/qInf, Q))
return cp

# Create mesh
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)

# Define boundaries
cylinder = "on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3"

# Define pressure field and conditions at infinity
pInf = 3.0
qInf = 4.0
p = Function(FunctionSpace(mesh, "P", 1))
p.interpolate(Expression("x[0]*x[0] + x[1]*x[1]", degree=2))

# Compute pressure coefficient and export to file
cp = pressure_coefficient(p, pInf, qInf)
np.savetxt("cp.dat", cp.vector().array())
Community: FEniCS Project
There are ways to avoid projection if you need this for a boundary condition, can you explain more about your problem or provide a small working test case?
written 12 months ago by pf4d
I updated my post with a MWE.
written 12 months ago by Miguel Goncalves

4
12 months ago by
For one thing, you don't need to project, just do the simple array operations :

def pressure_coefficient(p, pInf, qInf):
cp   = Function(p.function_space())
p_a  = p.vector().array()
cp_a = (p_a - pInf) / qInf
cp.vector().set_local(cp_a)
return cp​

If you want to apply the function only over the boundary :

from fenics import *
from mshr import *
import numpy as np

def pressure_coefficient(p, pInf, qInf):
cp   = Function(p.function_space())
p_a  = p.vector().array()
cp_a = (p_a - pInf) / qInf
cp.vector().set_local(cp_a)
return cp

# Create mesh
channel  = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05)
domain   = channel - cylinder
mesh     = generate_mesh(domain, 64)
P1       = FunctionSpace(mesh, "P", 1)

# Define boundaries
exterior = "on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3"

# Define pressure field and conditions at infinity
pInf = 3.0
qInf = 4.0
p    = Function(P1)
p.interpolate(Expression("x[0]*x[0] + x[1]*x[1]", degree=2))

# Compute pressure coefficient and export to file
cp = pressure_coefficient(p, pInf, qInf)

bc = DirichletBC(P1, cp, exterior)

cp_ext = Function(P1)
bc.apply(cp_ext.vector())

np.savetxt("cp.dat",     cp.vector().array())
np.savetxt("cp_ext.dat", cp_ext.vector().array())
​

If you want to compute this only over the exterior regions with the goal of saving CPU time :

from fenics import *
from mshr import *
import numpy as np

# Create mesh
channel  = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05)
domain   = channel - cylinder
mesh     = generate_mesh(domain, 64)
P1       = FunctionSpace(mesh, "P", 1)

# Define boundaries
class Boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] > 0.1 and x[0]<0.3 and x[1]>0.1 and x[1]<0.3

boundary = Boundary()
ff       = FacetFunction('size_t', mesh, 0)
boundary.mark(ff, 1)
boundary_verticies = np.where(ff.array() == 1)[0]

# Define pressure field and conditions at infinity
pInf = 3.0
qInf = 4.0
p    = interpolate(Expression("x[0]*x[0] + x[1]*x[1]", degree=2), P1)

p_a = p.vector().array()

cp_ext = (p_a[boundary_verticies] - pInf) / qInf

np.savetxt("cp_ext.dat", cp_ext)

Would something like
cp_i = project(cp_ext, P1)​
work if you wanted to update the expression on the boundary if this were to be implemented over multiple time steps?
written 6 months ago by Dan Sweeney