### How to restrict a Function to boundary dofs

100

views

0

**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
4 months ago by
pf4d

I updated my post with a MWE.

written
4 months ago by
Miguel Goncalves

### 1 Answer

3

For one thing, you don't need to project, just do the simple array operations :

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

```
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)
```

Please login to add an answer/comment or follow this question.