How to create an array of boundary conditions to later use in an if statement for applying analytic BCs


203
views
0
7 months ago by
Hey everyone I am trying to use Dirichlet BCs on two different boundaries but apply them to a single array

BCs:
def boundary_x0(x, on_boundary):
    return near(x[0],x_min) or near(x[0],x_max)

def boundary_x1(x, on_boundary):
    return near(x[1],y_min) or near(x[1],y_max)

​

I then apply the BC values to an array bc_array

id = DirichletBC(V, Constant(1.0), boundary_x0)
id = DirichletBC(V, Constant(2.0), boundary_x1)

bc_id = Function(V)
bc_id.vector()[:] = 0.0
id.apply(bc_id.vector())
bc_array = np.array(bc_id.vector()[:])
print('bc', bc_array)
But when I go to print the array it is always the second condition applied:
(fenicsproject) -bash-4.1$ python BC_neumann.py
bc [ 2.  2.  2.  2.  2.  2.  2.  2.]
​

I need these values (1 or 2) to use in an if statement in a separate c function to calculate the analytic heat flux on the boundary to impose adiabatic conditions.

Any help would be greatly appreciated!

Best,
Terrence

Community: FEniCS Project

1 Answer


3
7 months ago by

You basically reassign what's stored under the name id with these two lines:

id = DirichletBC(V, Constant(1.0), boundary_x0)
id = DirichletBC(V, Constant(2.0), boundary_x1)

The first BC never has a chance to be used.

Just give your second BC a different name.
MWE:

from fenics import *
import numpy as np

mesh = UnitSquareMesh(4,4)
V = FunctionSpace(mesh,FiniteElement('Lagrange',mesh.ufl_cell(),1))

x_min = 0.
x_max = 1.
y_min = 0.
y_max = 1.

def boundary_x0(x, on_boundary):
    return near(x[0],x_min) or near(x[0],x_max)

def boundary_x1(x, on_boundary):
    return near(x[1],y_min) or near(x[1],y_max)

id0 = DirichletBC(V, Constant(1.0), boundary_x0)
id1 = DirichletBC(V, Constant(2.0), boundary_x1)

bc_id = Function(V)


id0.apply(bc_id.vector())
id1.apply(bc_id.vector())

bc_array = np.array(bc_id.vector()[:])
print(bc_array)
Output:
[ 2.  1.  2.  1.  0.  2.  1.  0.  0.  2.  2.  0.  0.  0.  2.  2.  0.  0.
  1.  2.  0.  1.  2.  1.  2.]

Two notes:
  1. The statement bc_id.vector()[:] = 0.0 is redundant, since bc_id = Function(V) already initializes to an all zero field.
  2. Maybe refrain from using python reserved keywords such as id as variable names
Thank you for the response I copied your code:
#define domain
x_min = y_min = z_min = 0.0
x_max = y_max = 1.0
z_max = 1.98/5.0

#define number of elements
nx = ny = 1
nz = 1

#define number of plot steps
plot_steps=50


#create mesh
mesh = BoxMesh(Point(x_min, y_min, z_min), Point(x_max, y_max, z_max), nx, ny, nz)

# define function space
V = FunctionSpace(mesh, 'P', 1)

# Create boundary markers
def boundary_bottom(x, on_boundary):
    return near(x[2],z_min)


def boundary_x0(x, on_boundary):
    return near(x[0],x_min) or near(x[0],x_max)

def boundary_x1(x, on_boundary):
    return near(x[1],y_min) or near(x[1],y_max)



id0 = DirichletBC(V, Constant(1.0), boundary_x0)
id1 = DirichletBC(V, Constant(2.0), boundary_x1)

bc_id = Function(V)
id0.apply(bc_id.vector())
id1.apply(bc_id.vector())
bc_array = np.array(bc_id.vector()[:])
print('bc', bc_array)

# Initial condition and body forces
​

I still get the same thing:

(fenicsproject) -bash-4.1$ python BC_neumann.py
bc [ 2.  2.  2.  2.  2.  2.  2.  2.]

I really don't know what I am doing wrong except that I am in 3D and have a coarser mesh

written 7 months ago by Terrence  
1
I can't try it out myself at the moment, but try to increase your element numbers.
All your points lie in both your boundary definitions, since your BoxMesh only consists of 8 points.
written 7 months ago by klunkean  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »