Coverage for /usr/local/lib/python3.8/site-packages/spam/mesh/projection.py: 76%
59 statements
« prev ^ index » next coverage.py v7.2.3, created at 2023-11-22 13:26 +0000
1# Library of SPAM functions for projecting morphological field onto tetrahedral
2# meshes
3# Copyright (C) 2020 SPAM Contributors
4#
5# This program is free software: you can redistribute it and/or modify it
6# under the terms of the GNU General Public License as published by the Free
7# Software Foundation, either version 3 of the License, or (at your option)
8# any later version.
9#
10# This program is distributed in the hope that it will be useful, but WITHOUT
11# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
13# more details.
14#
15# You should have received a copy of the GNU General Public License along with
16# this program. If not, see <http://www.gnu.org/licenses/>.
18"""
19 The ``spam.mesh.projection`` module offers a two functions that enables the construction of three dimensional unstructured meshes (four noded tetrahedra) that embbed heterogeneous data
20 (namely, **subvolumes** and interface **orientation**) based on eihter:
21 - the distance field of an **segmented images** (see ``spaM.filters.distanceField``),
22 - ideal **geometrical objects** (see ``spam.filters.objects``).
24 >>> import spam.mesh
25 >>> spam.mesh.projectField(mesh, fields)
26 >>> spam.mesh.projectObjects(mesh, objects)
28 NOTE
29 ----
31 Objects array conventions:
33 - Sphere: ``[radius, centerX, centerY, centerZ, phase]``
34 - Cylinder: ``[radius, centerX, centerY, centerZ, directionX, directionY, directionZ, phase]``
36 Distance field and orientation convention:
38 - The distance fields inside connected components have positive values
39 - The normal vectors are pointing outside (to negative values)
40 - The subvolumes returned correspond to the outside (negative part)
42"""
44import numpy
45from spam.mesh.meshToolkit import projmorpho
48def _refactor_mesh(mesh):
49 """
50 Helper private function for projection.
51 Remove additional nodes and change connectivity matrix accordingly.
52 """
53 points = mesh["points"]
54 cells = mesh["cells"]
56 # DEBUG: Artificially add unused node
57 # points = numpy.insert(points, 10, [-1.0,-1.0,-1.0], axis=0)
58 # for i, node_number_x4 in enumerate(cells):
59 # for j, node_number in enumerate(node_number_x4):
60 # if cells[i, j] >= 10:
61 # cells[i, j] += 1
63 # test if unused nodes
64 cells_set = set(cells.ravel()) # ordered numbering in connectivity
65 n_points_used = len(cells_set)
67 if n_points_used != points.shape[0]:
68 print("Deleting points before projection")
70 # removing unused nodes
71 points_to_delete = []
72 n_deleted = 0
73 for new, old in enumerate(cells_set):
74 # check if holes in the connectivity matrix
75 if new != old - n_deleted:
76 points_to_delete.append(new + n_deleted)
77 n_deleted += 1
78 # print(new, old, n_deleted, old - n_deleted)
79 points = numpy.delete(points, points_to_delete, axis=0)
81 print("Renumbering connectivity matrix")
82 # renumbering connectivity matrix accordingly
83 new_node_numbering = {old: new + 1 for new, old in enumerate(cells_set)}
84 for i, node_number_x4 in enumerate(cells):
85 for j, node_number in enumerate(node_number_x4):
86 if new_node_numbering[cells[i, j]] != cells[i, j]:
87 cells[i, j] = new_node_numbering[cells[i, j]]
88 del cells_set
90 # check if cell number start by 1 or 0
91 if cells.min() == 0:
92 print("Shift connectivity by 1 to start at 1")
93 cells += 1
95 return points, cells
98def projectObjects(
99 mesh,
100 objects,
101 analyticalOrientation=True,
102 cutoff=1e-6,
103 writeConnectivity=None,
104 vtkMesh=False,
105):
106 """
107 This function project a set of objects onto an unstructured mesh.
108 Each objects can be attributed a phase.
110 Parameters
111 ----------
112 mesh: dict
113 Dictionary containing points coordinates and mesh connectivity.
115 .. code-block:: python
117 mesh = {
118 # numpy array of size number of nodes x 3
119 "points": points,
120 # numpy array of size number of elements x 4
121 "cells": connectivity,
122 }
124 objects: 2D array
125 The list of objects. Each line corresponds to an object encoded as follow:
127 .. code-block:: text
129 - spheres: radius, oZ, oY, oX, phase
130 - cylinders: radius, oZ, oY, oX, nZ, nY, nX, phase
132 analyticalOrientation: bool, default=True
133 Change how the interface's normals are computed:
135 - `True`: as the direction from the centroid of the tetrahedron to the closest highest value of the object distance field (ie, the center of a sphere, the axis of a cylinder).
136 - `False`: as the normals of the facets which points lie on the object surface.
138 cutoff: float, default=1e-6
139 Volume ratio below wich elements with interfaces are ignored and considered being part of a single phase.
141 writeConnectivity: string, default=None
142 When not None, it writes a text file called `writeConnectivity` the
143 list of node and the list of elements
144 which format is:
146 .. code-block:: text
148 COORdinates ! number of nodes
149 nodeId, 0, x, y, z
150 ...
152 ELEMents ! number of elemens
153 elemId, 0, elemType, n1, n2, n3, n4, subVol(-), normalX,
154 normalY, normalZ
155 ...
157 where:
159 - ``n1, n2, n3, n4`` is the connectivity (node numbers).
160 - ``subVol(-)`` is the sub volume of the terahedron inside the inclusion.
161 - ``normalX, normalY, normalZ`` are to components of the interface normal vector.
162 - ``elemType`` is the type of element. Their meaning depends on the their phase.
163 Correspondance can be found in the function output
164 after the key word **MATE** like:
166 .. code-block:: text
168 <projmorpho::set_materials
169 . field 1
170 . . MATE,1: background
171 . . MATE,2: phase 1
172 . . MATE,3: interface phase 1 with
173 background
174 . field 2
175 . . MATE,1: background
176 . . MATE,4: phase 2
177 . . MATE,5: interface phase 2 with
178 background
179 >
181 Sub volumes and interface vector are only relevant for interfaces. Default value is 0 and [1, 0, 0].
183 vtkMesh: bool, default=False
184 Writes the VTK of interpolated fields and materials. The files take the same name as ``writeConnectivity``.
186 Returns
187 -------
188 (nElem x 5) array
189 For each element: ``elemType, subVol(-), normalX, normalY, normalZ`` (see outputFile for details).
191 NOTE
192 ----
193 Only four noded tetrahedra meshes are supported.
195 >>> import spam.mesh
196 >>> # unstructured mesh
197 >>> points, cells = spam.mesh.createCuboid([1, 1, 1], 0.1)
198 >>> mesh = {"points": points, "cells": cells}
199 >>> # one centered sphere (0.5, 0.5, 0.5) of radius 0.2 (phase 1)
200 >>> sphere = [[0.2, 0.8, 0.1, 0.1, 1]]
201 >>> # projection
202 >>> materials = spam.mesh.projectObjects(mesh, sphere, writeConnectivity="mysphere", vtkMesh=True)
204 """
206 # init projmorpho
207 pr = projmorpho(
208 name="spam" if writeConnectivity is None else writeConnectivity,
209 cutoff=cutoff,
210 )
212 # STEP 1: objects
213 objects = numpy.array(objects)
215 # if objects.shape[1] == 5:
216 # swaps = [
217 # [1, 3], # swap ox <-> oz
218 # ]
219 # for a, b in swaps:
220 # # swap axis for origin point
221 # tmp = objects[:, a].copy()
222 # objects[:, a] = objects[:, b]
223 # objects[:, b] = tmp
224 # elif objects.shape[1] == 7:
225 # # swap axis ellipsoids
226 # tmp = objects[:, 0].copy()
227 # objects[:, 0] = objects[:, 2]
228 # objects[:, 2] = tmp
229 # tmp = objects[:, 3].copy()
230 # objects[:, 3] = objects[:, 5]
231 # objects[:, 5] = tmp
232 # elif objects.shape[1] == 8: # cylinders
233 # swaps = [
234 # [1, 3], # swap ox <-> oz
235 # [4, 6], # swap nx <-> nz
236 # ]
237 # for a, b in swaps:
238 # # swap axis for origin point
239 # tmp = objects[:, a].copy()
240 # objects[:, a] = objects[:, b]
241 # objects[:, b] = tmp
243 pr.setObjects(objects)
245 # STEP 2: Mesh
246 points, cells = _refactor_mesh(mesh)
247 pr.setMesh(points, cells.ravel())
249 # STEP 3: projection
250 pr.computeFieldFromObjects()
251 pr.projection(analytical_orientation=analyticalOrientation)
253 # STEP 4: write VTK
254 if writeConnectivity:
255 pr.writeFEAP()
256 if vtkMesh:
257 pr.writeVTK()
258 pr.writeInterfacesVTK()
260 return numpy.array(pr.getMaterials())
263def projectField(mesh, fields, thresholds=[0.0], writeConnectivity=None, vtkMesh=False):
264 """
265 This function project a set of distance fields onto an unstructured mesh.
266 Each distance field corresponds to a phase and the interface between
267 the two phases is set by the thresholds.
269 Parameters
270 ----------
271 mesh: dict
272 Dictionary containing points coordinates and mesh connectivity.
274 .. code-block:: python
276 mesh = {
277 # numpy array of size number of nodes x 3
278 "points": points,
279 # numpy array of size number of elements x 4
280 "cells": connectivity,
281 }
283 fields: list of dicts
284 The fields should be continuous (like a distance
285 field) for a better projection. They are discretised over a regular
286 mesh (lexicographical order) and eahc one corresponds to a phase.
287 The dictionary containing the fields data is defined as follow
289 .. code-block:: python
291 fields = {
292 # coordinates of the origin of the field (3 x 1)
293 "origin": origin,
294 # lengths of fields domain of definition (3 x 1)
295 "lengths": lengths,
296 # list of fields values (list of 3D arrays)
297 "values": [phase1, phase2],
298 }
300 thresholds: list of floats, default=[0]
301 The list of thresholds.
303 writeConnectivity: string, default=None
304 When not None, it writes a text file called `writeConnectivity` the
305 list of node and the list of elements
306 which format is:
308 .. code-block:: text
310 COORdinates ! number of nodes
311 nodeId, 0, x, y, z
312 ...
314 ELEMents ! number of elemens
315 elemId, 0, elemType, n1, n2, n3, n4, subVol(-), normalX,
316 normalY, normalZ
317 ...
319 where:
321 - ``n1, n2, n3, n4`` is the connectivity (node numbers).
322 - ``subVol(-)`` is the sub volume of the terahedron inside the inclusion.
323 - ``normalX, normalY, normalZ`` are to components of the interface normal vector.
324 - ``elemType`` is the type of element. Their meaning depends on the their phase.
325 Correspondance can be found in the function output
326 after the key word **MATE** like:
328 .. code-block:: text
330 <projmorpho::set_materials
331 . field 1
332 . . MATE,1: background
333 . . MATE,2: phase 1
334 . . MATE,3: interface phase 1 with
335 background
336 . field 2
337 . . MATE,1: background
338 . . MATE,4: phase 2
339 . . MATE,5: interface phase 2 with
340 background
341 >
343 Sub volumes and interface vector are only relevant for interfaces. Default value is 0 and [1, 0, 0].
345 vtkMesh: bool, default=False
346 Writes the VTK of interpolated fields and materials. The files take the same name as ``writeConnectivity``.
348 Returns
349 -------
350 (nElem x 5) array
351 For each element: ``elemType, subVol(-), normalX, normalY, normalZ`` (see outputFile for details).
353 NOTE
354 ----
355 Only four noded tetrahedra meshes are supported.
357 >>> import spam.mesh
358 >>> import spam.kalisphera
359 >>> # unstructured mesh
360 >>> points, cells = spam.mesh.createCuboid([1, 1, 1], 0.1)
361 >>> mesh = {"points": points, "cells": cells}
362 >>> # create an image of a sphere and its distance field
363 >>> sphereBinary = numpy.zeros([101] * 3, dtype=float)
364 >>> spam.kalisphera.makeSphere(sphereBinary, [50.5] * 3, 20)
365 >>> sphereField = spam.mesh.distanceField(one_sphere_binary.astype(bool).astype(int))
366 >>> # format field object
367 >>> fields = {
368 >>> # coordinates of the origin of the field (3 x 1)
369 >>> "origin": [0] * 3,
370 >>> # lengths of fields domain of definition (3 x 1)
371 >>> "lengths": [1] * 3,
372 >>> # list of fields
373 >>> "values": [sphereField],
374 >>> }
375 >>> # projection
376 >>> materials = spam.mesh.projectField(mesh, fields, writeConnectivity="mysphere", vtkMesh=True)
378 """
379 # init projmorpho
380 pr = projmorpho(
381 name="spam" if writeConnectivity is None else writeConnectivity,
382 thresholds=thresholds,
383 )
385 # STEP 1: Fields
387 # origin
388 # origin = [_ for _ in reversed(fields["origin"])]
389 origin = [_ for _ in fields["origin"]]
390 # nPoints
391 # nPoints = [n for n in reversed(fields["values"][0].shape)]
392 nPoints = [n for n in fields["values"][0].shape]
393 # field size
394 # lengths = [_ for _ in reversed(fields["lengths"])]
395 lengths = [_ for _ in fields["lengths"]]
396 # ravel fields values (ravel in F direction to swap zyx to xyz)
397 values = [v.flatten(order="F") for v in fields["values"]]
398 # values = [v.flatten(order='C') for v in fields["values"]]
400 pr.setImage(values, lengths, nPoints, origin)
402 # STEP 2: Mesh
403 points, cells = _refactor_mesh(mesh)
404 pr.setMesh(points, cells.ravel())
406 # STEP 3: projection
407 pr.computeFieldFromImages()
408 pr.projection()
410 # STEP 4: write VTK
411 if writeConnectivity:
412 pr.writeFEAP()
413 if vtkMesh:
414 pr.writeVTK()
415 pr.writeInterfacesVTK()
417 return numpy.array(pr.getMaterials())