Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Hydrostatics for fully submerged body. #19

Open
akeeste opened this issue Nov 23, 2020 · 4 comments
Open

Hydrostatics for fully submerged body. #19

akeeste opened this issue Nov 23, 2020 · 4 comments

Comments

@akeeste
Copy link

akeeste commented Nov 23, 2020

This issue relates exactly to #7 . I am having the same issue and would like to work around it. I am also using fully submerged body and have dug into the error a bit. When the hydrostatics class is initialized, it splits the full mesh into the lower portion and the 'clipped mesh' (crown mesh) that represents the portion of the body above the water. The crown mesh is initialized as empty with no vertices / faces, and they are added in mesh_clipper._clip_crown_by_plane(). If the mesh is fully submerged, no vertices/faces end up being added and meshmagick attempts to create a meshmagick.Mesh() with no vertices/faces. This ultimately fails with an assert error. Can meshmagick catch this error and avoid creating a clipped mesh altogether? then it can return 0 stiffness and the rest of the desired properties from the Hydrostatics module.

I am currently interested in the displaced volume, center of buoyancy, and hydrostatic stiffness. When calling meshmagick.Hydrostatics, I can catch this error and assume that the stiffness is zero for a submerged body. The displaced volume can be calculated with the meshmagick.mesh.Mesh.volume module. However there seems to be no way to get the center of buoyancy outside of the Hydrostatics module. Can I calculate this with another meshmagick function?

File "C:\Users\akeeste\Anaconda3\envs\py37\lib\site-packages\meshmagick\hydrostatics.py", line 165, in init
self._update_hydrostatic_properties()
File "C:\Users\akeeste\Anaconda3\envs\py37\lib\site-packages\meshmagick\hydrostatics.py", line 564, in _update_hydrostatic_properties
clipper = MeshClipper(self.mesh, assert_closed_boundaries=True, verbose=False)
File "C:\Users\akeeste\Anaconda3\envs\py37\lib\site-packages\meshmagick\mesh_clipper.py", line 37, in init
self._update()
File "C:\Users\akeeste\Anaconda3\envs\py37\lib\site-packages\meshmagick\mesh_clipper.py", line 114, in _update
self._clip()
File "C:\Users\akeeste\Anaconda3\envs\py37\lib\site-packages\meshmagick\mesh_clipper.py", line 751, in _clip
self._clip_crown_by_plane()
File "C:\Users\akeeste\Anaconda3\envs\py37\lib\site-packages\meshmagick\mesh_clipper.py", line 665, in _clip_crown_by_plane
clipped_crown_mesh = Mesh(vertices, crown_faces)
File "C:\Users\akeeste\Anaconda3\envs\py37\lib\site-packages\meshmagick\mesh.py", line 409, in init
assert np.array(faces).shape[1] == 4
IndexError: tuple index out of range

@ryancoe
Copy link
Contributor

ryancoe commented Mar 6, 2021

@akeeste - Not sure if you already solved this, but I recently stumbled across a workable solution using the methods in meshmagick that generate RigidBodyInertia objects. For example, you can use the eval_plain_mesh_inertias method as follows.

def eval_plain_mesh_inertias(self, rho_medium=1023.):

from meshmagick.mmio import load_mesh
from meshmagick.mesh import Mesh
import pygmsh
import pytest

r = 1
geom = pygmsh.opencascade.Geometry()
sphere = geom.add_ball([0,0,-1], r)
mshArgs = ['-clscale', str(0.25),
           '-clcurv', str(360/50)]

fbname = 'mmtest'
_ = pygmsh.generate_mesh(geom,
                            dim=2,
                            extra_gmsh_arguments=mshArgs,
                            remove_lower_dim_cells=True,
                            geo_filename=fbname + '.geo',
                            msh_filename=fbname + '.stl',
                            mesh_file_type="stl")

mesh = Mesh(*load_mesh(fbname + '.stl', 'stl'))


inertia = mesh.eval_plain_mesh_inertias()
cob = inertia.gravity_center

assert pytest.approx(-1, rel=1e-4) == cob[2]

image

@akeeste
Copy link
Author

akeeste commented Mar 8, 2021

@ryancoe Thanks for pointing this out! I hadn't found a solution for this yet, this is very helpful. Now in the instance where the body is fully submerged, I can catch the hydrostatics error and just use the mesh to get the cb.

@akeeste akeeste closed this as completed Mar 8, 2021
@ryancoe
Copy link
Contributor

ryancoe commented Mar 8, 2021

@akeeste - Happy to help. Would you mind posting an example of your workflow here? Or (even better) submit a PR that would make this seamless?

@akeeste
Copy link
Author

akeeste commented Mar 8, 2021

@ryancoe Good idea, the python workflow I use to get around this issue is below. The 'body' object is a Capytaine Floating Body, though of course these parameters could be input however required. For reference, my entire Capytaine workflow can be found in the WEC-Sim repo here.

I'll work on a PR to implement this in meshmagick and resolve this and issue #7 simultaneously.

import meshmagick.mesh as mmm
import meshmagick.hydrostatics as mmhs

# use meshmagick to compute hydrostatic stiffness matrix
       # NOTE: meshmagick currently has issue if a body is completely submerged (e.g. OSWEC base)
       # See issues #7, #19 on meshmagick GitHub page.
       # use try-except statement to catch this error and use alternate function for cb/vo
body_mesh = mmm.Mesh(body.mesh.vertices, body.mesh.faces, name=body.mesh.name)
try:
    body_hs = mmhs.Hydrostatics(working_mesh=body_mesh,
                                cog=body.center_of_mass,
                                rho_water=1023.0,
                                grav=9.81)
    vo = body_hs.displacement_volume
    cb = body_hs.buoyancy_center
    khs = body_hs.hydrostatic_stiffness_matrix
except:
    # Exception if body is fully submerged.
    # if fully submerged:
        # stiffness is 0 as small displacements do not change the displaced volume or hydrostatic force
        # displaced volume is the total mesh volume
        # center of buoyancy is equal to the center of gravity of the mesh with a constant density 
    vo = body_mesh.volume
    khs = np.zeros((3,3))
    inertia = body_mesh.eval_plain_mesh_inertias()
    cb = inertia.gravity_center

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants