Problem with passing a matrix to a Dirichlet boundary condition

1

Hello everyone,

I am facing some trouble passing a matrix to a Dirichlet boundary condition. Here is the description of the problem followed by an MWE:

Consider a disk of radius 1. I want to apply a matrix

gamma2 = as_matrix((10.0*(boundary_coordinates[:][0]),0.0*(boundary_coordinates[:][1]))

where the boundary coordinates are the coordinates of the boundary of the disk. Finally, I pass this matrix as a Dirichlet BC to my problem with the following code:

c_out = DirichletBC(V, gamma2, boundary_markers, 1)
bcs=[c_out]

where the boundary of the disk has been tagged 1. This produces the error:

c_out = DirichletBC(V, gamma2, boundary_markers, 1)
    ^
SyntaxError: invalid syntax

Here is an MWE:

C1 = Circle(Point(0, 0),1)


class b1(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 1, eps=1e+06) and on_boundary

class s1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]**2 + x[1]**2 < 1 

mesh = generate_mesh(C1, 5)

boundary_nodes = BoundaryMesh(mesh, "exterior",  True)
boundary_coordinates = boundary_nodes.coordinates()
n = len(boundary_coordinates)

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1,0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(),0)
b1().mark(boundary_markers, 1)
s1().mark(surface_markers, 8)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)

n = FacetNormal(mesh)

V = VectorFunctionSpace(mesh, "Lagrange", 1)

# Define boundary condition
# Define functions
du = TrialFunction(V)            
v  = TestFunction(V)             
u  = Function(V)     

