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)