unvhandler fixes and freeform performances improved (#91)

* improved performance in freeform (~25%)
* fixed bug and improved flexibility in unvhandler
* new write comparison in unv tests

Author: sc-ita

pygem/freeform.py

@@ -67,8 +67,8 @@ class FFD(object):
 	>>> original_mesh_points = np.load('tests/test_datasets/meshpoints_sphere_orig.npy')
 	>>> free_form = ffd.FFD(ffd_parameters, original_mesh_points)
 	>>> free_form.perform()
-	>>> new_mesh_points = free_form.modified_mesh_points		
-    """
+	>>> new_mesh_points = free_form.modified_mesh_points
+	"""
 	def __init__(self, ffd_parameters, original_mesh_points):
 		self.parameters = ffd_parameters
 		self.original_mesh_points = original_mesh_points
@@ -81,22 +81,31 @@ def perform(self):
 		it sets `self.modified_mesh_points`.
 
 		.. todo::
-			In order to improve the performances, we need to perform the FFD only on the points inside
-			the lattice.
-		"""
-		(n_rows_mesh, n_cols_mesh) = self.original_mesh_points.shape
 
+		"""
 		# translation and then affine transformation
 		translation = self.parameters.origin_box
 
-		fisical_frame = np.array([self.parameters.position_vertex_1 - translation, \
-			  self.parameters.position_vertex_2 - translation, \
-			  self.parameters.position_vertex_3 - translation, \
-			  [0, 0, 0]])
+		physical_frame = np.array([self.parameters.position_vertex_1 - translation,
+								   self.parameters.position_vertex_2 - translation,
+								   self.parameters.position_vertex_3 - translation,
+								   [0, 0, 0]])
 		reference_frame = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1], [0, 0, 0]])
 
-		transformation = at.affine_points_fit(fisical_frame, reference_frame)
-		inverse_transformation = at.affine_points_fit(reference_frame, fisical_frame)
+		transformation = at.affine_points_fit(physical_frame, reference_frame)
+		inverse_transformation = at.affine_points_fit(reference_frame, physical_frame)
+
+		# apply transformation to original mesh points
+		reference_frame_mesh_points = self._transform_points(self.original_mesh_points - translation, transformation)
+
+		# select mesh points inside bounding box
+		mesh_points = reference_frame_mesh_points[(reference_frame_mesh_points[:,0] >= 0.) &
+												  (reference_frame_mesh_points[:,0] <= 1.) &
+												  (reference_frame_mesh_points[:,1] >= 0.) &
+												  (reference_frame_mesh_points[:,1] <= 1.) &
+												  (reference_frame_mesh_points[:,2] >= 0.) &
+												  (reference_frame_mesh_points[:,2] <= 1.)]
+		(n_rows_mesh, n_cols_mesh) = mesh_points.shape
 
 		# Initialization. In order to exploit the contiguity in memory the following are transposed
 		(dim_n_mu, dim_m_mu, dim_t_mu) = self.parameters.array_mu_x.shape
@@ -105,22 +114,19 @@ def perform(self):
 		bernstein_z = np.zeros((dim_t_mu, n_rows_mesh))
 		shift_mesh_points = np.zeros((n_cols_mesh, n_rows_mesh))
 
-		reference_frame_mesh_points = self._transform_points(self.original_mesh_points - translation, \
-															 transformation)
-
 		for i in range(0, dim_n_mu):
-			aux1 = np.power((1-reference_frame_mesh_points[:, 0]), dim_n_mu-1-i)
-			aux2 = np.power(reference_frame_mesh_points[:, 0], i)
+			aux1 = np.power((1-mesh_points[:, 0]), dim_n_mu-1-i)
+			aux2 = np.power(mesh_points[:, 0], i)
 			bernstein_x[i, :] = special.binom(dim_n_mu-1, i) * np.multiply(aux1, aux2)
 
 		for i in range(0, dim_m_mu):
-			aux1 = np.power((1-reference_frame_mesh_points[:, 1]), dim_m_mu-1-i)
-			aux2 = np.power(reference_frame_mesh_points[:, 1], i)
+			aux1 = np.power((1-mesh_points[:, 1]), dim_m_mu-1-i)
+			aux2 = np.power(mesh_points[:, 1], i)
 			bernstein_y[i, :] = special.binom(dim_m_mu-1, i) * np.multiply(aux1, aux2)
 
 		for i in range(0, dim_t_mu):
-			aux1 = np.power((1-reference_frame_mesh_points[:, 2]), dim_t_mu-1-i)
-			aux2 = np.power(reference_frame_mesh_points[:, 2], i)
+			aux1 = np.power((1-mesh_points[:, 2]), dim_t_mu-1-i)
+			aux2 = np.power(mesh_points[:, 2], i)
 			bernstein_z[i, :] = special.binom(dim_t_mu-1, i) * np.multiply(aux1, aux2)
 
 
@@ -138,18 +144,22 @@ def perform(self):
 		shift_mesh_points[0, :] += aux_x
 		shift_mesh_points[1, :] += aux_y
 		shift_mesh_points[2, :] += aux_z
