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
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.
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.