Arbitrary geometry and mesh

[ryancoe]

We discussed this a few weeks ago... currently we can generate axisymmetric meshes, but we will need to find a means of generating arbitrary geometries and meshes. Some ideas:

• capytaine
• Gmsh
• meshio
• meshmagick
• dual-mesh

[ryancoe]

I was tinkering with Gmsh this morning for another project and I think this may be a very good option. I was working on a TLP and was able to create a parametric geometry fairly easily. I used the GUI to construct the basic form, and then parameterized this using the automatically generated ASCII input file (tlpBase.geo):

//+
SetFactory("OpenCASCADE");
lambda = 10;
r_spar = 0.125 * lambda;
w_pontoon = 0.1 * lambda;
l_pontoon = 0.75 * lambda;
T = 1 * lambda;
ofst = 0.1 * lambda;
K = ofst - T;
Cylinder(1) = {0, 0, ofst, 0, 0, K, r_spar, 2*Pi};
//+
Box(2) = {-w_pontoon/2, 0, K+w_pontoon, w_pontoon, l_pontoon, w_pontoon};
//+
Rotate {{0, 0, 1}, {0, 0, 0}, -2*Pi/3} {
  Duplicata { Volume{2}; }
}
//+
Rotate {{0, 0, 1}, {0, 0, 0}, 2*Pi/3} {
  Duplicata { Volume{2}; }
}
//+
BooleanUnion{ Volume{1}; Delete; }{ Volume{2}; Volume{3}; Volume{4}; Delete; }

To generate the mesh file, you can call: gmsh tlpBase.geo -clscale 0.25 -2 -o tlpBase.stl where clscale is a mesh refinement factor. This produces the mesh rendered below.

We could write a wrapper to generate the geo files. Alternatively, gmsh has a python interface that I we use. Another similar option is pygmsh, which I ran into problems with last time I looked at it (LHEEA/meshmagick#9).

Gmsh can output many mesh formats (STL, VTK, etc.). The native format (.msh) can be read by meshmagick (I'm having some trouble with this... but STL works ok) and then converted to NEMOH's or WAMIT's formats.

[ryancoe]

I've worked a bit more with pygmsh. Its a very nice tool and gives you somewhat cleaner programmatic access to the functionality of gmsh (programming the gmsh .geo files is pretty challenging).

However, one advantage of using the .geo files is that you can create a geometry via the gmsh GUI and then export a .geo file, which can latter be edited and/or parameterized. This is a potentially interesting approach for WecOptTool. The procedure would be

1. Create your geometry (probably via the gmsh GUI)
2. Export a .geo file (see example in my previous comment)
3. Name the parameters you want to use as design variables in a memorable way
4. Have a MATLAB function that goes in and edits the parameter values
5. Create the mesh via a system call (system('gmsh <geofile> -clscale <meshScaleFactor> -1 -2 -o <stlFileOutName>'))
6. Convert/check mesh (this may or may not be necessary, meshmagick is a great python tool to do this if needed)
7. Call Nemoh, etc.

[H0R5E]

Yeah, that's pretty much the procedure I had in mind. Mesh quality is definitely an issue. Is there some way we can check that the mesh is good enough?

[ryancoe]

We must ensure

1. Mesh normals all point outwards. I know meshmagick can check this. We can also doe this within MATLAB as well:

https://www.mathworks.com/help/matlab/ref/surfnorm.html
https://www.mathworks.com/matlabcentral/fileexchange/33226-compute-stl-vertex-normals

2. The mesh is open at the free surface -- don't know exactly how to do this, but certainly seems possible.

Mesh quality is another issue... in BEM, one easy thing I know we can check is the size of the cells relative to the wave lengths.

[H0R5E]

@ryancoe, this is going to sit outside of the architecture proposed in #126, which will take a mesh in some sort of general definition, so would you like to put something together to handle to procedure you outlined above? Maybe put it into a new package like +WecOptTool\+mesh?

I think the starting point is a function that takes a parameterised .geo file, the parameters to be modified, and some values for those parameters and spits out a mesh format (preferably a portable/generic one), which we can pass on to the BEM solver.

Would you mind putting that together?

[ryancoe]

@H0R5E - Yes, happy to take a cut at this. I will likely prototype it as simple functions and then look at your work and ask for help to make it OOP, etc.

[H0R5E]

Great. If you can use WaveBot as a case study then this could support #123, and we can pop this into the v0.1 milestone.

[ryancoe]

Messing more with Gmsh...

I can create the WaveBot in Gmsh by sweeping a series of curves about the vertical axis (analogous to what we're doing with MATLAB at the moment) using the following input (.geo file):

SetFactory("OpenCASCADE");

Mesh.Algorithm = 6; // 1: MeshAdapt, 2: Automatic, 5: Delaunay, 6: Frontal-Delaunay, 7: BAMG, 8: Frontal-Delaunay for Quads, 9: Packing of Parallelograms (default 6)
Mesh.CharacteristicLengthMin = 0.1;
Mesh.CharacteristicLengthMax = 0.1;

fb = 0;
r1 = 0.88;
Tc = 0.16;
r2 = 0.35;
Tk = 0.53;

Point(1) = {0,0,fb};
Point(2) = {r1,0,fb};
Point(3) = {r1,0,-Tc};
Point(4) = {r2,0,-Tk};
Point(5) = {0,0,-Tk};


Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,5};

Extrude {{0, 0, 1}, {0, 0, 0}, -2*Pi} {
  Curve{2}; Curve{3}; Curve{4}; 
}

Alternatively, you could build the geometry by combining a cylinder and a cone with something like:

SetFactory("OpenCASCADE");
vol0 = newv;
Cylinder(vol0) = {0, 0, 0, 0.0, 0.0, -Tc, r1, -6.283185307179586};
pts_vol0[] = PointsOf{Volume{vol0};};
Characteristic Length{pts_vol0[]} = 1;
vol1 = newv;
Cone(vol1) = {0, 0, -Tc, 0.0, 0.0, -Tk, r1, r2};
pts_vol1[] = PointsOf{Volume{vol1};};
Characteristic Length{pts_vol1[]} = 1;
bo1[] = BooleanUnion{ Volume{vol0}; Delete; } { Volume{vol1}; Delete;};

From there, the sky is the limit... for example you can use Bezier curves:

SetFactory("OpenCASCADE");

Mesh.Algorithm = 6; // 1: MeshAdapt, 2: Automatic, 5: Delaunay, 6: Frontal-Delaunay, 7: BAMG, 8: Frontal-Delaunay for Quads, 9: Packing of Parallelograms (default 6)
//Mesh.CharacteristicLengthMin = 0.05;
//Mesh.CharacteristicLengthMax = 0.1;
//Mesh.CharacteristicLengthFromCurvature = 360;
Mesh.CharacteristicLengthFactor = 0.2;


r1 = 0.88;
Tc = 0.16;
r2 = 0.35;
Tk = 0.53;


Point(1) = {r1, 0, 0};
//Point(2) = {r1, 0, -Tc/2};
Point(3) = {r1, 0, -Tc};
Point(4) = {0, 0, -Tc};
Point(5) = {r2, 0, -Tk*0.8};
Point(6) = {0, 0, -Tk};


Bezier(1) = {1, 3, 4, 5, 6};

Extrude {{0, 0, 1}, {0, 0, 0}, -2 * Pi} { // use negative to get normals right
  Curve{1}; 
}