gamma2 = as_matrix((10.0*(boundary_coordinates[:][0]),0.0*(boundary_coordinates[:][1]))            

c_out = DirichletBC(V, gamma2, boundary_markers, 1)

bcs=[c_out]

d = mesh.geometry().dim()

print(d) 
I = Identity(d)             # Identity tensor
F = I + grad(u) 
J_e = det(F)
C_lumen = F.T*F
Ic_lumen = tr(C_lumen)

alpha, mu_lumen = 0.00001, 0.000001
psi_lumen = (mu_lumen)*(Ic_lumen - 3) + (alpha)*(J_e-1)**2 + mu_lumen*ln(J_e) 
Pi_lumen = (psi_lumen*dx(8)) #+ dot(gamma2,u)*ds(1) #- dot(P*(n), u)*ds(1) 

F1 = derivative(Pi_lumen, u, v)

# Compute Jacobian of F
Ju = derivative(F1, u, du)

# Solve variational problem
solve(F1 == 0, u, bcs,
form_compiler_parameters=ffc_options)
    
#mesh_lumen = Mesh(mesh)
ALE.move(mesh,u)

mesh_moved = Mesh(mesh)

plot(mesh_moved)
plt.show()

Any suggestions why this is happening and what should be done?

2

You are missing a closing bracket, i.e.

gamma2 = as_matrix((10.0*(boundary_coordinates[:][0]),0.0*(boundary_coordinates[:][1])))

resolves this error message, but gives you a new one, i.e.:

UFLException: Expecting nested list or tuple of Exprs.

This is because you are trying to use an numpy array as a row in the matrix, and you need to write out each component of the matrix.

3

I am a little confused as what this means. By the above do you mean that I explicitly need to create the components of the matrix using a for loop for each of the nodes?

4

If we look at what boundary_coordinates is, we observe:

In [2]: boundary_coordinates[:][0]                                              
Out[2]: array([ 1.,  0.])

You cannot mix as_matrix and numpy arrays, so you would have to:

gamma_2 = as_matrix(((10.0*boundary_coordinates[0][0], 10.0*boundary_coordinates[1][0]),(0,0)))
5

I have a question on the above code suggested by you.

According to my understanding, in the above code snippet, gamma2 only registers the x and the y coordinate of one node on the boundary and stores it as the (0,0) entry. Is that correct?
So if I want to extend it over the whole boundary nodes, I need a for loop. Isn’t that true?

6

You could do something like

x = SpatialCoordinate(mesh)
gamma2 = as_matrix(((10.0*x[0], 10*x[0]),(0.0*x[1], 0.0*x[1])))

However, I do not understand why you are passing a matrix to the a boundary condition for a VectorFunctionSpace (which has shape (2, )).

7

Here is the reason:

Suppose I have a disk of radius 1. Further suppose, on the boundary of the disk, i.e., the circle, I have 20 nodes. For each of these nodes, I want to define a constant displacement at each of those nodes. For example, for the first node, I want to define the displacement (1.0,1.0) in the x and y direction, similarly for the second node and so on. So for each of the 20 nodes, I have a displacement vector (1.0,1.0). So for 20 nodes, I will have 20 displacement vectors of (1.0,1.0). Next, I want to pass this to a Dirichlet BC as a matrix that records the displacement for each of those nodes. So the matrix will be of size (number of nodes on the circle, 2), like gamma2 = [[1.0,1.0],[1.0,1.0],…,[1.0,1.0]], and solve my variational form.

Please let me know if it’s clear.

8

This seems a bit off. Are you just exemplifying or does each node on your boundary have a displacement of (1., 1.)? If it’s the latter then I don’t think the problem is well posed since all you are doing is a rigid body translation by (1., 1.) unless, of course, you are pinning some other part on your domain, i.e. specifying a zero displacement.

1 Like
9

The input to the boundary condition constructor DirichletBC cannot be the matrix that you are trying to specify. It should be a function in the appropriate function space (or something that can be projected into the function space), which you for instance can create using an UserExpression.

10

Hi @dokken ,

The problem is I only have the displacement at the nodes and not a function that gives the displacement.
I am not sure what you mean by

"which you for instance can create using an UserExpression"

Will it be possible to give an example?

11

Hi @bhaveshshrimali ,

I think the well-posedness of the problem depends on the presence of essential boundary conditions and the properties of the material on which the mesh is defined.

For instance, in the original code, if we replace gamma2 by

gamma2 = Expression(('0.3*sin(x[0])','0.01*cos(x[1])'),degree=0)

and pass it as a Dirichlet BC as

c_out = DirichletBC(V, gamma2, boundary_markers, 1)
bcs=[c_out]

the problem becomes well posed.

That’s what I think. Correct me if I am wrong here.

12

Yes, exactly! So you were exemplifying when you said the displacement is (1., 1.) at all the nodes. Applying a displacement of (1., 1.) everywhere on your boundary will lead to an ill-posed problem. You can clearly see that it corresponds to a rigid translation by (1., 1.) of the entire body. Also, what do you mean by prescribing displacement at each node and not as a function? You always have a function for applied boundary condition. The function can very well be a constant, but it always is a function.

If this is just the vector (function) [10 * x[0], 0. * x[1]], which you can create with a UserExpression as suggested above.

13

Maybe there is something I am not following here. But can you post your complete problem, maybe then it would be easier for everyone to help you

14

Hi @bhaveshshrimali , Thanks for the clarification.

Here is the complete MWE:

C1 = Circle(Point(0, 0),1)


class b1(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 1, eps=1e+06) and on_boundary

class s1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]**2 + x[1]**2 < 1 

mesh = generate_mesh(C1, 5)

boundary_nodes = BoundaryMesh(mesh, "exterior",  True)
boundary_coordinates = boundary_nodes.coordinates()

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1,0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(),0)
b1().mark(boundary_markers, 1)
s1().mark(surface_markers, 8)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)

plot(mesh, linewidth=4, color='#67001f')

AB = assemble(Constant(1)*dx(domain=mesh))
print('Area before =', AB)

# Normal vector to the intima
n = FacetNormal(mesh)

# Create mesh and define function space

V = VectorFunctionSpace(mesh, "Lagrange", 1)

