Unable to specify BC on the mesh in .inp format

mesh
1

Hi,
I have a .inp format mesh file with me which consists of triangular elements. I have successfully imported it into FEniCS. But when I am creating the mesh function for specifying boundary conditions, I am getting the error as “Unable to recognize cell type”. The same thing I have done on .msh format of the same mesh, which worked fine, but in .inp format, I am getting this error.
Here is the .inp file of my mesh

*Heading
 F:\Intelimek\FEniCS\loaded_bar
*NODE
1, 0, 0, 0
2, 2, 0, 0
3, 2, 1, 0
4, 0, 1, 0
5, 0, 0.8, 0
6, 0, 0.6, 0
7, 0, 0.4, 0
8, 0, 0.2, 0
9, 0.22222222222222, 0, 0
10, 0.44444444444445, 0, 0
11, 0.66666666666667, 0, 0
12, 0.88888888888889, 0, 0
13, 1.1111111111111, 0, 0
14, 1.3333333333333, 0, 0
15, 1.5555555555556, 0, 0
16, 1.7777777777778, 0, 0
17, 2, 0.2, 0
18, 2, 0.4, 0
19, 2, 0.6, 0
20, 2, 0.8, 0
21, 1.7777777777778, 1, 0
22, 1.5555555555556, 1, 0
23, 1.3333333333333, 1, 0
24, 1.1111111111111, 1, 0
25, 0.88888888888889, 1, 0
26, 0.66666666666667, 1, 0
27, 0.44444444444444, 1, 0
28, 0.22222222222222, 1, 0
29, 0.16318629087875, 0.49991689477637, 0
30, 1.8371082135119, 0.4997925291569, 0
31, 0.9999732758533, 0.18583150772237, 0
32, 0.77777777777778, 0.80754991027013, 0
33, 1.2223584742912, 0.80787928969703, 0
34, 0.55018403961601, 0.20151035272084, 0
35, 1.44803071747, 0.19866320702855, 0
36, 0.32461259571583, 0.79836728586279, 0
37, 1.6751106883342, 0.79812371677888, 0
38, 1.4458744900672, 0.80678721726964, 0
39, 1.3338119195254, 0.61639853254751, 0
40, 1.1119286235248, 0.61557126164576, 0
41, 1.2232005180586, 0.42447148876223, 0
42, 1.443953855655, 0.42519162532407, 0
43, 1.0008157665817, 0.43705418827036, 0
44, 1.5432390862639, 0.60754271114358, 0
45, 0.55377677312343, 0.80696430844115, 0
46, 0.6651856131325, 0.61712049178473, 0
47, 0.45607529930052, 0.60845503134092, 0
48, 0.55263861999883, 0.42877673329863, 0
49, 0.77308567420635, 0.43269283724509, 0
50, 0.77653324195303, 0.19173155663644, 0
51, 0.31427629654814, 0.17337582499222, 0
52, 1.6857262401802, 0.17276543122905, 0
53, 1.223466758047, 0.19173155663644, 0
54, 0.88813224253719, 0.61965221182112, 0
55, 1.0000328530218, 0.80844211223901, 0
56, 0.16986240997462, 0.3, 0
57, 0.35659805168811, 0.4020747954424, 0
58, 1.8456129213862, 0.69312367342191, 0
59, 0.15439491583137, 0.69318509790067, 0
60, 1.8309828585347, 0.3, 0
61, 1.6427295429972, 0.40092051752582, 0
62, 0.14263487707956, 0.85831824790823, 0
63, 1.8573267116729, 0.85827467739927, 0
64, 0.14866987469895, 0.14173195954204, 0
65, 1.8588973752985, 0.13455308624581, 0
66, 0.28996340549991, 0.60248873096781, 0
67, 1.7100658293649, 0.60248811879862, 0
68, 0.88734502039672, 0.31120219505938, 0
69, 1.3339563924019, 0.30858848017805, 0
70, 0.67291061911725, 0.31225437769288, 0
71, 1.1114596218899, 0.3098132361048, 0
******* E L E M E N T S *************
*ELEMENT, type=T3D2, ELSET=Line1
1, 4, 5
2, 5, 6
3, 6, 7
4, 7, 8
5, 8, 1
*ELEMENT, type=T3D2, ELSET=Line2
6, 1, 9
7, 9, 10
8, 10, 11
9, 11, 12
10, 12, 13
11, 13, 14
12, 14, 15
13, 15, 16
14, 16, 2
*ELEMENT, type=T3D2, ELSET=Line3
15, 2, 17
16, 17, 18
17, 18, 19
18, 19, 20
19, 20, 3
*ELEMENT, type=T3D2, ELSET=Line4
20, 3, 21
21, 21, 22
22, 22, 23
23, 23, 24
24, 24, 25
25, 25, 26
26, 26, 27
27, 27, 28
28, 28, 4
*ELEMENT, type=CPS3, ELSET=Surface1
29, 48, 57, 34
30, 35, 61, 42
31, 34, 57, 51
32, 52, 61, 35
33, 42, 61, 44
34, 47, 57, 48
35, 43, 54, 49
36, 42, 44, 39
37, 15, 52, 35
38, 22, 38, 37
39, 38, 44, 37
40, 49, 54, 46
41, 46, 54, 32
42, 34, 51, 10
43, 22, 37, 21
44, 28, 36, 27
45, 45, 46, 32
46, 36, 47, 45
47, 39, 44, 38
48, 26, 45, 32
49, 16, 52, 15
50, 45, 47, 46
51, 36, 45, 27
52, 11, 34, 10
53, 15, 35, 14
54, 10, 51, 9
55, 41, 42, 39
56, 31, 50, 12
57, 13, 53, 31
58, 47, 48, 46
59, 23, 38, 22
60, 33, 38, 23
61, 11, 50, 34
62, 35, 53, 14
63, 33, 39, 38
64, 40, 41, 39
65, 54, 55, 32
66, 27, 45, 26
67, 12, 50, 11
68, 14, 53, 13
69, 33, 55, 40
70, 40, 55, 54
71, 24, 33, 23
72, 26, 32, 25
73, 13, 31, 12
74, 33, 40, 39
75, 24, 55, 33
76, 32, 55, 25
77, 25, 55, 24
78, 40, 43, 41
79, 40, 54, 43
80, 48, 49, 46
81, 47, 66, 57
82, 61, 67, 44
83, 7, 29, 6
84, 19, 30, 18
85, 8, 56, 7
86, 6, 59, 5
87, 20, 58, 19
88, 18, 60, 17
89, 7, 56, 29
90, 29, 59, 6
91, 19, 58, 30
92, 56, 57, 29
93, 30, 60, 18
94, 30, 61, 60
95, 51, 57, 56
96, 9, 64, 1
97, 4, 62, 28
98, 21, 63, 3
99, 2, 65, 16
100, 60, 61, 52
101, 34, 70, 48
102, 1, 64, 8
103, 5, 62, 4
104, 3, 63, 20
105, 17, 65, 2
106, 36, 66, 47
107, 57, 66, 29
108, 44, 67, 37
109, 30, 67, 61
110, 42, 69, 35
111, 37, 63, 21
112, 28, 62, 36
113, 51, 64, 9
114, 16, 65, 52
115, 49, 68, 43
116, 59, 66, 36
117, 37, 67, 58
118, 68, 70, 50
119, 49, 70, 68
120, 31, 71, 68
121, 68, 71, 43
122, 31, 68, 50
123, 53, 71, 31
124, 50, 70, 34
125, 35, 69, 53
126, 41, 69, 42
127, 69, 71, 53
128, 41, 71, 69
129, 56, 64, 51
130, 48, 70, 49
131, 58, 63, 37
132, 36, 62, 59
133, 43, 71, 41
134, 52, 65, 60
135, 8, 64, 56
136, 20, 63, 58
137, 59, 62, 5
138, 60, 65, 17
139, 29, 66, 59
140, 58, 67, 30
*ELSET,ELSET=fixed
1, 2, 3, 4, 5, 
*ELSET,ELSET=free
15, 16, 17, 18, 19, 
*ELSET,ELSET=top
20, 21, 22, 23, 24, 25, 26, 27, 28, 
*ELSET,ELSET=bottom
6, 7, 8, 9, 10, 11, 12, 13, 14, 
*ELSET,ELSET=fluid
29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 
39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 
49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 
59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 
69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 
79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 
89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 
109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 
119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 
129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 
139, 140, 
*NSET,NSET=fixed
1, 4, 5, 6, 7, 8, 
*NSET,NSET=free
2, 3, 17, 18, 19, 20, 
*NSET,NSET=top
3, 4, 21, 22, 23, 24, 25, 26, 27, 28, 
*NSET,NSET=bottom
1, 2, 9, 10, 11, 12, 13, 14, 15, 16, 
*NSET,NSET=fluid
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 
11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 
31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 
51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 
71,

and here is the code attached

from fenics import *
import matplotlib.pyplot as plt

import meshio
mesh_from_file = meshio.abaqus.read("loaded_bar.inp")
import h5py

import numpy
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    #cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}) #, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
    mvc = MeshValueCollection("size_t", mesh, 2) 

with XDMFFile("facet_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-9-a898e4cdcd2c> in <module>
      1 with XDMFFile("facet_mesh.xdmf") as infile:
----> 2     infile.read(mvc, "name_to_read")
      3 mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to recognise cell type.
*** Reason:  Unknown value "".
*** Where:   This error was encountered inside XDMFFile.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

Thanks

2

You don’t have any subdomain information extracted into your mesh. In fact you have commented out the cell_data field that is required to get information about subdomains/boundaries which you can then pass on to MeshValueCollection/MeshFunction and so on…

See Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #104 by dokken and the comment above it to get an idea on how to do if for your Abaqus mesh.

1 Like
3

I am not getting the solution from the above posts you have given.
The reason why I commented out the cell_data field because there is no Cell data: gmsh:physical, gmsh:geometrical in my mesh_from_file because it is the Abaqus file, that data will be present in the GMSH mesh file.
The other thing which I am not able to follow from the given post Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #104 by dokken
is facet_data = mesh.cell_data_dict["gmsh:physical"]["line"], because this dictionary is empty in my case.
I request you to please provide me the code if you have it in which the Abaqua mesh file is deal with and following line is worked successfully.

with XDMFFile("facet_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

ds = Measure("ds", domain=mesh, subdomain_data=mf)

So that, I will be able to apply boundary conditions.
Thanks in advance

4

This is the exact reason I referenced the above post. Abaqus does not create markers for you to be able to import the mesh in dolfin. You will to construct cell_data by yourself which is of shape say num_line_cells which can then be imported into a MeshFunction which further can be used to assign boundary conditions.

5

I wasn’t on a computer. The following should provide you some help. It seems there are a lot of duplicate cell_sets in your Abaqus mesh.

In [1]: import meshio, numpy as np

In [2]: msh = meshio.abaqus.read("foo.inp")

In [3]: msh
Out[3]: 
<meshio mesh object>
  Number of points: 71
  Number of cells:
    line: 5
    line: 9
    line: 5
    line: 9
    triangle: 112
  Point sets: fixed, free, top, bottom, fluid
  Cell sets: fixed, free, top, bottom, fluid, Line1, Line2, Line3, Line4, Surface1

In [4]: point, tri_cells = msh.points, msh.get_cells_type("triangle")

In [5]: line_cells = msh.get_cells_type("line")

In [6]: num_line_cells = line_cells.shape[0]

In [7]: msh.cell_sets_dict
Out[7]: 
{'fixed': {'line': array([0, 1, 2, 3, 4])},
 'free': {'line': array([14, 15, 16, 17, 18])},
 'top': {'line': array([19, 20, 21, 22, 23, 24, 25, 26, 27])},
 'bottom': {'line': array([ 5,  6,  7,  8,  9, 10, 11, 12, 13])},
 'fluid': {'triangle': array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,
          13,  14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,
          26,  27,  28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,
          39,  40,  41,  42,  43,  44,  45,  46,  47,  48,  49,  50,  51,
          52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,  64,
          65,  66,  67,  68,  69,  70,  71,  72,  73,  74,  75,  76,  77,
          78,  79,  80,  81,  82,  83,  84,  85,  86,  87,  88,  89,  90,
          91,  92,  93,  94,  95,  96,  97,  98,  99, 100, 101, 102, 103,
         104, 105, 106, 107, 108, 109, 110, 111])},
 'Line1': {'line': array([0, 1, 2, 3, 4])},
 'Line2': {'line': array([ 5,  6,  7,  8,  9, 10, 11, 12, 13])},
 'Line3': {'line': array([14, 15, 16, 17, 18])},
 'Line4': {'line': array([19, 20, 21, 22, 23, 24, 25, 26, 27])},
 'Surface1': {'triangle': array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,
          13,  14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,
          26,  27,  28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,
          39,  40,  41,  42,  43,  44,  45,  46,  47,  48,  49,  50,  51,
          52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,  64,
          65,  66,  67,  68,  69,  70,  71,  72,  73,  74,  75,  76,  77,
          78,  79,  80,  81,  82,  83,  84,  85,  86,  87,  88,  89,  90,
          91,  92,  93,  94,  95,  96,  97,  98,  99, 100, 101, 102, 103,
         104, 105, 106, 107, 108, 109, 110, 111])}}

In [8]: boundary_tags = np.zeros(num_line_cells, int)

In [9]: cell_info = msh.cell_sets_dict

In [10]: boundary_tags[cell_info["fixed"]["line"]]=1 # for fixed boundary

In [11]: boundary_tags[cell_info["free"]["line"]]=2 # for free boundary

In [12]: boundary_tags[cell_info["top"]["line"]]=3 # for top boundary

In [13]: boundary_tags[cell_info["bottom"]["line"]]=4 # for bottom  boundary

In [14]: boundary_tags
Out[14]: 
array([1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 2, 2, 2, 2, 2, 3, 3, 3,
       3, 3, 3, 3, 3, 3])

In [15]: line_mesh = meshio.Mesh(point, cells={"line": line_cells}, cell_data={"boundaryMarks": [boundary 
    ...: _tags]})

In [16]: meshio.xdmf.write("linemesh.xdmf", line_mesh)

In [17]: meshio.xdmf.write("mesh.xdmf", meshio.Mesh(point, cells={"triangle": tri_cells}))

You should then be simply able to do as before

from dolfin import *
comm = MPI.comm_world
mesh = Mesh()
with XDMFFile(comm, "mesh.xdmf")  as infile:
    infile.read(mesh)
    mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim() - 1) # you are looking for `1D` entities

