From a54366bcfd054d8f8ea3046bb71a0df67b6e60ad Mon Sep 17 00:00:00 2001 From: Mathew Topper Date: Mon, 12 Sep 2022 17:34:52 +0100 Subject: [PATCH] Catch some bugs when bad mesh files are used in dtocean-wec (#14) - Renamed `dtocean_wec.submodule.utils.mesh_util` module as `dtocean_wec.submodule.utils.mesh` and refactored - Fixed bugs that occur if opening wrong type of mesh file in dtocean-wec - Bump to version 3.2.0 --- CHANGELOG.md | 11 + README.md | 2 +- appveyor.yml | 4 +- dtocean_wec/submodule/nemoh_run.py | 75 ++-- dtocean_wec/submodule/utils/mesh.py | 460 ++++++++++++++++++++++ dtocean_wec/submodule/utils/mesh_util.py | 479 ----------------------- dtocean_wec/tab2.py | 8 +- dtocean_wec/utils/mesh_class.py | 89 ++--- setup.py | 4 +- test_data/cube.GDF | 9 + test_data/cube.dat | 28 ++ test_data/cube2.GDF | 24 ++ test_data/cube2.dat | 27 ++ tests/conftest.py | 5 + tests/test_wec_submodules_nemoh_run.py | 56 +++ tests/test_wec_submodules_utils_mesh.py | 392 +++++++++++++++++++ 16 files changed, 1096 insertions(+), 577 deletions(-) create mode 100644 dtocean_wec/submodule/utils/mesh.py delete mode 100644 dtocean_wec/submodule/utils/mesh_util.py create mode 100644 test_data/cube.GDF create mode 100644 test_data/cube.dat create mode 100644 test_data/cube2.GDF create mode 100644 test_data/cube2.dat create mode 100644 tests/test_wec_submodules_nemoh_run.py create mode 100644 tests/test_wec_submodules_utils_mesh.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 63c7e72..594c316 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,17 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](http://keepachangelog.com/) and this project adheres to [Semantic Versioning](http://semver.org/). +## [3.2.0] - 2022-09-12 + +### Changed + +- Renamed dtocean_wec.submodule.utils.mesh_util module as + dtocean_wec.submodule.utils.mesh and refactored. + +### Fixed + +- Fixed bugs that occur if opening wrong type of mesh file in dtocean-wec. + ## [3.1.1] - 2022-07-15 ### Fixed diff --git a/README.md b/README.md index cfeec92..4e4c144 100644 --- a/README.md +++ b/README.md @@ -129,7 +129,7 @@ $ conda activate _dtocean_hydro Install pytest to the environment (one time only): ``` -$ conda install -y mock pytest pytest-mock pytest-qt +$ conda install -y mock pytest pytest-catchlog pytest-mock pytest-qt ``` Run the tests: diff --git a/appveyor.yml b/appveyor.yml index 9992682..88f3d20 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -8,7 +8,7 @@ #---------------------------------# # version format -version: 3.1.1.build{build} +version: 3.2.0.build{build} image: - Visual Studio 2015 @@ -41,7 +41,7 @@ install: - python setup.py bootstrap - conda install --file requirements-conda-dev.txt - pip install -e . - - conda install mock "pytest>=3.6,<4" pytest-cov pytest-mock pytest-qt + - conda install mock "pytest>=3.6,<4" pytest-catchlog pytest-cov pytest-mock pytest-qt build: off diff --git a/dtocean_wec/submodule/nemoh_run.py b/dtocean_wec/submodule/nemoh_run.py index 5d5fb4b..08d75f0 100644 --- a/dtocean_wec/submodule/nemoh_run.py +++ b/dtocean_wec/submodule/nemoh_run.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- # Copyright (C) 2016 Pau Mercadez Ruiz -# Copyright (C) 2017-2018 Mathew Topper +# Copyright (C) 2017-2022 Mathew Topper # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -33,9 +33,10 @@ import subprocess import numpy as np +from numpy import linalg as LA from hydrostatics import Hydrostatics_Nemohcal -from utils.mesh_util import MeshBem +from utils.mesh import MeshBem from utils.multibody_analysis import MultiBodyAnalysis # Start logging @@ -198,29 +199,13 @@ def __init__(self, self.mesh_obj += [msh_obj] # identify the size of the circumscribing cylinder - if self.Cyl_Nzeta == 0 or self.Cyl_Ntheta == 0 or (not self.__get_array_mat): + if (self.Cyl_Nzeta == 0 or + self.Cyl_Ntheta == 0 or + not self.__get_array_mat): self.Cyl_radius = 0 - else : - coord = np.empty((0, 3)) - for b_mesh in self.mesh_obj: - coord = np.vstack((coord, b_mesh.Vertex)) - area = 0.0 - for panel in b_mesh.panels: - if panel.area > area: - area = panel.area - d = max([(panel.x[0]-panel.x[s])**2+(panel.y[0]-panel.y[s])**2+ - (panel.z[0]-panel.z[s])**2 for s in [1,-1]]) - - centroid = coord.mean(0) - if not np.sqrt(np.sum(centroid[:2]**2))/coord.max() < 1e-3: - module_logger.warning("WARNING: the centroid of the mesh file is not centered in the" - " origin of the mesh coordinate system. \nConsider to regenerate" - " a mesh that satisfy this condition.") - - self.Cyl_radius = (np.sqrt(coord[:,0]**2+coord[:,1]**2)).max() - self.Cyl_radius += np.sqrt(d) - - + else: + self.Cyl_radius = _get_cylinder_radius(self.mesh_obj) + def load_folder_tree(self): self.path_prj_hdy = os.path.join(self.prj_folder,'hydrodynamic') self.path_prj_dyn_res = os.path.join(self.path_prj_hdy,'results') @@ -254,14 +239,8 @@ def gen_mesh_files(self, show_norm=False): for nb in range(self.n_bodies): msh_obj = self.mesh_obj[nb] if self._debug: msh_obj.visualise_mesh() - - if self.bodies[nb]['mesh'][-3:].lower()=="gdf": - msh_obj.mesh_generation("nemoh", output_path=self.path_prj_dyn_mesh) - - else: - msh_obj.mesh_generation("nemoh", output_path=self.path_prj_dyn_mesh) - - if self._debug and self.n_bodies>1: msh_obj.visualise_mesh() + msh_obj.mesh_generation("nemoh", + output_path=self.path_prj_dyn_mesh) def gen_multibody_structure(self): mb_obj = MultiBodyAnalysis(self.bodies, self.shared_dof_binary, self.point_application) @@ -376,7 +355,39 @@ def run_nemoh(self, nemoh_folder, start=True): os.path.join(nemoh_folder, 'postprocessor.exe') ) os.chdir(actual_dir) + + +def _get_cylinder_radius(meshes): + + coord = np.empty((0, 3)) + d = None + + for mesh in meshes: + coord = np.vstack((coord, mesh.vertices)) + + for panel in mesh.connectivity: + + d_panel = max([LA.norm(mesh.vertices[panel[0]] - + mesh.vertices[panel[s]]) for s in [1, -1]]) + + if d is None or d_panel > d: + d = d_panel + + if d is None: + raise RuntimeError("Can not determine mesh extents") + + centroid = coord.mean(0) + + if not np.sqrt(np.sum(centroid[:2] ** 2)) / coord.max() < 1e-3: + module_logger.warning( + "WARNING: the centroid of the mesh file is not centered " + "at the origin of the mesh coordinate system.\n" + "Consider regenerating a mesh to satisfy this condition.") + + return (np.sqrt(coord[:, 0] ** 2 + coord[:, 1] ** 2)).max() + d + + def rot_matrix(ang): Rz = lambda angle: np.array([[np.cos(angle), -np.sin(angle), 0], [np.sin(angle), np.cos(angle), 0], [0, 0, 1]]) Ry = lambda angle: np.array([[np.cos(angle), 0, np.sin(angle)], [0, 1, 0], [-np.sin(angle), 0, np.cos(angle)]]) diff --git a/dtocean_wec/submodule/utils/mesh.py b/dtocean_wec/submodule/utils/mesh.py new file mode 100644 index 0000000..8ccf33a --- /dev/null +++ b/dtocean_wec/submodule/utils/mesh.py @@ -0,0 +1,460 @@ +# -*- coding: utf-8 -*- + +# Copyright (C) 2016 Francesco Ferri +# Copyright (C) 2022 Mathew Topper +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import os +import re + +import numpy as np +import matplotlib.pyplot as plt +from numpy import linalg as LA +from mpl_toolkits.mplot3d import axes3d # pylint: disable=unused-import + + +class MeshBem(object): + """ + Mesh_BEM: class used to open, visualise, transform and save structured + meshes + + Args: + file_name (str): file name of the mesh to be read + + Optional args: + path (str): location of the mesh file + + Attributes: + file_name (str): file name of the mesh to be read + path (str): location of the mesh file + xsim (int): index identifying the symmetry around the x-axis + connectivity (numpy.ndarray): vertex ID of each panel, listed in a + counter-clockwise direction from the fluid perspective + vertices (numpy.ndarray): x,y,z coordinates of the mesh vertex + nP (int): number of panel of the mesh + nV (int): number of vertex of the mesh + """ + + def __init__(self, file_name, path=""): + + self.file_name = file_name + self.path = path + + mesh_fn = os.path.join(path, file_name) + _, extension = os.path.splitext(mesh_fn) + + if extension.lower() == ".gdf": + xsim, vertices, connectivity = read_WAMIT(mesh_fn) + elif extension.lower() == ".dat": + xsim, vertices, connectivity = read_NEMOH(mesh_fn) + else: + raise IOError("Mesh file type not supported. Use GDF or dat.") + + self.xsim = xsim + self.vertices = vertices + self.connectivity = np.int_(connectivity) + self.nP = len(connectivity) + self.nV = len(vertices) + + def translate(self, x, y, z): + """ + translate: translates the mesh to the given point + + Args: + x (float) [m]: x-axis translation + y (float) [m]: y-axis translation + z (float) [m]: z-axis translation + """ + self.vertices += np.array([x,y,z]) + + def rotate(self, rotZ, pivot): + """ + rotate: rotates the mesh around the given pivoting point + + Args: + rotZ (float)[rad]: rotation angle + pivot (numpy.ndarray) [m]: pivoting point coordinates + """ + + v = self.vertices.copy() + v -= pivot + + R = np.array([[np.cos(rotZ), -np.sin(rotZ), 0], + [np.sin(rotZ), np.cos(rotZ), 0], + [0, 0, 1]]) + + self.vertices = np.dot(R,v.T).T + pivot + + def invert_norm(self, index=None): + """ + invert_norm: used to invert the direction of the specified panel norm + + Args: + index: index of the panel subject to the transformation. If None + the panels will be transformed + """ + + if index is None: + self.connectivity = np.fliplr(self.connectivity) + return + + self.connectivity[index, :] = self.connectivity[index, ::-1] + + def visualise_mesh(self, ax=None, fig=None): + """ + visualise_mesh: plot the mesh grid + + Optional args: + ax (matplotlib axes): parent figure axes pointer + fig (matplotlib figure): parent figure pointer + """ + + panels = [Panel(self.vertices[panel, :]) + for panel in self.connectivity] + + if ax is None: + fig, ax = plt.subplots(2, 2) + ax[1, 1] = fig.add_subplot(224, projection="3d") + + ax[0, 0].margins(0.05) + + for elm in panels: + elm.show(ax[0, 0], dimension=1) + + ax[0, 0].set_aspect('equal') + plt.grid() + + ax[0, 1].margins(0.05) + + for elm in panels: + elm.show(ax[0, 1], dimension=2) + + ax[0, 1].set_aspect('equal') + plt.grid() + + ax[1, 0].margins(0.05) + + for elm in panels: + elm.show(ax[1, 0], dimension=3) + + ax[1, 0].set_aspect('equal') + plt.grid() + + for elm in panels: + elm.show(ax[1, 1]) + + ax[1, 1].set_aspect('equal') + plt.show() + + def visualise_norm(self, scale=1): + """ + visualise_mesh: plot the mesh grid + + Optional args: + scale (float): scaling factor for the norm. Only for visualisation + """ + + panels = [Panel(self.vertices[panel, :]) + for panel in self.connectivity] + + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + + for elm in panels: + elm.show_norm(ax) + + ax.set_aspect('equal') + plt.show() + + def mesh_generation(self, output_format, output_path=None): + """ + mesh_generation: generate the mesh file specified in the output_format + + Args: + output_format (str): defines the output format of the saved files. + gdf: wamit standard + nemoh: nemoh dat file + mesh: nenoh dat file for hydrostatic calculation + + Optional args: + output_path (str): name of the folder where to save the file + """ + + if output_path is None: output_path = self.path + + if output_format == "gdf": + + file_n = '{}gdf'.format(self.file_name[0:-3]) + file_p = os.path.join(output_path, file_n) + + with open(file_p, 'w') as f: + + f.write('WAMIT mesh file.\n') + f.write('1 9.80665 ULEN GRAV\n') + f.write('{} 0 ISX ISY\n'.format(self.xsim)) + f.write('{}\n'.format(self.nP)) + + for panel in self.connectivity: + + for vertex in panel[:-1]: + f.write('{:15.6f} {:15.6f} {:15.6f}'.format( + self.vertices[vertex, 0], + self.vertices[vertex, 1], + self.vertices[vertex, 2])) + f.write(' ') + + f.write('{:15.6f} {:15.6f} {:15.6f}'.format( + self.vertices[panel[-1], 0], + self.vertices[panel[-1], 1], + self.vertices[panel[-1], 2])) + f.write('\n') + + elif output_format == "nemoh": + + file_n = '{}dat'.format(self.file_name[0:-3]) + file_p = os.path.join(output_path, file_n) + float_format = '{:>10d} {:15.6f} {:15.6f} {:15.6f}\n' + int_format = '{:>10d} {:>10d} {:>10d} {:>10d}\n' + + with open(file_p, 'w') as f: + + f.write('2 0\n') + + for vertex in range(self.nV): + f.write(float_format.format(vertex + 1, + self.vertices[vertex, 0], + self.vertices[vertex, 1], + self.vertices[vertex, 2])) + + f.write(float_format.format(0, 0, 0, 0)) + + for panel in range(self.nP): + f.write(int_format.format( + *(self.connectivity[panel,:] + 1))) + + f.write(int_format.format(0, 0, 0, 0)) + + elif output_format == "mesh": + + file_n = '{}dat'.format(self.file_name[0:-3]) + file_p = os.path.join(output_path, file_n) + float_format = '{:15.6f} {:15.6f} {:15.6f}\n' + int_format = '{:>10d} {:>10d} {:>10d} {:>10d}\n' + + with open(file_p, 'w') as f: + + f.write('{}\n'.format(self.nV)) + f.write('{}\n'.format(self.nP)) + + for vertex in range(self.nV): + f.write(float_format.format(self.vertices[vertex, 0], + self.vertices[vertex, 1], + self.vertices[vertex, 2])) + + for panel in range(self.nP): + f.write(int_format.format( + *(self.connectivity[panel,:] + 1))) + + else: + + raise ValueError("Unsupported output format. Use gdf, nemoh or " + "mesh") + + +class Panel(object): + """ + Panel: a class for visualizing a 4-node BEM panel + + Args: + vertex (numpy.ndarray): coordinates of the vertexs + + Attributes: + x (numpy.ndarray): x coordinates of the vertexs + y (numpy.ndarray): y coordinates of the vertexs + z (numpy.ndarray): z coordinates of the vertexs + centroid (numpy.ndarray): panel centroid + n (numpy.ndarray): normal vector at the panel centroid + """ + def __init__(self, vertex): + + self.x = vertex[[0, 1, 2, 3], 0] + self.y = vertex[[0 ,1 ,2, 3], 1] + self.z = vertex[[0, 1, 2, 3], 2] + self.centroid = np.mean(vertex, axis=0) + self.n = _get_panel_norm(self.x, self.y, self.z) + + def show_norm(self, ax): + """ + show_norm: plot the norm or the panel into the given axes + + Args: + ax (matplotlib axes): axes of the parent figure pointer + """ + + x = self.x + y = self.y + z = self.z + + v24 = np.array([x[3] - x[1], + y[3] - y[1], + z[3] - z[1]]) + v31 = np.array([x[0] - x[2], + y[0] - y[2], + z[0] - z[2]]) + + scale = (LA.norm(v24) + LA.norm(v31)) / 8 + d = self.n * scale + self.centroid + + ax.plot_wireframe(x, y, z, color="#000000") + ax.plot_wireframe([self.centroid[0], d[0]], + [self.centroid[1], d[1]], + [self.centroid[2], d[2]], + color="red") + ax.scatter(self.centroid[0], + self.centroid[1], + self.centroid[2], + c="r", + marker="o") + + def show(self, ax, dimension=None): + """ + show: plot the panel into the given axes + + Args: + ax (matplotlib axes): axes of the parent figure pointer + + Optional args: + dimension (int): + """ + + if not dimension is None: + + V = [self.x,self.y,self.z] + Vred = [el for ind, el in enumerate(V) if ind != dimension - 1] + ax.plot(Vred[0], Vred[1], 'k') + + else: + + ax.plot_wireframe(self.x[[0, 1, 2, 3, 0]], + self.y[[0, 1, 2, 3, 0]], + self.z[[0, 1, 2, 3, 0]]) + + +def _get_panel_norm(x, y, z): + + v24 = np.array([x[3] - x[1], + y[3] - y[1], + z[3] - z[1]]) + v31 = np.array([x[0] - x[2], + y[0] - y[2], + z[0] - z[2]]) + + dr = np.cross(v24, v31) + + return dr / LA.norm(dr) + + +def read_NEMOH(f_n): + + with open(f_n,'r') as mesh_f: + msg = mesh_f.read() + + msg = strip_comments(msg) + lines = msg.split('\n') + + first_line = np.array(lines.pop(0).split(), dtype=int) + + if len(first_line) == 1: + + xsim = 0 + nV = first_line[0] + nP = int(lines.pop(0)) + vertices = np.zeros((nV, 3)) + + for vertex in range(nV): + vertices[vertex, :] = np.array(lines.pop(0).split(), dtype=float) + + connectivity = np.zeros((nP, 4)) + + for panel in range(nP): + connectivity[panel, :] = np.array(lines.pop(0).split(), + dtype=float) - 1 + + elif len(first_line) == 2: + + xsim = first_line[1] + vertices = [] + connectivity = [] + pass_to_connectivity = False + + for line in lines: + + if not line: continue + + temp = np.array(line.split(), dtype=float) + + if not int(temp[0]) == 0 and not pass_to_connectivity: + vertices.append(temp[1:].tolist()) + else: + pass_to_connectivity = True + connectivity.append([v - 1 for v in temp.tolist()]) + + vertices = np.array(vertices) + connectivity.pop(0) + connectivity.pop(-1) + + connectivity = np.array(connectivity, dtype=int) + + return xsim, vertices, connectivity + + +def read_WAMIT(f_n): + + with open(f_n,'r') as mesh_f: + msg = mesh_f.read() + + msg = strip_comments(msg) + lines = msg.split('\n') + lines = lines[2:] + + xsim = np.array(lines.pop(0).split()[0], dtype=int) + nP = np.array(lines.pop(0), dtype=int) + nV = nP * 4 + + points = np.array([float(val) for entries in lines + for val in entries.split()]) + vertices = np.reshape(points, (-1, 3)) + + assert vertices.shape[0] == nV + + connectivity = np.zeros((nP, 4)) + vertex = 0 + + for panel in range(nP): + connectivity[panel, :] = np.array([vertex, + vertex + 1, + vertex + 2, + vertex + 3]) + vertex +=4 + + connectivity = np.array(connectivity, dtype=int) + + return xsim, vertices, connectivity + + +def strip_comments(code): + msg = re.sub('\s*!.*', '', code, re.MULTILINE) + msg = re.sub(r'^\n', '', msg) + return re.sub(r'\n\s*\n', '\n', msg, re.MULTILINE) diff --git a/dtocean_wec/submodule/utils/mesh_util.py b/dtocean_wec/submodule/utils/mesh_util.py deleted file mode 100644 index 04d9474..0000000 --- a/dtocean_wec/submodule/utils/mesh_util.py +++ /dev/null @@ -1,479 +0,0 @@ -# -*- coding: utf-8 -*- - -# Copyright (C) 2016 Francesco Ferri -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - -""" -This module contains the classes used read and visualise the mesh of the WEC structure - -.. module:: input - :platform: Windows - :synopsis: Input module to DTOcean WP2 - -.. moduleauthor:: Francesco Ferri -""" -import os -import matplotlib.animation as animation -import matplotlib.pyplot as plt -from numpy import array, zeros, cos, int_, cross, mean, dot -from mpl_toolkits.mplot3d import axes3d -from numpy import linalg as LA -import numpy as np - - -class MeshBem(): - """ - Mesh_BEM: class used to open, visualise, transform and save structured meshes - - Args: - file_name (str): file name of the mesh to be read - - Optional args: - path (str): location of the mesh file - - Attributes: - file_name (str): file name of the mesh to be read - path (str): location of the mesh file - mesh_fn (str): full path name to the mesh file - xsim (int): index identifying the symmetry around the x-axis - Connectivity (numpy.ndarray): vertex ID of each panel, listed in a counter-clockwise direction from the fluid perspective - Vertex (numpy.ndarray): x,y,z coordinates of the mesh vertex - panels (list): list of panel objects - nP (int): number of panel of the mesh - nV (int): number of vertex of the mesh - """ - def __init__(self, file_name, path=""): - self.file_name = file_name - self.path = path - self.mesh_fn = os.path.join(path,file_name) - if self.mesh_fn[-3:]=="gdf" or self.mesh_fn[-3:]=="GDF": - with open(self.mesh_fn,'r') as mesh_f: - burnheader= mesh_f.readline(); del(burnheader) - burnheader= mesh_f.readline(); del(burnheader) - self.xsim= array(mesh_f.readline().split()[0],dtype=int) - nP = array(mesh_f.readline().partition('!')[0],dtype=int) - Vertex = zeros((nP*4,3)) - for vertex in range(nP*4): - Vertex[vertex,:] = array(mesh_f.readline().partition('!')[0].split(),dtype=float) - Connectivity = zeros((nP,4)) - vertex = 1 - for panel in range(nP): - Connectivity[panel,:] = array([vertex,vertex+1,vertex+2,vertex+3]) - vertex +=4 - elif self.mesh_fn[-3:]=="dat": - f = open(self.mesh_fn,'r') - first_line = array(f.readline().split(),dtype = int) - if len(first_line) == 1: - nV = first_line - nP = array(f.readline().split(),dtype = int) - Vertex = zeros((nV,3)) - for vertex in range(nV): - Vertex[vertex,:] = array(f.readline().split(),dtype=float) - Connectivity = zeros((nP,4)) - for panel in range(nP): - Connectivity[panel,:] = array(f.readline().split(),dtype=float) - elif len(first_line) == 2: - self.xsim = first_line[1] - Vertex = [] - Connectivity = [] - pass_to_connectivity = False - for line in f: - temp = array(line.split(),dtype = float) - if not int(temp[0]) == 0 and not pass_to_connectivity: - Vertex.append(temp[1:].tolist()) - else: - pass_to_connectivity = True - Connectivity.append(temp.tolist()) - - Vertex = array(Vertex) - Connectivity.pop(0) - Connectivity.pop(-1) - Connectivity = array(Connectivity,dtype=int) - - else: - pass - else: - pass - self.Connectivity = int_(Connectivity) - self.Vertex = Vertex - self.nP = len(Connectivity) - self.nV = len(Vertex) - - self.panels = [ Panel(self.Vertex[panel-1,:]) for panel in self.Connectivity] - - def translate(self, x, y, z): - """ - translate: translates the mesh to the given point - - Args: - x (float) [m]: x-axis translation - y (float) [m]: y-axis translation - z (float) [m]: z-axis translation - """ - self.Vertex += np.array([x,y,z]) - self.panels = [ Panel(self.Vertex[panel-1,:]) for panel in self.Connectivity] - - #for pn in self.panels: - # pn.translate(delt=np.array([x, y, z])) - - def rotate(self, rotZ, pivot): - """ - rotate: rotates the mesh around the given pivoting point - - Args: - rotZ (float)[rad]: rotation angle - pivot (numpy.ndarray) [m]: pivoting point coordinates - """ - v = self.Vertex.copy() - v -= pivot - - R = np.array([[np.cos(rotZ), -np.sin(rotZ), 0], - [np.sin(rotZ), np.cos(rotZ), 0], - [0, 0, 1]]) - - self.Vertex = np.dot(R,v.T).T+pivot - self.panels = [ Panel(self.Vertex[panel-1,:]) for panel in self.Connectivity] - #for pn in self.panels: - # pn.rotate(rotZ, pivot) - - - def invertNorm(self,index=-1): - """ - invertNorm: used to invert the direction of the specified panel norm - - Args: - index: index of the panel subject to the transformation. If -1 all the panels will be transformed - """ - if index==-1: - index = range(self.nP) - - for pn in index: - self.panels[pn].invert_direction() - - Vertex = zeros((self.nP*4,3)) - ind_v = -1 - for _p in self.panels: - for _v in range(4): - ind_v +=1 - Vertex[ind_v,:] = array([_p.x[_v],_p.y[_v],_p.z[_v]],dtype=float) - - Connectivity = zeros((self.nP,4)) - vertex = 1 - for panel in range(self.nP): - Connectivity[panel,:] = array([vertex,vertex+1,vertex+2,vertex+3]) - vertex +=4 - - self.Connectivity = int_(Connectivity) - self.Vertex = array(Vertex) - - def visualise_mesh(self, a=None, f=None): - """ - visualise_mesh: plot the mesh grid - - Optional args: - a (matplotlib axes): parent figure axes pointer - f (matplotlib figure): parent figure pointer - """ - if a is None: - f, a = plt.subplots(2, 2) - a[1,1] = f.add_subplot(224, projection="3d") - - a[0,0].margins(0.05) - for elm in self.panels: - elm.show(a[0,0],dimension=1) - - a[0,0].set_aspect('equal') - plt.grid() - - a[0,1].margins(0.05) - for elm in self.panels: - elm.show(a[0,1],dimension=2) - - a[0,1].set_aspect('equal') - plt.grid() - - a[1,0].margins(0.05) - for elm in self.panels: - elm.show(a[1,0],dimension=3) - - a[1,0].set_aspect('equal') - plt.grid() - - - for elm in self.panels: - elm.show(a[1,1]) - - a[1,1].set_aspect('equal') - plt.show() - - def visualise_norm(self, scale=1): - """ - visualise_mesh: plot the mesh grid - - Optional args: - scale (float): scaling factor for the norm. Only for visualisation - """ - f = plt.figure() - a = f.add_subplot(111, projection='3d') - for elm in self.panels: - #elm.show(a) - elm.norm(scale) - elm.show_norm(a) - a.set_aspect('equal') - plt.show() - - def __animate(self,i,a): - """ - Unused too slow - :param i: - :param a: - :return: - """ - delta = [0.1*cos(i/10.*3.14),0,0] - a.clear() - a.set_aspect('equal') - a.set_xlim3d(-1, 1) - a.set_ylim3d(-1,1) - a.set_zlim3d(-1,1) - - for elm in self.panels: - elm.update(delta) - elm.show(a) - - def animate_mesh(self): - """ - Unused too slow - :return: - """ - f = plt.figure() - a = f.add_subplot(111, projection='3d') - ani = animation.FuncAnimation(f, self.__animate , fargs=[a], repeat = "none", frames = 2 ,interval=10) - plt.show() - plt.close(f) - - def mesh_generation(self, output_format, output_path=""): - """ - mesh_generation: generate the mesh file specified in the output_format - - Args: - output_format (str): defines the output format of the saved files. - gdf: wamit standard - nemoh: nemoh dat file - mesh: nenoh dat file for hydrostatic calculation - - Optional args: - output_path (str): name of the location where to save the file - """ - if output_path=="": - output_path = self.path - if output_format == "gdf": - file_n = '{}gdf'.format(self.file_name[0:-3]) - f = open(os.path.join(output_path,'{}gdf'.format(self.file_name[0:-3])), 'w') - f.write('WAMIT mesh file.\n') - f.write('1 9.82 ULEN GRAV\n') - f.write('{} 0 ISX ISY\n'.format(self.xsim)) - f.write('{}\n'.format(self.nP)) - for vertex in range(self.nV): - f.write('{} {} {}\n'.format(self.Vertex[vertex,0],self.Vertex[vertex,1],self.Vertex[vertex,2])) - - f.close() - elif output_format == "nemoh": - file_n = '{}dat'.format(self.file_name[0:-3]) - f = open(os.path.join(output_path,'{}dat'.format(self.file_name[0:-3])), 'w') - f.write('2 0\n') - for vertex in range(self.nV): - f.write('{} {} {} {}\n'.format(vertex+1,self.Vertex[vertex,0],self.Vertex[vertex,1],self.Vertex[vertex,2])) - - f.write('0 0.00 0.00 0.00\n') - for panel in range(self.nP): - f.write('{} {} {} {}\n'.format(self.Connectivity[panel,0],self.Connectivity[panel,1],self.Connectivity[panel,2],self.Connectivity[panel,3])) - f.write('0 0 0 0\n') - f.close() - elif output_format == "mesh": - file_n = '{}_mesh.dat'.format(self.file_name[0:-4]) - f = open(os.path.join(output_path,'{}_mesh.dat'.format(self.file_name[0:-4])), 'w') - f.write('{}\n'.format(self.nV)) - f.write('{}\n'.format(self.nP)) - for vertex in range(self.nV): - f.write('{} {} {}\n'.format(self.Vertex[vertex,0],self.Vertex[vertex,1],self.Vertex[vertex,2])) - - for panel in range(self.nP): - f.write('{} {} {} {}\n'.format(self.Connectivity[panel,0],self.Connectivity[panel,1],self.Connectivity[panel,2],self.Connectivity[panel,3])) - f.close() - - else: - #print "Error: wrong output_format!" - pass - - #print("File {} saved in {}".format(file_n,output_path)) - - - -class Panel(): - """ - Panel: the class is the base object of the mesh class. Each panel is defined by the four vertexs - - Args: - vertex (numpy.ndarray): coordinates of the vertexs - - Attributes: - x (numpy.ndarray): x coordinates of the vertexs - y (numpy.ndarray): y coordinates of the vertexs - z (numpy.ndarray): z coordinates of the vertexs - centroid (numpy.ndarray): panel centroid - area (float): panel area - dr (numpy.ndarray): normal direction - d (numpy.ndarray): normal direction applied at the panel centroid - """ - def __init__(self,vertex): - self.x = vertex[[0,1,2,3],0] - self.y = vertex[[0,1,2,3],1] - self.z = vertex[[0,1,2,3],2] - self.centroid = mean(vertex,axis=0) - self.norm() - - def invert_direction(self): - """ - invert_direction: invert the direction of the panel vertex, inverting the normal direction - - """ - self.x = self.x[-1::-1] - self.y = self.y[-1::-1] - self.z = self.z[-1::-1] - - self.norm() - - def norm(self, scale=1): - """ - norm: calculated the panel norm - - Optional args: - scale: scaling factor of the normal direction - """ - x = self.x - y = self.y - z = self.z - - v12 = array([x[1]-x[0],y[1]-y[0],z[1]-z[0]]) - v14 = array([x[3]-x[0],y[3]-y[0],z[3]-z[0]]) - v34 = array([x[3]-x[2],y[3]-y[2],z[3]-z[2]]) - v32 = array([x[1]-x[2],y[1]-y[2],z[1]-z[2]]) - - - n1 = cross(v12,v14) - n2 = cross(v34,v32) - d1 = n1/LA.norm(n1) - d2 = n2/LA.norm(n2) - - self.area = (LA.norm(n1)+dot(d1,d2)*LA.norm(n2))/2. - - self.dr = (((d1+d2)/2.*self.area))*scale - self.d = self.dr+self.centroid - - def show_norm(self,a): - """ - show_norm: plot the norm or the panel into the given axes - - Args: - a (matplotlib axes): axes of the parent figure pointer - """ - a.plot_wireframe(self.x,self.y,self.z,color="#000000") - a.plot_wireframe([self.centroid[0],self.d[0]],[self.centroid[1],self.d[1]],[self.centroid[2],self.d[2]],color="red") - a.scatter(self.centroid[0],self.centroid[1],self.centroid[2],c="r", marker="o") - #a.scatter(self.d[0],self.d[1],self.d[2],c="g", marker=">") - - def show(self,a,dimension=None): - """ - show: plot the panel into the given axes - - Args: - a (matplotlib axes): axes of the parent figure pointer - - Optional args: - dimension (int): - """ - if not dimension is None: - V = [self.x,self.y,self.z] - Vred = [el for ind,el in enumerate(V) if not ind==dimension-1] - a.plot(Vred[0],Vred[1],'k') - else: - a.plot_wireframe(self.x[[0,1,2,3,0]],self.y[[0,1,2,3,0]],self.z[[0,1,2,3,0]]) - - def translate(self, delt): - """ - translate: translates the panel by the amount specified in delt - - Args: - delt (numpy.ndarray): transformation vector - """ - self.x = self.x+delt[0] - self.y = self.y+delt[1] - self.z = self.z+delt[2] - self.centroid = array([mean(self.x), - mean(self.y), - mean(self.z)]) - self.norm(scale=1) - - def rotate(self,rotZ, pivot): - """ - rotate: rotates the panel by the amount specified in rotZ, arount the pivoting point - - Args: - rotZ (float) [rad]: angle defining the transformation - pivot (numpy.ndarray): pivoting point coordinates - """ - vert = np.vstack((self.x, self.y, self.z)).T - R = np.array([[np.cos(rotZ), -np.sin(rotZ), 0], - [np.sin(rotZ), np.cos(rotZ), 0], - [0, 0, 1]]) - vert_r = np.dot(R,(vert-pivot).T).T+pivot - - self.x = vert_r[:,0] - self.y = vert_r[:,1] - self.centroid = array([mean(self.x), - mean(self.y), - mean(self.z)]) - self.norm(scale=1) - - def update(self,delt): # not used now - """ - unused - :param delt: - :return: - """ - self.x = self.x+delt[0] - self.y = self.y+delt[1] - self.z = self.z+delt[2] - - -if __name__== "__main__": - - path = r"C:\Users\francesco\Desktop\Nemoh_run" - name = "Cylinder.GDF" - - m = MeshBem(name, path=path) - m.visualise_mesh() - m.translate(5,0,0) - m.rotate(90./180*np.pi, np.array([0,0,0])) - m.visualise_mesh() - m.rotate(-90./180*np.pi, np.array([0,0,0])) - m.rotate(90./180*np.pi, m.Vertex.mean(0)) - m.visualise_mesh() - - -#path = "C:\\Users\\fferri\\Documents\\DTOcean\\Nemoh_module\\version0\\nemoh\\Cylsurface\\Mesh" -#mesh = "Flap.dat" -#mesh_obj = ConvertMesh(mesh,path) -#mesh_obj.mesh_generation("gdf",output_path="..\\..\\..\\..\\PMfitting") diff --git a/dtocean_wec/tab2.py b/dtocean_wec/tab2.py index dd987d9..0f4bd96 100644 --- a/dtocean_wec/tab2.py +++ b/dtocean_wec/tab2.py @@ -150,9 +150,13 @@ def load_nemoh_solution(): """ """ - + def browse_folder(self): - folder = QFileDialog.getOpenFileName(self, "Select the mesh file of the body") + folder = QFileDialog.getOpenFileName( + self, + "Select the mesh file of the body", + "", + "mesh files (*.dat *.gdf *.GDF)") self.mesh_f_t2.setText(folder) def submit_input_data(self): diff --git a/dtocean_wec/utils/mesh_class.py b/dtocean_wec/utils/mesh_class.py index 8648c7b..afdeaf3 100644 --- a/dtocean_wec/utils/mesh_class.py +++ b/dtocean_wec/utils/mesh_class.py @@ -1,8 +1,30 @@ -import numpy as np -from numpy import array, zeros +# -*- coding: utf-8 -*- + +# Copyright (C) 2016 Francesco Ferri +# Copyright (C) 2022 Mathew Topper +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + import os +import numpy as np + +from ..submodule.utils.mesh import read_NEMOH, read_WAMIT + + class Mesh(): + def __init__(self, path, fn): self.fn = os.path.join(path, fn) self.path = path @@ -11,11 +33,10 @@ def __init__(self, path, fn): mesh_type = self.fn[-3:] if mesh_type.lower() == "gdf": - v, p = self.readWAMIT(self.fn) - p -= 1 + _, v, p = read_WAMIT(self.fn) elif mesh_type.lower() == "dat": # try to see if the specified filename is Nemoh format compliant - v, p = self.readNemoh(self.fn) + _, v, p = read_NEMOH(self.fn) else: raise IOError('Mesh type not understood') @@ -24,7 +45,7 @@ def __init__(self, path, fn): n, c = self.__evalnorm() self.n = n self.c = c - + def __evalnorm(self): pan = self.p+1 ver = self.v @@ -78,58 +99,6 @@ def __evalnorm(self): # n /= np.linalg.norm(n) # print(np.linalg.norm(n)) return n, centroid - - def readNemoh(self, f_n): - with open(f_n,'r') as mesh_f: - first_line = array(mesh_f.readline().split(),dtype = int) - if len(first_line) == 1: - nV = first_line - nP = array(mesh_f.readline().split(),dtype = int) - Vertex = zeros((nV,3)) - for vertex in range(nV): - Vertex[vertex,:] = array(mesh_f.readline().split(),dtype=float) - Connectivity = zeros((nP,4)) - for panel in range(nP): - Connectivity[panel,:] = array(mesh_f.readline().split(),dtype=float) - elif len(first_line) == 2: - self.xsim = first_line[1] - Vertex = [] - Connectivity = [] - pass_to_connectivity = False - for line in mesh_f: - temp = array(line.split(),dtype = float) - if not int(temp[0]) == 0 and not pass_to_connectivity: - Vertex.append(temp[1:].tolist()) - else: - pass_to_connectivity = True - Connectivity.append(temp.tolist()) - - Vertex = array(Vertex) - Connectivity.pop(0) - Connectivity.pop(-1) - Connectivity = array(Connectivity,dtype=int) - else: - raise IOError('Mesh type not understood') - - return Vertex, Connectivity - - def readWAMIT(self, f_n): - - with open(f_n,'r') as mesh_f: - burnheader= mesh_f.readline(); del(burnheader) - burnheader= mesh_f.readline(); del(burnheader) - xsim= array(mesh_f.readline().split()[0],dtype=int) - nP = array(mesh_f.readline().partition('!')[0],dtype=int) - Vertex = zeros((nP*4,3)) - for vertex in range(nP*4): - Vertex[vertex,:] = array(mesh_f.readline().partition('!')[0].split(),dtype=float) - Connectivity = zeros((nP,4),dtype=int) - vertex = 1 - for panel in range(nP): - Connectivity[panel,:] = array([vertex,vertex+1,vertex+2,vertex+3]) - vertex +=4 - - return Vertex, Connectivity def gen_vtpxml(self): points = self.v @@ -202,10 +171,10 @@ def gen_vtpxml(self): print("vtk xml file (__Mesh.vtp) generated.") return 1 - + + if __name__ == "__main__": fn = "mesh_test.gdf" path = "" ms = Mesh(fn, path) ms.generate_visualiser() - diff --git a/setup.py b/setup.py index 261cb7b..c2362f1 100644 --- a/setup.py +++ b/setup.py @@ -256,7 +256,9 @@ def get_appveyor_version(): zip_safe=False, # Important for reading config files tests_require=['mock', 'pytest', - 'pytest-mock'], + 'pytest-catchlog', + 'pytest-mock', + 'pytest-qt'], cmdclass={'bootstrap': Bootstrap, 'test': PyTest, 'cleantest': CleanTest} diff --git a/test_data/cube.GDF b/test_data/cube.GDF new file mode 100644 index 0000000..5862252 --- /dev/null +++ b/test_data/cube.GDF @@ -0,0 +1,9 @@ +Cube +1 9.80665 ULEN GRAV +0 0 ISX ISY +5 + 1.0 -1.0 0.0 -1.0 -1.0 0.0 -1.0 -1.0 -2.0 1.0 -1.0 -2.0 + -1.0 -1.0 0.0 -1.0 1.0 0.0 -1.0 1.0 -2.0 -1.0 -1.0 -2.0 + -1.0 1.0 0.0 1.0 1.0 0.0 1.0 1.0 -2.0 -1.0 1.0 -2.0 + 1.0 1.0 0.0 1.0 -1.0 0.0 1.0 -1.0 -2.0 1.0 1.0 -2.0 + 1.0 1.0 -2.0 1.0 -1.0 -2.0 -1.0 -1.0 -2.0 -1.0 1.0 -2.0 diff --git a/test_data/cube.dat b/test_data/cube.dat new file mode 100644 index 0000000..b949b1a --- /dev/null +++ b/test_data/cube.dat @@ -0,0 +1,28 @@ + 2 0 + 1 1.000000 -1.000000 0.000000 + 2 -1.000000 -1.000000 0.000000 + 3 -1.000000 -1.000000 -2.000000 + 4 1.000000 -1.000000 -2.000000 + 5 -1.000000 -1.000000 0.000000 + 6 -1.000000 1.000000 0.000000 + 7 -1.000000 1.000000 -2.000000 + 8 -1.000000 -1.000000 -2.000000 + 9 -1.000000 1.000000 0.000000 + 10 1.000000 1.000000 0.000000 + 11 1.000000 1.000000 -2.000000 + 12 -1.000000 1.000000 -2.000000 + 13 1.000000 1.000000 0.000000 + 14 1.000000 -1.000000 0.000000 + 15 1.000000 -1.000000 -2.000000 + 16 1.000000 1.000000 -2.000000 + 17 1.000000 1.000000 -2.000000 + 18 1.000000 -1.000000 -2.000000 + 19 -1.000000 -1.000000 -2.000000 + 20 -1.000000 1.000000 -2.000000 + 0 0 0 0 + 2 3 4 1 + 6 7 8 5 + 10 11 12 9 + 14 15 16 13 + 18 19 20 17 + 0 0 0 0 diff --git a/test_data/cube2.GDF b/test_data/cube2.GDF new file mode 100644 index 0000000..d4fc9a7 --- /dev/null +++ b/test_data/cube2.GDF @@ -0,0 +1,24 @@ +Cube +1 9.80665 ULEN GRAV +0 0 ISX ISY +5 + 1.0 -1.0 0.0 + -1.0 -1.0 0.0 + -1.0 -1.0 -2.0 + 1.0 -1.0 -2.0 + -1.0 -1.0 0.0 + -1.0 1.0 0.0 + -1.0 1.0 -2.0 + -1.0 -1.0 -2.0 + -1.0 1.0 0.0 + 1.0 1.0 0.0 + 1.0 1.0 -2.0 + -1.0 1.0 -2.0 + 1.0 1.0 0.0 + 1.0 -1.0 0.0 + 1.0 -1.0 -2.0 + 1.0 1.0 -2.0 + 1.0 1.0 -2.0 + 1.0 -1.0 -2.0 + -1.0 -1.0 -2.0 + -1.0 1.0 -2.0 diff --git a/test_data/cube2.dat b/test_data/cube2.dat new file mode 100644 index 0000000..31ebe13 --- /dev/null +++ b/test_data/cube2.dat @@ -0,0 +1,27 @@ + 20 + 5 + 1.000000 -1.000000 0.000000 + -1.000000 -1.000000 0.000000 + -1.000000 -1.000000 -2.000000 + 1.000000 -1.000000 -2.000000 + -1.000000 -1.000000 0.000000 + -1.000000 1.000000 0.000000 + -1.000000 1.000000 -2.000000 + -1.000000 -1.000000 -2.000000 + -1.000000 1.000000 0.000000 + 1.000000 1.000000 0.000000 + 1.000000 1.000000 -2.000000 + -1.000000 1.000000 -2.000000 + 1.000000 1.000000 0.000000 + 1.000000 -1.000000 0.000000 + 1.000000 -1.000000 -2.000000 + 1.000000 1.000000 -2.000000 + 1.000000 1.000000 -2.000000 + 1.000000 -1.000000 -2.000000 + -1.000000 -1.000000 -2.000000 + -1.000000 1.000000 -2.000000 + 2 3 4 1 + 6 7 8 5 + 10 11 12 9 + 14 15 16 13 + 18 19 20 17 diff --git a/tests/conftest.py b/tests/conftest.py index 1b5dd51..e26dd9b 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -566,3 +566,8 @@ def main_window(mocker, qtbot, tmpdir, install_lines): qtbot.addWidget(window) return window + + +@pytest.fixture +def test_data_folder(): + return os.path.join(this_dir, "..", "test_data") diff --git a/tests/test_wec_submodules_nemoh_run.py b/tests/test_wec_submodules_nemoh_run.py new file mode 100644 index 0000000..33517ea --- /dev/null +++ b/tests/test_wec_submodules_nemoh_run.py @@ -0,0 +1,56 @@ + +import numpy as np +import pytest +import logging + +from dtocean_wec.submodule.nemoh_run import _get_cylinder_radius +from dtocean_wec.submodule.utils.mesh import MeshBem + + +def test_get_cylinder_radius_no_meshs(): + + with pytest.raises(RuntimeError) as excinfo: + _get_cylinder_radius([]) + + assert "Can not determine mesh extents" in str(excinfo) + + +def test_get_cylinder_radius_cube(test_data_folder): + + mesh = MeshBem("cube.GDF", test_data_folder) + test = _get_cylinder_radius([mesh]) + expected = 2 + np.sqrt(2) + + assert np.isclose(test, expected) + + +def test_get_cylinder_radius_cube_offset(caplog, test_data_folder): + + mesh = MeshBem("cube.GDF", test_data_folder) + mesh.translate(1, 0, 0) + + with caplog.at_level(logging.WARNING): + test = _get_cylinder_radius([mesh]) + + expected_d = 2 + np.sqrt(5) + expected_log = "centroid of the mesh file is not centered" + + assert np.isclose(test, expected_d) + assert expected_log in caplog.record_tuples[0][2] + + +def test_get_cylinder_radius_two_cubes(caplog, test_data_folder): + + mesh1 = MeshBem("cube.GDF", test_data_folder) + mesh1.translate(1, 0, 0) + + mesh2 = MeshBem("cube.GDF", test_data_folder) + mesh2.translate(-1, 0, 0) + + with caplog.at_level(logging.WARNING): + test = _get_cylinder_radius([mesh1, mesh2]) + + expected_d = 2 + np.sqrt(5) + + assert np.isclose(test, expected_d) + assert not caplog.record_tuples diff --git a/tests/test_wec_submodules_utils_mesh.py b/tests/test_wec_submodules_utils_mesh.py new file mode 100644 index 0000000..3dc372d --- /dev/null +++ b/tests/test_wec_submodules_utils_mesh.py @@ -0,0 +1,392 @@ + +import os + +import numpy as np +import pytest +import matplotlib.pyplot as plt + +from dtocean_wec.submodule.utils.mesh import (strip_comments, + read_NEMOH, + read_WAMIT, + MeshBem, + _get_panel_norm, + Panel) + + +@pytest.fixture +def cube_xsim(): + return 0 + + +@pytest.fixture +def cube_vertices(): + return np.array([[ 1., -1., 0.], + [-1., -1., 0.], + [-1., -1., -2.], + [ 1., -1., -2.], + [-1., -1., 0.], + [-1., 1., 0.], + [-1., 1., -2.], + [-1., -1., -2.], + [-1., 1., 0.], + [ 1., 1., 0.], + [ 1., 1., -2.], + [-1., 1., -2.], + [ 1., 1., 0.], + [ 1., -1., 0.], + [ 1., -1., -2.], + [ 1., 1., -2.], + [ 1., 1., -2.], + [ 1., -1., -2.], + [-1., -1., -2.], + [-1., 1., -2.]]) + + +def test_strip_comments(): + + test = ("This is a multi-line string\n" + "! With some whole line comments\n" + " ! And some blank space before a whole line comment\n" + "And some comments that ! come after the token\n" + "But there are not comments here.") + expected = ("This is a multi-line string\n" + "And some comments that\n" + "But there are not comments here.") + + assert strip_comments(test) == expected + + +@pytest.mark.parametrize("file_name", ["cube.dat", "cube2.dat"]) +def test_read_nemoh(test_data_folder, + file_name, + cube_xsim, + cube_vertices): + + cube_connectivity = np.array([[ 2, 3, 4, 1], + [ 6, 7, 8, 5], + [10, 11, 12, 9], + [14, 15, 16, 13], + [18, 19, 20, 17]], dtype=int) - 1 + + file_path = os.path.join(test_data_folder, file_name) + xsim, vertices, connectivity = read_NEMOH(file_path) + + assert xsim == cube_xsim + assert (vertices == cube_vertices).all() + assert (connectivity == cube_connectivity).all() + + +@pytest.mark.parametrize("file_name", ["cube.GDF", "cube2.GDF"]) +def test_read_wamit(test_data_folder, + file_name, + cube_xsim, + cube_vertices): + + cube_connectivity = np.array([[ 0, 1, 2, 3], + [ 4, 5, 6, 7], + [ 8, 9, 10, 11], + [12, 13, 14, 15], + [16, 17, 18, 19]], dtype=int) + + file_path = os.path.join(test_data_folder, file_name) + xsim, vertices, connectivity = read_WAMIT(file_path) + + assert xsim == cube_xsim + assert (vertices == cube_vertices).all() + assert (connectivity == cube_connectivity).all() + + +@pytest.mark.parametrize("file_name", ["cube.GDF", "cube2.GDF"]) +def test_mesh_bem_init_gdf(file_name, + test_data_folder, + cube_xsim, + cube_vertices): + + cube_connectivity = np.array([[ 0, 1, 2, 3], + [ 4, 5, 6, 7], + [ 8, 9, 10, 11], + [12, 13, 14, 15], + [16, 17, 18, 19]], dtype=int) + mesh_bem = MeshBem(file_name, test_data_folder) + + assert mesh_bem.file_name == file_name + assert mesh_bem.path == test_data_folder + assert mesh_bem.xsim == cube_xsim + assert (mesh_bem.vertices == cube_vertices).all() + assert (mesh_bem.connectivity == cube_connectivity).all() + assert mesh_bem.nP == len(cube_connectivity) + assert mesh_bem.nV == len(cube_vertices) + + +@pytest.mark.parametrize("file_name", ["cube.dat", "cube2.dat"]) +def test_mesh_bem_init_dat(file_name, + test_data_folder, + cube_xsim, + cube_vertices): + + cube_connectivity = np.array([[ 2, 3, 4, 1], + [ 6, 7, 8, 5], + [10, 11, 12, 9], + [14, 15, 16, 13], + [18, 19, 20, 17]], dtype=int) - 1 + mesh_bem = MeshBem(file_name, test_data_folder) + + assert mesh_bem.file_name == file_name + assert mesh_bem.path == test_data_folder + assert mesh_bem.xsim == cube_xsim + assert (mesh_bem.vertices == cube_vertices).all() + assert (mesh_bem.connectivity == cube_connectivity).all() + assert mesh_bem.nP == len(cube_connectivity) + assert mesh_bem.nV == len(cube_vertices) + + +def test_mesh_bem_bad_ext(): + + with pytest.raises(IOError) as excinfo: + MeshBem("bad.bad") + + assert "file type not supported" in str(excinfo.value) + + +@pytest.fixture +def mesh_bem(test_data_folder): + return MeshBem("cube.GDF", test_data_folder) + + +@pytest.fixture +def mesh_bem_trans(mesh_bem): + mesh_bem.translate(1, 1, -1) + return mesh_bem + + +def test_mesh_bem_translate(mesh_bem_trans): + assert (mesh_bem_trans.vertices == np.array([[ 2., 0., -1.], + [ 0., 0., -1.], + [ 0., 0., -3.], + [ 2., 0., -3.], + [ 0., 0., -1.], + [ 0., 2., -1.], + [ 0., 2., -3.], + [ 0., 0., -3.], + [ 0., 2., -1.], + [ 2., 2., -1.], + [ 2., 2., -3.], + [ 0., 2., -3.], + [ 2., 2., -1.], + [ 2., 0., -1.], + [ 2., 0., -3.], + [ 2., 2., -3.], + [ 2., 2., -3.], + [ 2., 0., -3.], + [ 0., 0., -3.], + [ 0., 2., -3.]])).all() + + +def test_mesh_bem_rotate(mesh_bem_trans): + + mesh_bem_trans.rotate(-np.pi / 2, (0, 0, 0)) + expected = np.array([[ 0., -2., -1.], + [ 0., 0., -1.], + [ 0., 0., -3.], + [ 0., -2., -3.], + [ 0., 0., -1.], + [ 2., 0., -1.], + [ 2., 0., -3.], + [ 0., 0., -3.], + [ 2., 0., -1.], + [ 2., -2., -1.], + [ 2., -2., -3.], + [ 2., 0., -3.], + [ 2., -2., -1.], + [ 0., -2., -1.], + [ 0., -2., -3.], + [ 2., -2., -3.], + [ 2., -2., -3.], + [ 0., -2., -3.], + [ 0., 0., -3.], + [ 2., 0., -3.]]) + + assert np.isclose(mesh_bem_trans.vertices, expected).all() + + +def test_mesh_bem_invert_norm_all(mesh_bem): + + expected = np.array([[ 3, 2, 1, 0], + [ 7, 6, 5, 4], + [11, 10, 9, 8], + [15, 14, 13, 12], + [19, 18, 17, 16]]) + + mesh_bem.invert_norm() + assert (mesh_bem.connectivity == expected).all() + + +def test_mesh_bem_invert_norm_single(mesh_bem): + + expected = np.array([[ 3, 2, 1, 0], + [ 4, 5, 6, 7], + [ 8, 9, 10, 11], + [12, 13, 14, 15], + [16, 17, 18, 19]], dtype=int) + + mesh_bem.invert_norm(0) + assert (mesh_bem.connectivity == expected).all() + + +def test_get_panel_norm_zpos(): + x = np.array([-1., 1., 1., -1.]) + y = np.array([-1., -1., 1., 1.]) + z = np.array([ 0., 0., 0., 0.]) + norm = _get_panel_norm(x, y, z) + assert np.isclose(norm, (0, 0, 1)).all() + + +def test_get_panel_norm_zneg(): + x = np.array([-1., -1., 1., -1.]) + y = np.array([-1., 1., 1., -1.]) + z = np.array([ 0., 0., 0., 0.]) + norm = _get_panel_norm(x, y, z) + assert np.isclose(norm, (0, 0, -1)).all() + + +def test_get_panel_norm_yneg(): + x = np.array([ 1., 1., -1., -1.]) + y = np.array([-1., -1., -1., -1.]) + z = np.array([-2., 0., 0., -2.]) + norm = _get_panel_norm(x, y, z) + assert np.isclose(norm, (0, -1, 0)).all() + + +def test_panel(): + vertices = np.array([[-1., -1., 0.], + [ 1., -1., 0.], + [ 1., 1., 0.], + [-1., 1., 0.]]) + panel = Panel(vertices) + assert np.isclose(panel.centroid, (0, 0, 0)).all() + assert np.isclose(panel.n, (0, 0, 1)).all() + + +def test_panel_reverse(): + vertices = np.array([[-1., -1., 0.], + [-1., 1., 0.], + [ 1., 1., 0.], + [ 1., -1., 0.]]) + panel = Panel(vertices) + assert np.isclose(panel.centroid, (0, 0, 0)).all() + assert np.isclose(panel.n, (0, 0, -1)).all() + + +def test_panel_again(): + vertices = np.array([[ 1., -1., -2.], + [ 1., -1., 0.], + [-1., -1., 0.], + [-1., -1., -2.]]) + panel = Panel(vertices) + assert np.isclose(panel.centroid, (0, -1, -1)).all() + assert np.isclose(panel.n, (0, -1, 0)).all() + + +def test_mesh_bem_visualise_mesh(mocker, mesh_bem): + mocker.patch("dtocean_wec.submodule.utils.mesh.plt.show") + mesh_bem.visualise_mesh() + plt.close("all") + + +def test_mesh_bem_visualise_norm(mocker, mesh_bem): + mocker.patch("dtocean_wec.submodule.utils.mesh.plt.show") + mesh_bem.visualise_norm() + plt.close("all") + + +def test_mesh_bem_mesh_generation_gdf(tmpdir, test_data_folder, mesh_bem): + + mesh_bem.mesh_generation("gdf", str(tmpdir)) + files = tmpdir.listdir() + + assert len(files) == 1 + assert files[0].ext == ".gdf" + + expected = os.path.join(test_data_folder, "cube.GDF") + + with open(str(files[0])) as f1, open(expected) as f2: + + for _ in xrange(3): + next(f1) + next(f2) + + assert int(next(f1)) == int(next(f2)) + + for line1, line2 in zip(f1, f2): + vals1 = [float(x) for x in line1.split()] + vals2 = [float(x) for x in line2.split()] + assert np.isclose(vals1, vals2).all() + + +def test_mesh_bem_mesh_generation_nemoh(tmpdir, test_data_folder, mesh_bem): + + mesh_bem.mesh_generation("nemoh", str(tmpdir)) + files = tmpdir.listdir() + + assert len(files) == 1 + assert files[0].ext == ".dat" + + expected = os.path.join(test_data_folder, "cube.dat") + + with open(str(files[0])) as f1, open(expected) as f2: + + vals1 = [int(x) for x in next(f1).split()] + vals2 = [int(x) for x in next(f2).split()] + + assert vals1 == vals2 + + for line1, line2 in zip(f1, f2): + + vals1 = [float(x) for x in line1.split()] + vals2 = [float(x) for x in line2.split()] + + if np.isclose(vals1, (0, 0, 0, 0)).all(): + break + + assert np.isclose(vals1, vals2).all() + + for line1, line2 in zip(f1, f2): + vals1 = [int(x) for x in line1.split()] + vals2 = sorted([int(x) for x in line2.split()]) + assert vals1 == vals2 + + +def test_mesh_bem_mesh_generation_mesh(tmpdir, test_data_folder, mesh_bem): + + mesh_bem.mesh_generation("mesh", str(tmpdir)) + files = tmpdir.listdir() + + assert len(files) == 1 + assert files[0].ext == ".dat" + + expected = os.path.join(test_data_folder, "cube2.dat") + + with open(str(files[0])) as f1, open(expected) as f2: + + n_vertices = int(next(f1)) + + assert n_vertices == int(next(f2)) + assert int(next(f1)) == int(next(f2)) + + for _ in range(n_vertices): + vals1 = [float(x) for x in next(f1).split()] + vals2 = [float(x) for x in next(f2).split()] + assert np.isclose(vals1, vals2).all() + + for line1, line2 in zip(f1, f2): + vals1 = [int(x) for x in line1.split()] + vals2 = sorted([int(x) for x in line2.split()]) + assert vals1 == vals2 + + +def test_mesh_bem_mesh_generation_not_supported(mesh_bem): + + with pytest.raises(ValueError) as excinfo: + mesh_bem.mesh_generation("bad") + + assert "Unsupported output format" in str(excinfo.value)