# Extra argument in def (x,on\_boundary)

**URL:** <https://fenicsproject.discourse.group/t/extra-argument-in-def-x-on-boundary/5658>\
**Category:** Uncategorized\
**Created:** [April 23, 2021, 12:14pm UTC](https://fenicsproject.discourse.group/t/extra-argument-in-def-x-on-boundary/5658 "2021-04-23T12:14:04Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![giannoko](https://avatars.discourse-cdn.com/v4/letter/g/7c8e57/32.png) [@giannoko](https://fenicsproject.discourse.group/u/giannoko)\
**Post date:** [April 23, 2021, 12:14pm UTC](https://fenicsproject.discourse.group/t/extra-argument-in-def-x-on-boundary/5658/1 "2021-04-23T12:14:04Z")

</div>

Hello everyone,

I want to pass an extra argument in

```
def boundary_D(x, on_boundary):
    return on_boundary and (near(x[0], 0, tol) or near(x[0], 1, tol))

```

to finally have this

```
def boundary_D(x, on_boundary,index):
    return on_boundary and (near(x[0], list[index][0], tol) or near(x[0], list[index][1], tol))

```

How can I implement this?

Thank you in advance

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 25, 2021, 3:03pm UTC](https://fenicsproject.discourse.group/t/extra-argument-in-def-x-on-boundary/5658/2 "2021-04-25T15:03:25Z")

</div>

It would be beneficial if you could add a complete minimal example illustrating what you would like to use this function for. And you should also clarify what this index is supposed to be representing. Is it a cell index, facet index or something completely unrelated.

Most likely, you should use a MeshFunction to mark the corresponding facets of your domain.

---

<div class="post-metadata">

**Author:** ![giannoko](https://avatars.discourse-cdn.com/v4/letter/g/7c8e57/32.png) [@giannoko](https://fenicsproject.discourse.group/u/giannoko)\
**Post date:** [April 26, 2021, 9:44am UTC](https://fenicsproject.discourse.group/t/extra-argument-in-def-x-on-boundary/5658/3 "2021-04-26T09:44:29Z")

</div>

Dear dokken,

Firstly I would like to honestly thank you for your help. Let me explain my problem. I want to solve the 2D Navier-Stokes equations in a rectangular domain. I want to strongly impose the pressure at the inlet and the outlet through Dirichlet conditions. I use P2-P0 elements for the velocity and the pressure, respectively, because I will use this code later to an FSI problem where the pressure must be discontinuous. So I found the centroids of the mesh and I store them in the inlet\_DG and outlet\_DG lists for the inlet and outlet respectively. Then I make a loop for the inlet and outlet in order to build the Dirichlet conditions. The problem is when I apply the BCs the index (which is just an iterative parameter in inlet\_DG and outlet\_DG) has the last value of the loop and hence the datum function is called for only one index. I think it would be better to mark those centroids with a MeshFunction and then to apply the Dirichlet BCs on these boundaries, but I do not know if it possible to mark points which do not belong to the mesh (as the pressure dofs are in the center of each element). Here is my code.

```
from dolfin import *
import sys
import time
import matplotlib.pyplot as plt
from numericalMethod import *
import scipy.sparse as sp

#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# Compiler Settings
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

t0 = time.time()
comm = MPI.comm_world
rank = comm.Get_rank()
nprocs = comm.Get_size()
parameters['form_compiler']['optimize'] = True
parameters['form_compiler']['cpp_optimize'] = True
parameters["form_compiler"]["quadrature_degree"] = 3
parameters['form_compiler']["cpp_optimize_flags"] = '-O2 -funroll-loops'
mpi_comm_world=MPI.comm_world
PETScOptions.set("mat_mumps_icntl_14",40.0)
set_log_level(40)

#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# Import Meshes
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(),"./mesh/file.h5","r")
hdf.read(mesh,"/mesh",False)
subdomains = MeshFunction("size_t",mesh,2)
hdf.read(subdomains,"/subdomains")
boundaries = MeshFunction("size_t",mesh,1)
hdf.read(boundaries,"/boundaries")

#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# Boundaries & Subdomains ID
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

boundaryIndex = {}

with open("./mesh/boundaryIndex.txt") as bi:
	for i, line in enumerate(bi):
			fields = line.strip().split()
			boundaryIndex.update( { fields[0] : fields[1]} )				

fluid = int(boundaryIndex.get('fluid'))
inlet = int(boundaryIndex.get('inlet'))
outlet = int(boundaryIndex.get('outlet'))
wall = int(boundaryIndex.get('wall'))

#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# Find centroids of the mesh
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
centroid = []

for cell_no in range(mesh.num_cells()):
	coc = Cell(mesh, cell_no).midpoint()

	coo = []
	coo.append(coc.x())
	coo.append(coc.y())
	coo.append(coc.z())

	centroid.append(coo)

# Store the centroids of the inlet and the outlet

inlet_DG = []
outlet_DG = [] 

for i in range(len(centroid)):
	if (abs(centroid[i][0]-0.0041)<1.e-4):
		inlet_DG.append(i)

for i in range(len(centroid)):
	if (abs(centroid[i][0]-0.9958)<1.e-4):
		outlet_DG.append(i)

#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# Newton Raphson
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
 def NewtonRaphson(res,jac,bcs_,w,corr,rank):

	a_tol , r_tol = 1.0e-8, 1.0e-8
	nIter = 0
	bNorm = 1.0
	xNorm = 1.0
	relax = 1.0

	while (bNorm >= r_tol and xNorm >= a_tol and nIter < 10 ) :

		A = assemble(jac) 

		B = assemble(-res)

		[bc.apply(A,B) for bc in bcs_]

		Solver = LUSolver(A,"mumps")
		Solver.solve(corr.vector(),B)
		

		xNorm = corr.vector().norm('l2')
		bNorm = B.norm('l2')

		w.vector()[:] += relax*corr.vector()[:]
		nIter += 1

		if rank == 0:
			print( ' iter: {0:2d} Res: {1:3.2E} Corr: {2:3.2E}'.format(nIter,bNorm,xNorm))

#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# Function space
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

V = VectorElement('CG' ,mesh.ufl_cell(),2)
P = FiniteElement('DG' ,mesh.ufl_cell(),0)
VP = FunctionSpace(mesh,MixedElement([V,P]))
phi,eta = TestFunctions(VP)
 
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# pin points
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

# Set datum pressure
def datum(x, on_boundary ):
	tol = 1.e-4
	p = near( x[0], centroid[index][0] , tol) and near(x[1], centroid[index][1] , tol) 
	if ( p is True):
		return (p)

#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# BCs
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

pIn = Constant(1.0)
pOut = Constant(0.0)
ns2 = Constant((0.0,0.0))

vFluidWall = DirichletBC(VP.sub(0), ns2 ,boundaries, wall ) 

for index in inlet_DG:

	bcpIn = []
	bcpInH = []

	bcpIn.append( DirichletBC(VP.sub(1), pIn , datum0 ,'pointwise'))
	# homogeneous bcpIn
	bcpInH.append(DirichletBC(VP.sub(1), pOut , datum0 ,'pointwise')) 

for index in outlet_DG:

	bcpOut = []
	bcpOut.append(DirichletBC(VP.sub(1), pOut , datum1,'pointwise'))

bcs = [vFluidWall] + bcpIn + bcpOut
bcsH = [vFluidWall] + bcpInH + bcpOut # homegeneous bcs

W = Function(VP)
[bc.apply(W.vector()) for bc in bcs]
u,p = split(W)

#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
# stresses
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

def tau(X):
	nu = Constant(0.1)
	return nu*(grad(X)+grad(X).T)

def cauchy(X,P):
	return -P*Identity(len(X)) + tau(X)

#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
# weak forms
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

dx = Measure("dx")(subdomain_data=subdomains)
ds = Measure("ds")(subdomain_data=boundaries) # outer boundary
dS = Measure("dS")(subdomain_data=boundaries) # inner boundary

n = FacetNormal(mesh)
Diff = inner( cauchy(u,p) , grad(phi))*dx() 
outFBC = - inner( dot( cauchy(u,p),n ),phi )*ds(outlet) 
inFBC = - inner( dot( cauchy(u,p),n ),phi )*ds(inlet)

Cont = eta*div(u)*dx()
Conv = inner(grad(u)*u, phi)*dx()

Res = Diff + Cont + Conv + outFBC + inFBC
Jac = derivative(Res,W)
Corr = Function(VP)

#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# Solve
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

if rank == 0:
			print( "---------------------------------------------------------------")
			print( "\n")
			print( "Start solving...")
			print( "\n")

NewtonRaphson(Res,Jac,bcsH,W,Corr,rank)

#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# Export the solution
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

ufile = File("results/velocity.pvd")
pfile = File("results/pressure.pvd")

u,p = W.split(True)

u.rename("u","u")
p.rename("p","p")

ufile << u
pfile << p

```

Thank you very much  
G.K.