Commit b695f90e authored by Emmanuel Roubin's avatar Emmanuel Roubin
Browse files

fix holes in connectivity

parent 8e85f0a6
Pipeline #62643 failed with stages
in 2 minutes and 19 seconds
......@@ -125,6 +125,7 @@ For this reason the threshold is made on the image after the application of a va
Please note that it's also possible to use the variance filter in scipy through the *generic* interface as follows::
import scipy
standardDev = numpy.sqrt(scipy.ndimage.generic_filter(im.astype(float), numpy.var, footprint=structEl))
but it's more than 100 times slower!
......@@ -140,8 +141,8 @@ Then the aggregates are dilated 2 times in order to retrieve their original size
# now aggregates are dilated 2 times
for i in range(2):
aggregates = sfilters.binaryDilation(aggregates)
aggregates = spam.filters.binaryDilation(aggregates)
plt.imshow(aggregates[aggregates.shape[0]//2], cmap="Greys"); plt.title("Horizontal slice of the identified aggregates"); plt.show()
.. figure:: images/tutorial/06-projection/aggregates.png
......@@ -157,7 +158,7 @@ An array with all the phases can now be created where the value of each voxel wi
# mortar -> 1
# pores -> 2
# aggregates -> 3
phases = numpy.ones_like(im).astype('<u1')*sfilters.binaryErosion(mask)
phases = numpy.ones_like(im).astype('<u1')*spam.filters.binaryErosion(mask)
phases[pores==1] = 2
phases[aggregates==1] = 3
......@@ -217,7 +218,7 @@ Step 2.2: Create the Finite Element Mesh
For the projection we need an unstructured 3D mesh made of 4-node tetrahedra.
At this stage we have to consider the physical dimensions to the FE mesh.
We use the module `mesh.unstructured` which is a wrapper of `pygmsh` and `meshio` to create the needed `.msh` file
We use the module `mesh.unstructured` which is a wrapper of `meshio` to create the needed `.msh` file
and a vtk for visualisation::
# specimen dimensions in mm
......@@ -229,10 +230,10 @@ and a vtk for visualisation::
# create the gmsh file needed for the projection
points, connectivity = spam.mesh.createCylinder(centre, radius, height, lcar, gmshFile="cylinder")
# create the vtk file (the mesh can be visualised in paraview for example)
import spam.helpers
spam.helpers.writeUnstructuredVTK(points, connectivity, fileName="cylinder.vtk")
# Number of nodes: 124487
# Number of elements: 716525
......
......@@ -372,7 +372,7 @@ def projectField(mesh, fields, thresholds=[0.0], nSkip=1, writeConnectivity=None
If string: path to the file that contains the unstructured mesh. It currently takes ``gmsh`` format files.
If dict:
- mesh['nodes'] should be an array of the `n` node positions (n x 3)
- mesh['points'] should be an array of the `n` node positions (n x 3)
- mesh['cells'] should be the connectivity matrix of the `m` elements (m x 4)
fields: array of string or dictionary
......
......@@ -1141,12 +1141,12 @@ void projmorpho::parser( const std::string& config ) {
_n_field = n_field; // number of nodes
_d_field = d_field; // total length of the cube
_v_field = v_field; // fields values
fill_c_field();
std::cout << ".\t field size:\t " << _d_field[0] << " x " << _d_field[1] << " x " << _d_field[2] << std::endl;
std::cout << ".\t field origin:\t " << _o_field[0] << " x " << _o_field[1] << " x " << _o_field[2] << std::endl;
std::cout << ".\t field nodes:\t " << _n_field[0] << " x " << _n_field[1] << " x " << _n_field[2] << " = " << _n_field[0]*_n_field[1]*_n_field[2] << std::endl;
fill_c_field();
unsigned n = _v_field[0].size();
for( unsigned int f=0; f<_v_field.size(); f++ ) {
std::cout << ".\t field " << f+1 << std::endl;
......
......@@ -110,9 +110,6 @@ def createCuboid(lengths, lc, origin=[0., 0., 0.], verbosity=0, gmshFile=None, v
lx, ly, lz = lengths
ox, oy, oz = origin
print(f'lengths: {lengths}')
print(f'origin: {origin}')
# https://gmsh.info/doc/texinfo/gmsh.pdf
# initialize mesh
......@@ -156,7 +153,7 @@ def createCuboid(lengths, lc, origin=[0., 0., 0.], verbosity=0, gmshFile=None, v
if gmshFile is not None:
gmsh.write(f"{gmshFile}.msh")
if vtkFile is not None:
gmsh.write(f"{gmshFile}.vtk")
gmsh.write(f"{vtkFile}.vtk")
# finish here if no output
if skipOutput:
......@@ -165,7 +162,7 @@ def createCuboid(lengths, lc, origin=[0., 0., 0.], verbosity=0, gmshFile=None, v
# Get connectivity and node coordinates
element_tags, node_tags_element = gmsh.model.mesh.getElementsByType(4)
node_tags_nodes, coord, _ = gmsh.model.mesh.getNodes()
node_tags, coord, _ = gmsh.model.mesh.getNodesByElementType(4)
# DEBUG: gmsh GUI
# gmsh.fltk.run()
......@@ -177,18 +174,38 @@ def createCuboid(lengths, lc, origin=[0., 0., 0.], verbosity=0, gmshFile=None, v
# 1. the coordinates array must be switched back to zyx
# 2. the connectivity array must change so that the nodes of each tetrahedron are numbered
# in such a way that the Jacobian is positive. Which means a permutation of the node numbering.
# 3. in some cases the connectivity matrix has some holes in the node numerotations (some nodes are not used)
# it needs to be rebuild for the projection
# get number of nodes and elements
n_elements = element_tags.shape[0]
nodes_set = set(node_tags)
n_nodes = len(nodes_set)
# get new node numbering (without holes)
new_node_numbering = {}
for i, n in enumerate(nodes_set):
new_node_numbering[n] = i
# reshape coord in nx3 / view it ordered (because node_tags_nodes is not ordered) / switch x&z
nodes = coord.reshape(node_tags_nodes.shape[0], 3)[node_tags_nodes.argsort()][:, ::-1]
# reshape the connectivity matrix from flatten gmsh output
connectivity = node_tags.reshape((n_elements, 4))
# create nodes matrix
nodes = numpy.zeros((n_nodes, 3))
for i, (nodes_4x1, coord_4x3) in enumerate(zip(connectivity, coord.reshape((n_elements, 4, 3)))):
# change connectivity with new orderning
connectivity[i] = [new_node_numbering[n] for n in nodes_4x1]
# fill node vector with new connectivity numbering and switch x&z
for j, n in enumerate(connectivity[i]):
nodes[n] = coord_4x3[j][::-1]
# rearange connectivity
# /!\ node_tags_element starts at 1
connectivity = node_tags_element.reshape(element_tags.shape[0], 4) - 1
_ = connectivity.copy()
connectivity[:, 1] = _[:, 3]
connectivity[:, 3] = _[:, 1]
# DEBUG: to check connectivity
# DEBUG: to check connectivity
# import meshio
# meshio.write_points_cells("{}_2.vtk".format(vtkFile), nodes, {'tetra': connectivity}, file_format='vtk')
......@@ -284,11 +301,11 @@ def createCylinder(centre, radius, height, lc, verbosity=0, gmshFile=None, vtkFi
gmsh.option.setNumber("Mesh.Binary", binary)
# Create disk surface
p0 = gmsh.model.geo.addPoint(cx, cy, lc, 1)
p1 = gmsh.model.geo.addPoint(cx + r, cy, lc, 2)
p2 = gmsh.model.geo.addPoint(cx, cy + r, lc, 3)
p3 = gmsh.model.geo.addPoint(cx - r, cy, lc, 4)
p4 = gmsh.model.geo.addPoint(cx, cy - r, lc, 5)
p0 = gmsh.model.geo.addPoint(cx, cy, 0, lc, 1)
p1 = gmsh.model.geo.addPoint(cx + r, cy, 0, lc, 2)
p2 = gmsh.model.geo.addPoint(cx, cy + r, 0, lc, 3)
p3 = gmsh.model.geo.addPoint(cx - r, cy, 0, lc, 4)
p4 = gmsh.model.geo.addPoint(cx, cy - r, 0, lc, 5)
c1 = gmsh.model.geo.addCircleArc(p1, p0, p2)
c2 = gmsh.model.geo.addCircleArc(p2, p0, p3)
c3 = gmsh.model.geo.addCircleArc(p3, p0, p4)
......@@ -309,7 +326,7 @@ def createCylinder(centre, radius, height, lc, verbosity=0, gmshFile=None, vtkFi
if gmshFile is not None:
gmsh.write(f"{gmshFile}.msh")
if vtkFile is not None:
gmsh.write(f"{gmshFile}.vtk")
gmsh.write(f"{vtkFile}.vtk")
# finish here if no output
if skipOutput:
......@@ -317,8 +334,8 @@ def createCylinder(centre, radius, height, lc, verbosity=0, gmshFile=None, vtkFi
return None
# Get connectivity and node coordinates
element_tags, node_tags_element = gmsh.model.mesh.getElementsByType(4)
node_tags_nodes, coord, _ = gmsh.model.mesh.getNodes()
element_tags, node_tags = gmsh.model.mesh.getElementsByType(4)
node_tags, coord, _ = gmsh.model.mesh.getNodesByElementType(4)
# DEBUG: gmsh GUI
# gmsh.fltk.run()
......@@ -330,13 +347,33 @@ def createCylinder(centre, radius, height, lc, verbosity=0, gmshFile=None, vtkFi
# 1. the coordinates array must be switched back to zyx
# 2. the connectivity array must change so that the nodes of each tetrahedron are numbered
# in such a way that the Jacobian is positive. Which means a permutation of the node numbering.
# 3. in some cases the connectivity matrix has some holes in the node numerotations (some nodes are not used)
# it needs to be rebuild for the projection
# get number of nodes and elements
n_elements = element_tags.shape[0]
nodes_set = set(node_tags)
n_nodes = len(nodes_set)
# get new node numbering (without holes)
new_node_numbering = {}
for i, n in enumerate(nodes_set):
new_node_numbering[n] = i
# reshape the connectivity matrix from flatten gmsh output
connectivity = node_tags.reshape((n_elements, 4))
# create nodes matrix
nodes = numpy.zeros((n_nodes, 3))
for i, (nodes_4x1, coord_4x3) in enumerate(zip(connectivity, coord.reshape((n_elements, 4, 3)))):
# change connectivity with new orderning
connectivity[i] = [new_node_numbering[n] for n in nodes_4x1]
# reshape coord in nx3 / view it ordered (because node_tags_nodes is not ordered) / switch x&z
nodes = coord.reshape(node_tags_nodes.shape[0], 3)[node_tags_nodes.argsort()][:, ::-1]
# fill node vector with new connectivity numbering and switch x&z
for j, n in enumerate(connectivity[i]):
nodes[n] = coord_4x3[j][::-1]
# rearange connectivity
# /!\ node_tags_element starts at 1
connectivity = node_tags_element.reshape(element_tags.shape[0], 4) - 1
_ = connectivity.copy()
connectivity[:, 1] = _[:, 3]
connectivity[:, 3] = _[:, 1]
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment