Put everything in its right place

FEMOL/FEM.py

@@ -192,17 +192,20 @@ def solve(self, verbose=True, filtre=0, solve_from_zc=False, h_min=0, modal_solv
         :return: FEM Result class associated to the Problem kind
         """
         if self.physics == 'displacement':
-            return self._displacement_solve(verbose=verbose)
+            return self._displacement_solve(verbose=verbose, solve_from_zc=solve_from_zc, h_min=h_min)
 
         elif self.physics == 'modal':
             return self._modal_solve(verbose=verbose, solve_from_zc=solve_from_zc, filtre=filtre, h_min=h_min,
                                      modal_solver=modal_solver, sigma=sigma)
 
-    def _displacement_solve(self, verbose):
+    def _displacement_solve(self, verbose, solve_from_zc=False, h_min=0):
         """ General displacement solving function to allow multiple handling linear solvers"""
         # Assemble if not assembled
         if not hasattr(self, 'K'):
-            self.assemble('K')
+            if not solve_from_zc:
+                self.assemble('K')
+            else:
+                self._assemble_K(user_data=self._K_variable_core_laminate_data(h_min=h_min))
         # Solve with scipy
         U = self._scipy_displacement_solve(verbose=verbose)
         # add the solution to the mesh

FEMOL/TOM.py

@@ -152,16 +152,6 @@ def save_TOM_iteration(self, msh_file, eig_file, vec_file):
         else:
             np.save(eig_file, np.array(self.lmbds))
 
-    def save_guitar_modes(self, vecfile, eigfile):
-        """
-        Function saving the guitar modes from the modal solve result
-        :param eigfile: file base name for eigenfrequencies
-        :param vecfile: file base name for eigenvectors
-        :return: None
-        """
-        np.save(self.save_root + vecfile, np.array(self.guit_vecs))
-        np.save(self.save_root + eigfile, np.array(self.guit_freqs))
-
     def check_convergence(self, converge, min_iter, max_iter, convergence_criteria):
         """
         Method to check the convergence for the current iteration

FEMOL/domains.py

@@ -131,9 +131,7 @@ def ellipse(x):
 
     # Ellipse domain function
     def domain(x, y):
-
         Ri = (x - h) ** 2 / (a - eps) ** 2 + (y - k) ** 2 / (b - eps) ** 2
-
         if Ri < 1:
             return False
         else:
@@ -153,8 +151,8 @@ def outside_guitar(L):
 
     # Right ellipse
     elc2 = (0.8175 * L, 0.38 * L)
-    elsa2 = (0.8175 * L, 0.09)
-    elso2 = (1, 0.38 * L)
+    elsa2 = (0.8175 * L, 0.09*L)
+    elso2 = (1*L, 0.38 * L)
     ellipse2 = outside_ellipse(elc2, elsa2, elso2)
 
     # Top and bottom curves
@@ -249,3 +247,8 @@ def fixed_guitar(x, y):
             return False
 
     return fixed_guitar
+
+
+def top_brace(L=1):
+    domain = inside_box([[0.8175 * L - (L/30), 0.8175 * L + (L/20)]], [[0, L]])
+    return domain

FEMOL/materials.py

@@ -134,9 +134,6 @@ def G_mat(self, angle):
                          [0, Gyz]])
 
 
-
-
-
 class IsotropicMaterial(object):
 
     def __init__(self, E, mu, rho):

FEMOL/mesh.py

@@ -575,7 +575,7 @@ def point_data_to_STL(self, filename, which, hb=1, hc=0, symmetric=False, scale=
             else:
                 # if there is only a base
                 z_top = hb * self.point_data[which]
-                z_bot = np.ones((self.points.shape[0])) * self.point_data[which].min()
+                z_bot = np.zeros(self.points.shape[0])
 
             # points at the top
             points1[:, 2] = z_top
@@ -616,7 +616,6 @@ def point_data_to_STL(self, filename, which, hb=1, hc=0, symmetric=False, scale=
 
         tris = np.vstack([tris1, tris2, tris3])
         cells = [('triangle', tris)]
-        points *= scale
         meshio_mesh = meshio.Mesh(points, cells)
         meshio_mesh.write(filename + '.stl')
 
@@ -1054,7 +1053,6 @@ def guitar_sym(L=1, lcar=0.05, option='quad', algorithm=1):
         # Ellipse 2
         elc2 = (0.8175 * L, 0.38 * L)
         elh2 = 0.58 * L
-
         # Left top ellipse
         elsa1 = geom.add_point(p0, lcar)
         elce1 = geom.add_point(elc1, lcar)

FEMOL/problems.py

@@ -2,6 +2,59 @@
 import numpy as np
 
 
+class GuitarDeflexion(object):
+
+    # materials
+    n_plies_carbon = 2
+    n_plies_flax = 9
+    plies_flax = [0, 0, 0, 0, 0, 0, 0, 0, 0]
+    flax = FEMOL.materials.general_flax()  # material from library
+    carbon = FEMOL.materials.general_carbon()
+
+    def __init__(self, lcar=0.03,
+                 h_flax=0.0025, h_carbon=0.0002, hc_opt=0.010,
+                 plies_carbon=None, factor=0.480, core_data_file=None):
+
+        # mesh
+        self.mesh = FEMOL.mesh.guitar_sym(L=1, lcar=lcar)
+        self.factor = factor
+        self.core_file = core_data_file
+
+        # Ply thickness
+        self.flax.hi = h_flax / self.n_plies_flax
+        self.carbon.hi = h_carbon / self.n_plies_carbon
+        # layups
+        h = hc_opt + h_flax + h_carbon
+        if plies_carbon is None:
+            plies_carbon = [0, 90]
+        z_flax = -h / 2 + (self.n_plies_flax / 2) * self.flax.hi
+        self.flax_layup = FEMOL.laminate.Layup(material=self.flax, plies=self.plies_flax, symetric=False, z_core=z_flax)
+        z_carbon = h / 2 - (self.n_plies_carbon / 2) * self.carbon.hi
+        self.carbon_layup = FEMOL.laminate.Layup(material=self.carbon, plies=plies_carbon, symetric=False, z_core=z_carbon)
+
+    def construct_FEM_problem(self):
+
+        # Create a FEM Problem from the mesh (compute displacement with a plate bending model)
+        problem = FEMOL.FEM_Problem(mesh=self.mesh, physics='displacement', model='plate')
+
+        problem.add_fixed_domain(FEMOL.domains.outside_guitar(L=1))
+        problem.define_materials(self.flax, self.carbon)
+        problem.define_tensors(self.flax_layup, self.carbon_layup)
+        x_values = [0.135 / self.factor, 0.165 / self.factor]
+        y_values = [0.095 / self.factor, 0.265 / self.factor]
+        bridge = FEMOL.domains.inside_box([x_values], [y_values])
+        # string tension
+        problem.add_forces([0, 0, 0, -21, 0, 0], bridge)
+        zc_mesh = FEMOL.mesh.load_vtk(self.core_file)
+        zc_mesh.cell_data['zc'] = FEMOL.utils.add_top_brace_to_zc(zc_mesh, zc_mesh.cell_data['zc'], 0.010)
+        problem.mesh.cell_data['zc'] = zc_mesh.cell_data['zc']
+        return problem
+
+    def solve(self, h_min=0):
+        problem = self.construct_FEM_problem()
+        return problem.solve(solve_from_zc=True, h_min=h_min)
+
+
 class GuitarSimpVibe(object):
     """ Class containing a guitar Topology optimization problem"""
     # problem data
@@ -67,42 +120,39 @@ def solve(self, save=True, plot=False, verbose=True, **kwargs):
 
 class GuitarModal(object):
     """ Class representing a guitar modal analysis problem"""
-    hc_opt = 0.010  # optimized core thickness
-    h_flax = 0.003
-    h_carbon = 0.000250
     n_plies_carbon = 2
     n_plies_flax = 9
-
-    # flax material definition
-    flax = FEMOL.materials.general_flax()  # material from library
-    flax.hi = h_flax / n_plies_flax
     # carbon material definition
-    carbon = FEMOL.materials.general_carbon()
-    carbon.hi = h_carbon/n_plies_carbon
-
-    # Laminates definitions
-    h = hc_opt + h_flax + h_carbon
-    # Reference layups
-    plies_flax = [0, 0, 0, 0, 0, 0, 0, 0, 0]
-    flax_layup = FEMOL.laminate.Layup(material=flax, plies=plies_flax, symetric=False,
-                                      z_core=-h / 2 + (n_plies_flax / 2) * flax.hi)
-    plies_carbon = [90, 90]
-    carbon_layup = FEMOL.laminate.Layup(material=carbon, plies=plies_carbon, symetric=False,
-                                        z_core=h / 2 - (n_plies_carbon / 2) * carbon.hi)
 
-    def __init__(self, mesh_lcar=0.04):
+    def __init__(self, mesh_lcar=0.03,
+                 hc_opt=0.010, h_flax=0.003, h_carbon=0.000250,
+                 plies_flax=None, plies_carbon=None):
+        # Mesh definition
         self.mesh = FEMOL.mesh.guitar_sym(lcar=mesh_lcar)
         self.lcar = mesh_lcar
-
+        # Materials and layups
+        flax = FEMOL.materials.general_flax()  # material from library
+        flax.hi = h_flax / self.n_plies_flax
+        carbon = FEMOL.materials.general_carbon()
+        carbon.hi = h_carbon / self.n_plies_carbon
+        h = hc_opt + h_flax + h_carbon
+        if plies_flax is None:
+            plies_flax = [0, 0, 0, 0, 0, 0, 0, 0, 0]
+        if plies_carbon is None:
+            plies_carbon = [90, 90]
+        flax_layup = FEMOL.laminate.Layup(material=flax, plies=plies_flax, symetric=False,
+                                          z_core=-h / 2 + (self.n_plies_flax / 2) * flax.hi)
+        carbon_layup = FEMOL.laminate.Layup(material=carbon, plies=plies_carbon, symetric=False,
+                                            z_core=h / 2 - (self.n_plies_carbon / 2) * carbon.hi)
         # FEM problems definition
         self.problem = FEMOL.FEM_Problem(mesh=self.mesh, physics='modal', model='plate')
-        self.problem.define_materials(self.flax, self.carbon)
-        self.problem.define_tensors(self.flax_layup, self.carbon_layup)
+        self.problem.define_materials(flax, carbon)
+        self.problem.define_tensors(flax_layup, carbon_layup)
         self.problem.add_fixed_domain(FEMOL.domains.outside_guitar(L=1))
 
-    def solve(self, mac_find=True):
+    def solve(self, mac_find=True, **kwargs):
         # Solve the reference problem
-        w_opt, v_opt = self.problem.solve()
+        w_opt, v_opt = self.problem.solve(**kwargs)
         if mac_find:
             # Find guitar reference eigenfrequencies and vectors
             modes = ['T11', 'T21', 'T12', 'T31']

FEMOL/utils.py

@@ -4,7 +4,7 @@
 import matplotlib.pyplot as plt
 from matplotlib import animation
 import datetime
-from scipy.interpolate import griddata
+from scipy.interpolate import griddata, interp1d
 
 """
 Plot/Validation functions
@@ -196,13 +196,13 @@ def guitar_outline2(L):
     y = p[0] * x ** 3 + p[1] * x ** 2 + p[2] * x + p[3]
     ax.plot(x, y, color='k')
     # Top 2
-    p = FEMOL.domains.create_polynomial(0.625 * L, 0.71225 - 0.1645 / 2, 0.8175 * L, 0.67, 0)
-    x = np.linspace(0.625, 0.8175, 100)
+    p = FEMOL.domains.create_polynomial(0.625 * L, 0.71225*L - 0.1645*L / 2, 0.8175 * L, 0.67*L, 0)
+    x = np.linspace(0.625*L, 0.8175*L, 100)
     y = p[0] * x ** 3 + p[1] * x ** 2 + p[2] * x + p[3]
     ax.plot(x, y, color='k')
     # Bot 1
-    p = FEMOL.domains.create_polynomial(0.625 * L, 0.04775 + 0.1645 / 2, 0.8175 * L, 0.09, 0)
-    x = np.linspace(0.625, 0.8175, 100)
+    p = FEMOL.domains.create_polynomial(0.625 * L, 0.04775*L + 0.1645*L / 2, 0.8175 * L, 0.09*L, 0)
+    x = np.linspace(0.625*L, 0.8175*L, 100)
     y = p[0] * x ** 3 + p[1] * x ** 2 + p[2] * x + p[3]
     ax.plot(x, y, color='k')
     # Bot 2
@@ -497,6 +497,17 @@ def interpolate_vector(old_vector, old_mesh, new_mesh, N_dof=6):
     return new_vector
 
 
+def interpolate_point_data(old_data, old_mesh, new_mesh):
+    """ interpolate point data onto a new mesh"""
+    # Interpolate at the new mesh points
+    data_linear = griddata(old_mesh.points[:, :2], old_data, new_mesh.points[:, :2], method='linear')
+    data_near = griddata(old_mesh.points[:, :2], old_data, new_mesh.points[:, :2], method='nearest')
+    new_data = data_linear
+    # Use the nearest value where the linear value is NaN (boundary)
+    new_data[np.isnan(data_linear)] = data_near[np.isnan(data_linear)]
+    return new_data
+
+
 def analyse_mesh(meshfile, eigvalfile=None, mode=None):
     # Load the data
     mesh = FEMOL.mesh.load_vtk(meshfile)
@@ -521,3 +532,105 @@ def analyse_mesh(meshfile, eigvalfile=None, mode=None):
         title = f'TOM results {key} eigfreq {int(np.round(eig))} Hz'
         ax.set_title(title)
 
+
+def plot_soundboard_deflexion(mesh):
+    """ Function to plot the soundboard deflexion analysis"""
+
+    def align_yaxis_np(ax1, ax2):
+        """Align zeros of the two axes, zooming them out by same ratio"""
+        axes = np.array([ax1, ax2])
+        extrema = np.array([ax.get_ylim() for ax in axes])
+        tops = extrema[:, 1] / (extrema[:, 1] - extrema[:, 0])
+        # Ensure that plots (intervals) are ordered bottom to top:
+        if tops[0] > tops[1]:
+            axes, extrema, tops = [a[::-1] for a in (axes, extrema, tops)]
+
+        # How much would the plot overflow if we kept current zoom levels?
+        tot_span = tops[1] + 1 - tops[0]
+
+        extrema[0, 1] = extrema[0, 0] + tot_span * (extrema[0, 1] - extrema[0, 0])
+        extrema[1, 0] = extrema[1, 1] + tot_span * (extrema[1, 0] - extrema[1, 1])
+        [axes[i].set_ylim(*extrema[i]) for i in range(2)]
+
+    x_points = mesh.points[np.isclose(mesh.points[:, 1], (0.76 / 2)), 0]
+    idxs = np.argsort(x_points)
+    displacement = mesh.point_data['Uz'][np.isclose(mesh.points[:, 1], (0.76 / 2))][idxs]
+    x_points = x_points[idxs]
+
+    f = interp1d(x_points, displacement, kind='quadratic')
+    x1 = np.linspace(0, 0.6)
+    x2 = np.linspace(0.78, 1)
+    T1 = -np.degrees(np.arctan2(f(x1)[1:] - f(x1)[:-1], x1[1:] - x1[:-1]))
+    T2 = -np.degrees(np.arctan2(f(x2)[1:] - f(x2)[:-1], x2[1:] - x2[:-1]))
+    T_max = np.max(np.abs(np.hstack([T1, T2])))
+
+    fig, ax1 = plt.subplots()
+    ax1.plot(x1, f(x1) * 1000, color='#743F0B', label='deflexion')
+    ax1.plot(x2, f(x2) * 1000, color='#743F0B')
+    ax1.plot([0, 1], [0, 0], '--', color='0.6')
+    ax1.set_ylim(-2, 5)
+    ax1.plot([], [], linestyle='--', color='k', label='angle')
+    ax1.set_xlabel('normalized distance across soundboard')
+
+    ax2 = ax1.twinx()
+    ax2.plot(x1[1:], T1, linestyle='--', color='k', label='angle')
+    ax2.plot(x2[1:], T2, linestyle='--', color='k')
+
+    ax1.legend()
+    ax1.set_ylabel('Deflexion (mm)')
+    ax2.set_ylabel('Normal angle (degrees)')
+    ax1.grid('on')
+    ax1.set_xlim(0, 1)
+    align_yaxis_np(ax1, ax2)
+
+    return T_max
+
+
+def topology_info(mesh, h_min=0, scale=0.480):
+    """ Function printing the topology info from a mesh"""
+    Mf = flax_mass(mesh)
+    Vcore = core_volume(mesh, h_min=h_min)
+    Mpla = Vcore * 1.25 * 0.2
+    print(f'pla mass: {np.around(Mpla, 1)} g')
+    Mcarbon = carbon_mass(mesh, h_min=h_min)
+    print(f'ratio brace/board {np.around((Mcarbon + Mpla) / Mf, 2)}')
+
+
+def core_volume(mesh, h_min=0.002, scale=0.480):
+    """ Function computing the volume of a core of a topology result"""
+    h = mesh.cell_data['zc']['quad']
+    h[h < h_min] = 0
+    A = mesh.element_areas()['quad'] * (scale ** 2)
+    print(f'core volume: {np.around(np.sum(h * A) * (100 ** 3), 2)} cm^3')
+    return np.sum(h * A) * (100 ** 3)
+
+
+def flax_mass(mesh, t=0.0025, scale=0.480):
+    """ Function computing the flax mass of a board from a mesh"""
+    A = mesh.element_areas()['quad'] * (scale ** 2)
+    flax = FEMOL.materials.general_flax()
+    V = np.sum(A * t)  # m3
+    M = V * flax.rho
+    return 1000 * M
+
+
+def carbon_mass(mesh, t=0.00025, h_min=0, scale=0.480):
+    """ Fonction computing the carbon mass from a topology result"""
+    A = mesh.element_areas()['quad'] * (scale ** 2)
+    h = mesh.cell_data['zc']['quad']
+    flax = FEMOL.materials.general_carbon()
+    V = np.sum(A[h > h_min] * t)  # m3
+    M = V * flax.rho
+    M *= 1000
+    print(f'carbon mass: {np.around(M, 1)} g')
+    return M
+
+
+def add_top_brace_to_zc(mesh, zc, hc):
+    """ Add a top brace to the core height cell data of a mesh"""
+    domain = FEMOL.domains.top_brace(L=1)
+    for i, element in enumerate(mesh.cells[mesh.contains[0]]):
+        if np.array([domain(*coord[:2]) for coord in mesh.points[element]]).all():
+            zc[mesh.contains[0]][i] = hc
+    return zc
+