diff --git a/models/boneload.py b/models/boneload.py
index e688ac6b5de983480398c58ec99fac817518d6e8..e55fd33f57443099610b0676b904e1cac85e62f9 100644
--- a/models/boneload.py
+++ b/models/boneload.py
@@ -15,6 +15,7 @@
import os
import math
import vtk
+import time
colors = vtk.vtkNamedColors()
# import numpy
# from vtk.util.numpy_support import vtk_to_numpy
@@ -129,7 +130,7 @@ def load_off(meshfile):
return nodes, tris
-def identify_nodes(coords, all_no, all_coords, eps=1e-3):
+def identify_nodes_brutal(coords, all_no, all_coords, eps=1e-3):
"""identify the nodes based on their coordinates "coords"
among the list of coordinates "all_coords"
all_no, given as argument, is an array of tags for the nodes defined in all_coords
@@ -141,6 +142,9 @@ def identify_nodes(coords, all_no, all_coords, eps=1e-3):
returns ntags, a list of node tags of the identified nodes
"""
+ cpu_start = time.perf_counter()
+ print(f'identify_nodes_brutal...')
+
notfound = []
ntags = []
for n in coords: # for all nodes
@@ -161,11 +165,97 @@ def identify_nodes(coords, all_no, all_coords, eps=1e-3):
if not found:
notfound.append(n)
+ cpu_stop = time.perf_counter()
+
+ # print stats
+ print(f'identify_nodes_brutal done in {cpu_stop-cpu_start} seconds')
+
if(len(notfound) != 0):
print('WARNING: these nodes have not been identified:')
for n in notfound:
print(f'\tn={n}')
- # raise Exception(f'ERROR: {len(notfound)} nodes have not been found!')
+ raise Exception(f'ERROR: {len(notfound)} nodes have not been found!')
+ return ntags
+
+
+def identify_nodes_SaP(coords, all_no, all_coords, eps=1e-3):
+ """identify the nodes based on their coordinates "coords"
+ among the list of coordinates "all_coords"
+ all_no, given as argument, is an array of tags for the nodes defined in all_coords
+
+ The algorithm used is "sweep & prune" along X.
+
+ Comparision with the brutal approach: (6567 nodes, 12696 triangles)
+ identify_nodes_brutal done in 26.600111499999997 seconds
+ sweep and prune done in 0.30412259999999947 seconds
+ (48500 tests instead of 818707890)
+
+ returns ntags, a list of node tags of the identified nodes
+ """
+ print('identify_nodes_SaP...')
+ cpu_start = time.perf_counter()
+
+ # define nod structure
+ class Nod:
+ def __init__(self, x, idx=-1): # selected points have a -1 index
+ self.x = x
+ self.idx = idx
+
+ def __str__(self):
+ return f'{self.idx}: {self.x}'
+
+ # build node list
+ nods = []
+ for i, p in enumerate(all_coords):
+ nods.append(Nod(p, all_no[i]))
+ for p in coords:
+ nods.append(Nod(p, -1))
+
+ # sort nodes according to their x coordinate
+ nods.sort(key=lambda nod: nod.x[0])
+
+ ntags = []
+
+ eps = 1e-3 # geometrical tolerance
+ ntests = 0
+ nnods = len(nods)
+ i = 0
+ while i < nnods:
+ n = nods[i]
+ x = n.x[0]
+ i2 = i
+ while True:
+ i2 += 1
+ if i2 == nnods:
+ break
+ n2 = nods[i2]
+ x2 = n2.x[0]
+ if abs(x2-x) > eps: # test X
+ break
+ ntests += 1
+ if n.idx >= 0 and n2.idx >= 0: # nodes from different sets
+ continue
+ if abs(n.x[1]-n2.x[1]) > eps: # test Y
+ continue
+ if abs(n.x[2]-n2.x[2]) > eps: # test Z
+ continue
+ no = max(n.idx, n2.idx)
+ ntags.append(no)
+ break # hyp: a max of 1 pair should be found per node
+ # print(f'{no} is selected')
+ i += 1
+
+ cpu_stop = time.perf_counter()
+
+ # print stats
+ print(f'sweep and prune done in {cpu_stop-cpu_start} seconds')
+ print(f'({ntests} tests instead of {len(coords)*len(all_coords)})')
+ print(f'{len(coords)} to be identified')
+ print(f'{len(ntags)} successfully identified')
+ if(len(coords) != len(ntags)):
+ print('WARNING: some nodes have not been identified!')
+ raise Exception(f'ERROR: in identify_nodes_SaP!')
+
return ntags
@@ -973,82 +1063,38 @@ if __name__ == "__main__":
# print(help(vtk.vtkCellArray))
# import sys; sys.exit()
- # basic implementation of "sweep & prune" algorithm.
+ # test of the "sweep & prune" algorithm.
+
+ npts = 10
# build a list of points with random coordinates
import random
- all_pts = []
- for i in range(10):
- all_pts.append(Pt([random.random(), random.random(), random.random()]))
+ all_coords = []
+ all_no = []
+ for i in range(npts):
+ all_coords.append([random.random(), random.random(), random.random()])
+ all_no.append(i)
# build a list of selected points
- pts = []
- pts.append(all_pts[1])
- pts.append(all_pts[2])
- pts.append(all_pts[5])
- pts.append(all_pts[8])
+ coords = []
+ coords.append(all_coords[1])
+ coords.append(all_coords[2])
+ coords.append(all_coords[5])
+ coords.append(all_coords[8])
# add some variants of the selected pts
- for p in pts[0:3]:
- all_pts.append(Pt([p.x[0], 0.0, 1.0]))
- all_pts.append(Pt([p.x[0], p.x[1], 1.0]))
-
+ for i, p in enumerate(all_coords[0:3]):
+ all_coords.append([p[0], 2.0, 2.0]); all_no.append(npts); npts+=1
+ all_coords.append([p[0], p[1], 2.0]); all_no.append(npts); npts+=1
+
+ # display all points
print(f'all points:')
- for no, p in enumerate(all_pts):
+ for no, p in zip(all_no, all_coords):
print(f'\t{no}: {p}')
- print(f'selected:')
- for p in pts:
+ # display selected coordinates
+ print(f'selected coordinates:')
+ for p in coords:
print(f'\t{p}')
- class Nod:
- def __init__(self, x, idx=-1): # selected points have a -1 index
- self.x = x
- self.idx = idx
-
- def __str__(self):
- return f'{self.idx}: {self.x}'
-
- # build node list
- nods = []
- for no, p in enumerate(all_pts):
- nods.append(Nod(p.x, no))
- for p in pts:
- nods.append(Nod(p.x, -1))
-
- # sort nodes according to their x coordinate
- nods.sort(key=lambda nod: nod.x[0])
-
- print('sorted list:')
- for n in nods:
- print(f'\t{n}')
-
- print('sweep and prune:')
- eps = 1e-3 # geometrical tolerance
- ntests = 0
- nnods = len(nods)
- i = 0
- while i < nnods:
- n = nods[i]
- x = n.x[0]
- i2 = i
- while True:
- i2 += 1
- if i2 == nnods:
- break
- n2 = nods[i2]
- x2 = n2.x[0]
- if abs(x2-x) > eps: # test X
- break
- ntests += 1
- if n.idx >= 0 and n2.idx >= 0: # nodes from different sets
- continue
- if abs(n.x[1]-n2.x[1]) > eps: # test Y
- continue
- if abs(n.x[2]-n2.x[2]) > eps: # test Z
- continue
- no = max(n.idx, n2.idx)
- print(f'{no} is selected')
- i += 1
- print(
- f'sweep and prune done ({ntests} tests instead of {len(pts)*len(all_pts)})')
+ identify_nodes_SaP(coords, all_no, all_coords, eps=1e-3)
diff --git a/models/mandiblemodel.py b/models/mandiblemodel.py
index e86639ed17e37162cf35348e778e45c306c0cfdb..ce67332ea4de8202075545e505630400fb39c835 100644
--- a/models/mandiblemodel.py
+++ b/models/mandiblemodel.py
@@ -314,7 +314,8 @@ def create_group(mesh_file_or_coord, all_no, all_coords, domain, groups, gname=N
else:
raise Exception(f'Data is not understood: {mesh_file_or_coord}')
# identify nodes of the mesh among the given vertices of the main mesh
- ntags = boneload.identify_nodes(nodes, all_no, all_coords)
+ # ntags = boneload.identify_nodes_brutal(nodes, all_no, all_coords)
+ ntags = boneload.identify_nodes_SaP(nodes, all_no, all_coords)
# create a metafor group
groupset = domain.getGeometry().getGroupSet()