Dear community,
While transitioning from dolfin to dolfinx, I’m trying to solve a simple 3D quasi-static nonlinear elasticity problem of a cube being fixed on one end and loaded on the other.
For testing, I just try to reproduce a zero solution by fixing all DOFs on one end to zero and doing nothing else.
I read my surface mesh with
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with XDMFFile(comm,'path/input_boundary.xdmf','r',encoding=encoding) as infile:
mt_s = infile.read_meshtags(mesh, "Grid_surf")
The problem: Setting a Dirichlet BC on Face 1 like
dbc = DirichletBC(zero, locate_dofs_topological(V_u, 2, mt_s.indices[mt_s.values == 1]))
with
zero = Function(V_u)
zero.vector.set(0.0)
yields an exploding solution! My Newton scheme output is
Predictor (disp 2-norm 0.0000e+00) yields residual 2-norm:
6.4099e-16
iter struct res 2-norm struct incr 2-norm
1 2.4222e-13 3.7008e-16
### TIME STEP 1 / 10 successfully completed, TIME: 0.1000
-----------------------------------------------------------
Predictor (disp 2-norm 3.7008e-16) yields residual 2-norm:
2.4222e-13
iter struct res 2-norm struct incr 2-norm
1 9.9914e-11 1.8821e-13
### TIME STEP 2 / 10 successfully completed, TIME: 0.2000
-----------------------------------------------------------
Predictor (disp 2-norm 1.8784e-13) yields residual 2-norm:
9.9914e-11
iter struct res 2-norm struct incr 2-norm
1 4.0749e-08 7.7227e-11
2 1.6678e-05 3.1718e-08
3 6.8392e-03 1.3034e-05
4 2.8155e+00 5.3572e-03
5 3.9269e+02 2.2093e+00
6 1.4504e+14 3.6851e+02
7 nan 1.1696e+14
8 nan 0.0000e+00
9 nan 0.0000e+00
10 nan 0.0000e+00
Important: Eliminating the DBC produces the intended behavior: nothing happens and the res is always zero to machine precision.
I have tried a lot already and have the feeling that there may be a problem in locate_dofs_topological or my specific DBC syntax.
Since my Newton scheme is well tested with dolfin I’m pretty sure that the problem does not lie within the nonlinear solver loop. The only thing I had to change was the assemble. Now, in dolfinx I’m calling
K = assemble_matrix(Jac, dbcs)
K.assemble()
r = assemble_vector(-Res)
compared to the dolfin version
K, r = assemble_system(Jac, -Res, dbcs)
(The DBC from above is appended to the array dbcs.)
For reproducibility, I will add my xdmf domain and boundary file.
This is the domain file:
<Xdmf Version="3.0">
<Domain>
<Grid Name="Grid">
<Geometry GeometryType="XYZ">
<DataItem DataType="Float" Dimensions="14 3" Format="XML" Precision="8">
0.0 0.0 1.0
0.0 0.0 0.0
0.0 1.0 1.0
0.0 1.0 0.0
1.0 0.0 1.0
1.0 0.0 0.0
1.0 1.0 1.0
1.0 1.0 0.0
0.0 0.5 0.5
1.0 0.5 0.5
0.5 0.0 0.5
0.5 1.0 0.5
0.5 0.5 0.0
0.5 0.5 1.0
</DataItem>
</Geometry>
<Topology NodesPerElement="4" NumberOfElements="24" TopologyType="tetrahedron">
<DataItem DataType="Int" Dimensions="24 4" Format="XML" Precision="4">
9 11 10 12
8 13 11 10
11 10 13 9
8 10 11 12
1 0 8 10
0 2 8 13
10 0 13 4
3 11 8 2
1 8 3 12
11 13 2 6
4 13 9 6
6 11 9 7
10 4 9 5
11 7 3 12
12 9 7 5
12 1 10 5
13 8 0 10
8 2 11 13
10 8 1 12
11 3 8 12
7 11 9 12
13 9 10 4
11 9 13 6
12 10 9 5
</DataItem>
</Topology>
</Grid>
</Domain>
</Xdmf>
This is the BC file with the mesh tags:
<Xdmf Version="3.0">
<Domain>
<Grid Name="Grid_surf">
<Topology NodesPerElement="3" NumberOfElements="24" TopologyType="triangle">
<DataItem DataType="Int" Dimensions="24 3" Format="XML" Precision="4">
1 0 8
0 2 8
3 1 8
2 3 8
0 1 10
4 0 10
1 5 10
5 4 10
1 3 12
5 1 12
3 7 12
7 5 12
5 9 4
4 9 6
7 9 5
6 9 7
2 11 3
6 11 2
3 11 7
7 11 6
0 13 2
4 13 0
2 13 6
6 13 4
</DataItem>
</Topology>
<Attribute AttributeType="Scalar" Center="Face" Name="boundaries_surf">
<DataItem DataType="Int" Dimensions="24" Format="XML" Precision="4">
1
1
1
1
2
2
2
2
3
3
3
3
4
4
4
4
5
5
5
5
6
6
6
6
</DataItem>
</Attribute>
</Grid>
<Grid Name="Grid_edge">
<Topology NodesPerElement="2" NumberOfElements="2" TopologyType="Polyline">
<DataItem DataType="Int" Dimensions="2 2" Format="XML" Precision="4">
1 3
1 0
</DataItem>
</Topology>
<Attribute AttributeType="Scalar" Center="Edge" Name="boundaries_edge">
<DataItem DataType="Int" Dimensions="2" Format="XML" Precision="4">
1
2
</DataItem>
</Attribute>
</Grid>
<Grid Name="Grid_point">
<Topology NodesPerElement="1" NumberOfElements="1" TopologyType="Polyvertex">
<DataItem DataType="Int" Dimensions="1 1" Format="XML" Precision="4">
1
</DataItem>
</Topology>
<Attribute AttributeType="Scalar" Center="Node" Name="boundaries_point">
<DataItem DataType="Int" Dimensions="1" Format="XML" Precision="4">
1
</DataItem>
</Attribute>
</Grid>
</Domain>
</Xdmf>
Hope someone has an idea what’s going wrong with the DBC application.
Bw,
Marc