How to apply DirichletBC to node set in mesh imported from Abaqus?

Question

I am trying to solve a simple linear elasticity problem using a mesh created
in Abaqus/CAE. The problem consists of a cylinder that is fixed at both ends
and loaded with a uniform internal pressure. I converted the Abaqus input file
to the FEniCS xml format using dolfin-convert.

Now I am having trouble applying the DirichletBC to the top and bottom node
sets of the mesh. Using dolfin-convert The *nset keyword in the Abaqus
input file gets exported as

<mesh_value_collection name="top" type="uint" dim="0" size="75">
  <value cell_index="3270" local_entity="0" value="2"/>
  etc.

How can I apply the DirichletBC to the vertices defined by this node set?

The code I have so far is as follows:

from dolfin import *

# Load mesh
mesh_file = "abaqus/cylinder_mod"
mesh = Mesh(mesh_file+".xml")
interior = MeshFunction("size_t", mesh, mesh_file+"_facet_region.xml")

plot(mesh, title="Mesh")
plot(interior, title="Interior surface")

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

# Define Dirichlet boundary conditions
# DOESN'T WORK
nodesets = MeshFunction("size_t", mesh, 0)
nset_bot = 1
nset_top = 2
c = Constant(("0.0", "0.0", "0.0"))
bc_bot = DirichletBC(V, c, nodesets, nset_bot, "pointwise")
bc_top = DirichletBC(V, c, nodesets, nset_top, "pointwise")
bcs = [bc_bot, bc_top]

# Internal pressure
f = Constant(1.0e6)
ds = Measure("ds")[interior]
n = FacetNormal(mesh)

# Elasticity parameters
E = 2.0e11
nu = 0.3
mu = E/(2.0*(1.0 + nu))
lmbda = E*nu/((1.0 + nu)*(1.0 - 2.0*nu))

# Stress
def sigma(v):
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(len(v))

# Variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = inner(sigma(u), grad(v))*dx
L = inner(f*n, v)*ds

u = Function(V)
solve(a == L, u, bcs)

plot(u, mode="displacement", title="Displacement")
interactive()

When I try to run the code I get the following error:

Traceback (most recent call last):
  File "cylinder.py", line 19, in <module>
    bc_bot = DirichletBC(V, c, nodesets, nset_bot, "pointwise")
  File "/opt/fenics-1.6.0/lib/python2.7/site-packages/dolfin/fem/bcs.py", line 129, in __init__
    cpp.DirichletBC.__init__(self, *args)
  File "/opt/fenics-1.6.0/lib/python2.7/site-packages/dolfin/cpp/fem.py", line 1159, in __init__
    _fem.DirichletBC_swiginit(self,_fem.new_DirichletBC(*args))
RuntimeError:

*** -------------------------------------------------------------------------
*** Error:   Unable to create Dirichlet boundary condition.
*** Reason:  User MeshFunction is not a facet MeshFunction (dimension is wrong).
*** Where:   This error was encountered inside DirichletBC.cpp.
*** Process: 0
***
*** DOLFIN version: 1.6.0

P.S. The geometry for the actual problem that I'm trying to solve is slightly
more complicated so I can't define the subdomains using compile_subdomains
as, for example, in the hyperelasticity demo.

The Abaqus input file is listed below. The converted FEniCS xml files were too long to include here but can be created with dolfin-convert cylinder_mod.inp cylinder_mod.xml.

*Heading
** Job name: cylinder Model name: Model-1
** Generated by: Abaqus/CAE 6.14-1
*Preprint, echo=NO, model=NO, history=NO, contact=NO
**
** PARTS
**
*Part, name=Part-1
*Node
      1,           0., 0.0299999993,  0.100000001
      2,           0., 0.0299999993,           0.
      3,           0., 0.0500000007,           0.
      4,           0., 0.0500000007,  0.100000001
      5, -0.0299999993,           0.,  0.100000001
      6,           0., -0.0299999993,  0.100000001
      7, 0.0299999993,           0.,  0.100000001
      8, -0.0299999993,           0.,           0.
      9,           0., -0.0299999993,           0.
     10, 0.0299999993,           0.,           0.
     11, -0.0500000007,           0.,           0.
     12,           0., -0.0500000007,           0.
     13, 0.0500000007,           0.,           0.
     14, 0.0500000007,           0.,  0.100000001
     15,           0., -0.0500000007,  0.100000001
     16, -0.0500000007,           0.,  0.100000001
     17, -0.0299999993,           0., 0.0500000007
     18,           0., -0.0299999993, 0.0500000007
     19, 0.0299999993,           0., 0.0500000007
     20,           0., 0.0299999993, 0.0500000007
     21, 0.00858615153, 0.0290417653, 0.00903659221
     22, 0.00434918609, 0.0304067396,  0.023572538
*Element, type=C3D4
 1,  3, 19, 13, 14
 2,  1, 19, 20,  4
 3, 18, 16, 17, 11
 4,  6, 16, 17, 18
 5,  9, 18, 12, 13
 6, 13, 10, 21,  2
 7,  2, 17, 11,  8
 8,  6, 16, 18, 15
 9,  9,  8, 11, 18