with XDMFFile(comm, "linemesh.xdmf") as infile:
    infile.read(mvc, "boundaryMarks")

mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
1 Like
6

Thank you so much. It worked for me.
I just have some small doubts. The code provided by you is working fine but it is creating a 3D mesh, but as my domain is 2D only, I want 2D mesh.
I have found

line_mesh.prune_z_0()
triangle_mesh.prune_z_0()

writing these 2 lines will convert 3D mesh to 2D mesh. Can you just tell me what is the meaning of these lines?
The other thing I want to ask that, I am now moving towards FEniCSx because I want to deal with quad elements. Can you tell me, apart from this tutorial An overview of the FEniCS Project — FEniCSx tutorial , are there any additional resources to learn FEniCSx, and if I have similar doubts regarding importing of meshes, boundary conditions, where can I ask these doubts?

Thank you so much for helping me.

7

This is because the original mesh that you provided from Abaqus had a 0 z-component. The lines that you posted prune that component. It can also be ignored when creating the meshes by taking points[:, :2] instead of points.

As far as dolfinx goes, the tutorial is pretty comprehensive to begin with. For converting meshes, you can refer to the docstrings in meshio. Much of the package is self explanatory. If you have experience with dolfin then maybe trying to reproduce your work with dolfinx is also a good way to get a grasp of things. Good luck!