-
-		# Splitting points inside and outside the lattice: TODO not very efficient
-		for i in range(0, n_rows_mesh):
-			if (reference_frame_mesh_points[i, 0] < 0) or (reference_frame_mesh_points[i, 1] < 0) \
-			or (reference_frame_mesh_points[i, 2] < 0) or (reference_frame_mesh_points[i, 0] > 1) \
-			or (reference_frame_mesh_points[i, 1] > 1) or (reference_frame_mesh_points[i, 2] > 1):
-				shift_mesh_points[:, i] = 0
-
-		# Here shift_mesh_points needs to be transposed to be summed with reference_frame_mesh_points
-		self.modified_mesh_points = self._transform_points(np.transpose(shift_mesh_points) + \
-									reference_frame_mesh_points, inverse_transformation) + translation
-
+	
+		# shift_mesh_points needs to be transposed to be summed with mesh_points
+		# apply inverse transformation to shifted mesh points
+		new_mesh_points = self._transform_points(np.transpose(shift_mesh_points) +
+												 mesh_points, inverse_transformation) + \
+						  translation
+
+		# merge non-shifted mesh points with shifted ones
+		self.modified_mesh_points = np.copy(self.original_mesh_points)
+		self.modified_mesh_points[(reference_frame_mesh_points[:,0] >= 0.) &
+								  (reference_frame_mesh_points[:,0] <= 1.) &
+								  (reference_frame_mesh_points[:,1] >= 0.) &
+								  (reference_frame_mesh_points[:,1] <= 1.) &
+								  (reference_frame_mesh_points[:,2] >= 0.) &
+								  (reference_frame_mesh_points[:,2] <= 1.)] \
+			= new_mesh_points
 
 	@staticmethod
 	def _transform_points(original_points, transformation):
@@ -170,3 +180,4 @@ def _transform_points(original_points, transformation):
 			modified_points[i, :] = transformation(original_points[i])
 
 		return modified_points
+		

pygem/unvhandler.py

@@ -34,41 +34,26 @@ def parse(self, filename):
 
 		self.infile = filename
 
+		index = -9
+		mesh_points = []
 		with open(self.infile, 'r') as input_file:
-			nline = 0
-			while True:
-				line = input_file.readline()
-				nline += 1
-				if not line:
-					break
-				if line.startswith('    -1'):
-					section_id = input_file.readline().strip()
-					nline += 1
-					if section_id == '2411':
-						count = 0
-						while not input_file.readline().startswith('    -1'):
-							count += 1
-						start_line = nline + 2
-						last_line = start_line + count
+			for num, line in enumerate(input_file):
+				if line.startswith('  2411'):
+					index = num
+				if num == index + 2:
+					if line.startswith('    -1'):
+						break
 					else:
-						while not input_file.readline().startswith('    -1'):
-							nline += 1
-
-		n_points = count/2
-		mesh_points = np.zeros(shape=(n_points, 3))
-
-		nline = 0
-		i = 0
-		with open(self.infile, 'r') as input_file:
-			for line in input_file:
-				nline += 1
-				if nline % 2 == 1 and start_line < nline < last_line:
-					line = line.strip()
-					j = 0
-					for number in line.split():
-						mesh_points[i][j] = float(number)
-						j += 1
-					i += 1
+						line = line.replace('D', 'E')
+						l = []
+						for t in line.split():
+							try:
+								l.append(float(t))
+							except ValueError:
+								pass
+						mesh_points.append(l)
+						index = num
+			mesh_points = np.array(mesh_points)
 
 		return mesh_points
 
@@ -89,16 +74,23 @@ def write(self, mesh_points, filename):
 
 		self.outfile = filename
 
-		n_points = mesh_points.shape[0]
-		nrow = 0
+		index = -9
 		i = 0
-		with open(self.infile, 'r') as input_file, open(self.outfile, 'w') as output_file:
-			for line in input_file:
-				nrow += 1
-				if nrow % 2 == 1 and 20 < nrow <= (20 + n_points * 2):
-					for j in range(0, 3):
-						output_file.write('   ' + str(mesh_points[i][j]))
-					output_file.write('\n')
-					i += 1
-				elif nrow > 17:
-					output_file.write(line)
+		with open(self.outfile, 'w') as output_file:
+			with open(self.infile, 'r') as input_file:
+				for num, line in enumerate(input_file):
+					if line.startswith('  2411'):
+						index = num
+					if num == index + 2:
+						if line.startswith('    -1'):
+							index = -9
+							output_file.write(line)
+						else:
+							for j in range(0, 3):
+								output_file.write(3*' ' + '{:.16E}'.format(mesh_points[i][j]))
+							output_file.write('\n')
+							i += 1
+							index = num
+					else:
+						output_file.write(line)
+

tests/test_unvhandler.py

@@ -116,15 +116,9 @@ def test_unv_write_outfile(self):
 	def test_unv_write_comparison(self):
 		unv_handler = uh.UnvHandler()
 		mesh_points = unv_handler.parse('tests/test_datasets/test_square.unv')
-		mesh_points[0][0] = 2.2
-		mesh_points[5][1] = 4.3
-		mesh_points[9][2] = 0.5
-		mesh_points[45][0] = 7.2
-		mesh_points[132][1] = -1.2
-		mesh_points[255][2] = -3.6
 
 		outfilename = 'tests/test_datasets/test_square_out.unv'
-		outfilename_expected = 'tests/test_datasets/test_square_out_true.unv'
+		outfilename_expected = 'tests/test_datasets/test_square.unv'
 
 		unv_handler.write(mesh_points, outfilename)
 		self.assertTrue(filecmp.cmp(outfilename, outfilename_expected))