Hi! I’m a new learner of FEniCS and trying to apply pressure drop boundary conditions to check if the Stokes flow is Darcy flow. I haven’t found any examples in the FEniCS demos. Please help me with where I can start. Any help is much appreciated. Here is my code.

import numpy as np

import time

from dolfin import; from mshr import *4*(x[1]

import dolfin.common.plotting as fenicsplot

from matplotlib import pyplot as plt

L = 4

H = 2

XMIN = 0.0; XMAX = L

YMIN = 0.0; YMAX = H

VIN = 4.0

resolution = 32

mesh = generate_mesh(Rectangle(Point(XMIN, YMIN), Point(L, H)) - Circle(Point(1.5, 0.25 * H), 0.2)

- Circle(Point(0.5, 0.5 * H), 0.2) - Circle(Point(2.0, 0.75 * H), 0.2)

- Circle(Point(3.0, 0.3 * H), 0.15) - Circle(Point(2.5, 0.6 * H), 0.1), resolution)

no_levels = 0

for i in range(0,no_levels):

cell_marker = MeshFunction(“bool”, mesh, mesh.topology().dim())

for cell in cells(mesh):

cell_marker[cell] = False

p = cell.midpoint()

if p.distance(Point(0.5, 0.5 * H)) < 0.5:

cell_marker[cell] = True

mesh = refine(mesh, cell_marker)

plt.figure()

plot(mesh, title=“Mesh”)

plt.show()

VE = VectorElement(“CG”, mesh.ufl_cell(), 2)

QE = FiniteElement(“CG”, mesh.ufl_cell(), 1)

WE = VE * QE

W = FunctionSpace(mesh, WE, constrained_domain=pbc)

V = FunctionSpace(mesh, VE, constrained_domain=pbc)

Q = FunctionSpace(mesh, QE, constrained_domain=pbc)

w = Function(W)

(u, p) = (as_vector((w[0],w[1])), w[2])

(v, q) = TestFunctions(W)

uin = Expression(("VIN(YMAX-x[1]))/(YMAXYMAX)", “0.”), YMAX=YMAX, VIN=VIN, element = V.ufl_element())

ib = Expression(“near(x[0],XMIN) ? 1. : 0.”, XMIN=XMIN, element = Q.ufl_element())

ob = Expression(“near(x[0],XMAX) ? 1. : 0.”, XMAX=XMAX, element = Q.ufl_element())

wb = Expression(“x[0] > XMIN + DOLFIN_EPS && x[0] < XMAX - DOLFIN_EPS ? 1. : 0.”, XMIN=XMIN, XMAX=XMAX, element = Q.ufl_element())

h = CellDiameter(mesh)

C = 1.0e10

gamma = C/h

f = Expression((“0.0”,“0.0”), element = V.ufl_element())

residual = ( - pdiv(v)inner(u, v))*ds - inner(f, v)dx + inner(grad(u), grad(v))inner(u - uin, v) + wbdx + div(u)(ibqdx +

gammadx )

solve(residual == 0, w)

!rm results-NS/

plt.figure()

plot(u, title=“Velocity”)

plt.figure()

plot(p, title=“Pressure”)

plt.show()