# Define boundary condition
# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 

gamma2 = Expression(('0.3*sin(boundary_coordinates[:,0])','0.01*cos(boundary_coordinates[:,1])'),degree=0)

#gamma2 = Expression(('0.3*sin(x[0])','0.01*cos(x[1])'),degree=0)  # THIS ONE WORKS

c_out = DirichletBC(V, gamma2, boundary_markers, 1)

bcs=[c_out]

d = mesh.geometry().dim()
I = Identity(d)             # Identity tensor
F = I + grad(u) 
J_e = det(F)
C_lumen = F.T*F
Ic_lumen = tr(C_lumen)

alpha, mu_lumen = 0.000001, 0.0000001 

psi_lumen = (mu_lumen)*(Ic_lumen - 3) + (alpha)*(J_e-1)**2 + mu_lumen*ln(J_e) 
Pi_lumen = (psi_lumen*dx(8))              

F1 = derivative(Pi_lumen, u, v)

# Compute Jacobian of F
Ju = derivative(F1, u, du)

# Solve variational problem
solve(F1 == 0, u, bcs,
form_compiler_parameters=ffc_options)
    
ALE.move(mesh,u)

mesh_moved = Mesh(mesh)

plot(mesh_moved, linewidth=1, color='#053061')
plt.show() 

AF = assemble(Constant(1)*dx(mesh))
print('Area after =',AF)

This spits out the error:

gamma2 = Expression(('0.3*sin(boundary_coordinates[:,0])','0.01*cos(boundary_coordinates[:,1])'),degree=0)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/expression.py", line 400, in __init__
    self._cpp_object = jit.compile_expression(cpp_code, params)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/jit.py", line 158, in compile_expression
    expression = compile_class(cpp_data, mpi_comm=mpi_comm)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 173, in compile_class
    raise RuntimeError("Unable to compile C++ code with dijitso")
RuntimeError: Unable to compile C++ code with dijitso
15

what’s your reason for not using this? This is infact applying the boundary condition at each node (since Lagrange of degree 1 has dofs only at the vertices) of the mesh according to it’s coordinate x. Hence sin(x[0]) and cos(x[1]) etc. For a better understanding of the DirichletBC object, you can maybe look at the old docs: https://fenicsproject.org/olddocs/dolfin/2017.2.0/python/programmers-reference/fem/bcs/DirichletBC.html

1 Like
16

The reason I don’t use this is due to my problem. For my problem, each of the nodes on the boundary (my geometry is not necessarily a circle. It is a simple closed analytic curve) has a displacement associated to it and the displacements are different for each of the nodes.
So I was thinking if I can collect the displacements for each of the nodes and pass it as a Dirichlet BC to my solver.

17

One thing you need to understand is that you cannot send in the data in the format you propose.
You need to insert the values you have at each vertex into the function space the boundary condition lives in. There are many methods for doing so.

A method that will work for small problems (not very efficient) is the following:
Here we map the vector (0.1, 0.2) to the dofs at (0, 0), (0, 2) to (0.5, 0.) and (3, 1) to (1, 0)

from dolfin import *
import numpy as np
coord_to_value_map = {(0,0): np.array([0.1,0.2]), (0.5, 0): np.array([0,2]), (1, 0): np.array([3, 1])}


mesh = UnitSquareMesh(2, 2)
V = VectorFunctionSpace(mesh, "CG", 1)


u_bndry = Function(V)

dof_coords = V.tabulate_dof_coordinates()
for coord in coord_to_value_map.keys():
    coord_ = np.asarray(coord)
    indices = np.flatnonzero(np.isclose(coord_, dof_coords).all(axis=1))
    u_bndry.vector()[indices] = coord_to_value_map[coord]


File("u.pvd") << u_bndry

bndry = "near(x[1], 0) && on_boundary"
bc = DirichletBC(V, u_bndry, bndry)