10, 10, 19,  9, 13
11,  8, 17, 11, 18
12, 20, 17, 11,  2
13,  7, 14, 19,  1
14, 13, 14, 19, 18
15, 11, 18, 12,  9
16, 13, 19,  9, 18
17, 20, 11,  3,  2
18, 13, 21, 22,  3
19,  3,  4, 19, 14
20,  6, 15, 18, 14
21,  1,  4, 20, 16
22, 14, 15, 18, 13
23, 18, 14, 19,  7
24, 14,  6,  7, 18
25, 15, 12, 16, 18
26,  1, 14, 19,  4
27,  2, 21, 22, 10
28, 16, 18, 12, 11
29, 10, 22, 20, 19
30, 20, 16, 17,  5
31,  3, 13, 19, 22
32, 13, 21,  3,  2
33, 10, 21, 22, 13
34, 22, 13, 19, 10
35, 20,  4, 19,  3
36,  4, 11,  3, 20
37, 12, 15, 13, 18
38, 11, 20,  4, 16
39,  5, 16, 17,  6
40,  1,  5, 16, 20
41, 16, 11, 20, 17
42,  3, 21, 22,  2
43, 19, 22, 20,  3
44,  3, 22, 20,  2
45,  2, 22, 20, 10
*Nset, nset=all, generate
  1,  22,   1
*Elset, elset=all, generate
  1,  45,   1
*Nset, nset=EdgeSeeds, generate
  1,  16,   1
*Elset, elset=EdgeSeeds
  1,  5,  6,  7,  8,  9, 10, 13, 15, 17, 19, 20, 21, 22, 24, 25
 26, 28, 32, 36, 37, 38, 39, 40, 45
** Section: Section-1
*Solid Section, elset=all, material=Material-1
,
*End Part
**
**
** ASSEMBLY
**
*Assembly, name=Assembly
**
*Instance, name=Part-1-1, part=Part-1
*End Instance
**
*Nset, nset=bot, instance=Part-1-1
  2,  3,  8,  9, 10, 11, 12, 13
*Elset, elset=bot, instance=Part-1-1
  5,  6,  7,  9, 10, 15, 17, 32
*Nset, nset=top, instance=Part-1-1
  1,  4,  5,  6,  7, 14, 15, 16
*Elset, elset=top, instance=Part-1-1
  8, 13, 20, 21, 24, 26, 39, 40
*Elset, elset=_interior_S1, internal, instance=Part-1-1
  2, 10
*Elset, elset=_interior_S4, internal, instance=Part-1-1
  4, 13, 23, 29, 30, 39, 45
*Elset, elset=_interior_S2, internal, instance=Part-1-1
  7,  9, 11, 12, 40
*Elset, elset=_interior_S3, internal, instance=Part-1-1
 16, 24
*Surface, name=interior, type=ELEMENT
_interior_S1, S1
_interior_S4, S4
_interior_S2, S2
_interior_S3, S3
*End Assembly
**
** MATERIALS
**
*Material, name=Material-1
*Elastic
 2e+11, 0.3
**
** BOUNDARY CONDITIONS
**
** Name: bot Type: Symmetry/Antisymmetry/Encastre
*Boundary
bot, ENCASTRE
** Name: top Type: Symmetry/Antisymmetry/Encastre
*Boundary
top, ENCASTRE
** ----------------------------------------------------------------
**
** STEP: Step-1
**
*Step, name=Step-1, nlgeom=NO
*Static
1., 1., 1e-05, 1.
**
** LOADS
**
** Name: pressure   Type: Pressure
*Dsload
interior, P, 1e+06
**
** OUTPUT REQUESTS
**
*Restart, write, frequency=0
**
** FIELD OUTPUT: F-Output-1
**
*Output, field, variable=PRESELECT
**
** HISTORY OUTPUT: H-Output-1
**
*Output, history, variable=PRESELECT
*End Step

Comments

Hi, I am getting the following error when running dolfin-convert with your input. Dolfin version is 1.6.0.

if l[0].startswith('**'): # Pass over comments
IndexError: list index out of range

---

That's interesting. When I run dolfin-convert with the data listed in my question using

dolfin-convert cylinder_mod.inp cylinder_mod.xml

I get the following output and two files cylinder_mod.xml and cylinder_mod_facet_region.xml are created.

