diff --git a/pygem/freeform.py b/pygem/freeform.py index 81d2655..45176fe 100644 --- a/pygem/freeform.py +++ b/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 + diff --git a/pygem/unvhandler.py b/pygem/unvhandler.py index 8262b4b..0f9af98 100644 --- a/pygem/unvhandler.py +++ b/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) + diff --git a/tests/test_unvhandler.py b/tests/test_unvhandler.py index 42f6b27..e6e99a3 100644 --- a/tests/test_unvhandler.py +++ b/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))