WARNING: unrecognised Abaqus input keyword: *Preprint
WARNING: generation of *elsets not tested.
WARNING: unrecognised Abaqus input keyword: *Solid Section
WARNING: unrecognised Abaqus input keyword: *Assembly
WARNING: unrecognised Abaqus input keyword: *Instance
WARNING: unrecognised Abaqus input keyword: *End Instance
WARNING: unrecognised Abaqus input keyword: *End Assembly
WARNING: unrecognised Abaqus input keyword: *Material
WARNING: unrecognised Abaqus input keyword: *Elastic
WARNING: unrecognised Abaqus input keyword: *Boundary
WARNING: unrecognised Abaqus input keyword: *Boundary
WARNING: unrecognised Abaqus input keyword: *Step
WARNING: unrecognised Abaqus input keyword: *Static
WARNING: unrecognised Abaqus input keyword: *Dsload
WARNING: unrecognised Abaqus input keyword: *Restart
WARNING: unrecognised Abaqus input keyword: *Output
WARNING: unrecognised Abaqus input keyword: *Output
WARNING: unrecognised Abaqus input keyword: *End Step
Expecting 22 vertices
Found all vertices
Expecting 45 cells
Found all cells

---

OK, it is really not that important. What matters is whether the set of vertices can be described as "all vertices connected to certain facets". Is that the case?

---

Yes, in this model the BC could equivalently be applied to a set of facets that include the vertices specified in the node set. After some more searching I found this related question on stackexchange:

http://scicomp.stackexchange.com/questions/7330/fenics-meshfunction-usage

and created surfaces in Abaqus for the ends of the cylinder. Following this I could apply the Dirichlet BC with

facet_regions = MeshFunction("size_t", mesh, mesh_file+"_facet_region.xml")
bc_bot = DirichletBC(V, c, facet_regions, 2)
bc_top = DirichletBC(V, c, facet_regions, 3)

However, in my "real" model the geometry is a bit more complex and the node sets are already defined. It might be possible to create surfaces that include only these nodes but this would require a lot of effort and is not possible in general (e.g. if not all of the nodes belonging to a facet are to be included). Therefore, it would still be preferable to apply the BC directly to the nodes/vertices in the mesh rather than to the facets/surfaces.

Answer 1

Hi, based on our discussion, the following should start you in the right direction

from dolfin import *
import numpy as np

manual = True  # How to set bcs

assert MPI.size(mpi_comm_world()) == 1, 'The code below is intended to run in serial'

mesh = UnitSquareMesh(40, 40)

# Build the set of vertices where bcs should be prescribed
node_set = VertexFunction('size_t', mesh, 0)
bc_boundary = CompiledSubDomain('near(x[0], 0)')
bc_boundary.mark(node_set, 1)
# Only keep vertices marked as 1
node_set = [v.index() for v in SubsetIterator(node_set, 1)]

# Get dofs corresponsing to vertices
V = FunctionSpace(mesh, 'CG', 1)
dof_set = np.array(vertex_to_dof_map(V)[node_set], dtype='intc')

# Assemble the system
u, v = TrialFunction(V), TestFunction(V)
a = inner(grad(u), grad(v))*dx
L = inner(Constant(1), v)*dx
A, b = assemble_system(a, L)

bc_f = Constant(2)  # Boundary value to be prescribed
# Manual application of bcs
if manual:
    # Modif A: zero bc row & set diagonal to 1
    A.ident_local(dof_set)
    A.apply('insert')

    # Modif b: entry in the bc row is taken from bc_f
    bc_values = interpolate(bc_f, V).vector().array()
    b_values = b.array()
    b_values[dof_set] = bc_values[dof_set]
    b.set_local(b_values)
    b.apply('insert')
# Auto at same domain which determined node_set
else:
    bc = DirichletBC(V, bc_f, 'near(x[0], 0)')
    bc.apply(A)
    bc.apply(b)

# Check that auto matches manual
uh = Function(V)
solve(A, uh.vector(), b)
plot(uh)
interactive()

print uh.vector().norm('l2')  

The modification to hadle VectorFunctionSpace is trivial. Extending to parallel requires some work - have a look into DirichletBC and GenericMatrix class to see how to do it.
Note that computing bc_values by interpolation is not very efficient. Should that be the problem you can compute the values at boundary dofs in a same way as DirichletBC does it. Finally, setting bcs with user supplied dof-value pair will become easier when this issue is resolved.

Comments

Thank you MiroK for your very comprehensive answer. I guess there is no easy/built-in way to do this in FEniCS then. In Abaqus I can specify a list of node numbers and apply a constant value BC to all of them. Do you know if something like this could be implemented in FEniCS in the future?

I haven't tested your code yet as I am now using facet regions defined using surfaces in Abaqus. However, it will definitely be useful for other models. In most of these model the boundary values are constant so interpolation should not be required. I will have a look at the classes you mentioned to see if the values can be applied directly to the boundary dofs.

---

I guess there is no easy/built-in way to do this in FEniCS then

I think the above code is easy enough to use

Do you know if something like this could be implemented in FEniCS in the future?

I imagine that once the linked issue is resolved you will be able to specify your dof->value map and call bc.apply. But even then you will have to use something like vertex to dof map to transform list of node numbers.

---

Yes, it is quite easy to use. I just thought that there might be a built-in way for reading the dofs from the file and applying the BCs that way. It appears that your code does that with only a couple of extra lines of code which is very nice.

Thank you again MiroK. Your answer was very helpful.