Standard mesh for Concave Hexagons with two Mouths?

Question

I am planning a visualization of flows though Concave Bisymmetric hexagons with two mouths.

Example where the length of the side d1 equals the other length of the side d2:

which naming I discussed initially here about Irregular hexagons.

There is the standard Mesh tool where you can draw your own meshes but I would like have some standard library such that I can cooperate better with others with the flow simulations later. I did not find any Mesh library for Hexagons in MathCentral File Exchange here.

Is there any standard Mesh library for irregular hexagon shapes? I am open to any other languages too because I could read the code and convert it those standards to Matlab library.

Answer 1

I discussed this issue with my colleque. He motivates to create a triangle network but also points out general algorithms for nonconvex polygons.

To create a Mesh by own algorithm

1. divide the hexagon into two trapezoids by the dashed line as illustrated in the body of the question
2. divide trapezoids into triangles: assemmbly and construct trimesh structures by inittri
3. load vector f is built by bilin_assembly (because this worked before but not necessarily the optimal choice)

Division of hexagon with two mouths into two trapezoids

1. take two nodes which index difference is three and which distance is smallest among six nodes.
2. form two trapezoids such that one trapezoid has numbers 1-4 while the other 1,6,5,4.

The following syntaxes are based on my editions of codes in 2013 but mostly on descriptions of the book Numerical solutions of partial differential equations by the finite element method by Claes Johnson.

Syntax of assembly where considered now for simplicity a Poisson problem (must be later refined here)

% [S,f]=assemblyGlobal(loadfunction,mesh)
%
%   loadfunction = function on the right-hand side of the Poisson equation
%   mesh         = mesh structure on which the assembly is done
%
%   S            = stiffness matrix of the Poisson problem
%   f            = load vector corresponding to loadfunction

where we need bilinear assembly which syntax

% function B=bilin_assembly(bilin,mesh)
%
%  bilin = function handle to integrand of the bilinear form
%          given as bilin(u,v,ux,uy,vx,vy,x,y), where
%          u  - values of function u
%          v  - values of function v
%          ux - x-derivative of function u
%          uy - y-derivative of function u
%          vx - x-derivative of function v
%          vy - y-derivative of function v
%          x  - global x-coordinate
%          y  - global y-coordinate
%  mesh  = mesh structure on which the matrix will be assembled
%
%  B     = matrix related to bilinear form bilin

Syntax to init triangular mesh

% function mesh = inittri(p,t)
%
%   p       = nodes
%   t       = triangles
%     
%   mesh    = trimesh structure corresponding to (p,t)
%
% TRIMESH STRUCTURE :
% 
%   p       = nodes in a 2xN-matrix
%
%   t       = triangles in a 3xN- matrix (nodes of the triangles as indeces to the p- matrix).
% 
%   edges   = a matrix of all edges in the mesh. Each column is an edge :
%             [n1 ; n2] where n1 < n2;
%
%   t2e     = a matrix connecting  triangle's and edges's. 
%             Each column corresponds to a triangle and has
%             triangle's edges in the order n1->n2, n2->n3,
%             n1->n3. 
%
%   e2t     = inverse of t2e.

I did no post the complete source codes here because of copyright issues and they would make the answer rather long. However, all algorithms are based on the first source so can be created by any researcher here. I think this kind of FEM method is the only way to find the optimum mesh here.

Algorithms for general nonconvex polygons

There are also algorithms which can create meshes for general nonconvex polygons. It is possible that the mesh provided by Margus' answer is produced by such an algorithm.

Ansys

I was experiencing different products about visualizing holes, different geometries and different materials. Ansys is promising solution here.

Sources

1. Numerical solutions of partial differential equations by the finite element method by Claes Johnson.
2. Class notes in Finite Element Methods at my university, 2013-2014

Answer 2

If this is to stand-in as purely the presentation layer, then it really depends on the output of the simulation and the representation you select for it.

For example, if the simulation outputs one single variable (e.g. the volume of one of the compartments), then this could be assigned directly to an attribute of the visualisation such as the location of the "top" cap of a cylinder. A "complementary" cylinder could stand on top of the previous one with its bottom-cap coordinates assigned to (totalVolume-lowerCompartmentVolume).

In this case, the visualisation is just a dashboard and does not feed back into the simulation (e.g. object collisions or proximity in general are not taken into account by the simulation).

To generalise this, we are talking about a solution where a set of quantities is assigned to a set of visualisation attributes.

From this point of view, a set of complex objects could be created using VRML (or X3D), their attributes could be binded directly to simulation outputs and the rendering be triggered every n-th simulation timestep.

To create the "scene" or visualisation objects, you can use a software like blender which can export VRML scenes, or write the VRML manually (it is really a straightforward task).

In terms of infrastructure, MATLAB has a VRML toolkit and Python has a really wide selection of modules by which VRML can be handled (please see this and this link for example).

For a more concrete example:

Given some simulation output y and a template VRML file like:

#VRML V2.0 utf8
Transform {
  translation 0 0 0
  children [
    Shape { geometry Box {2,2,zSize} }
  ]
}

You could do something like:

data = (Read contents of VRML file as string data).
for n in [0..1000]:
    y = getSimulationOutput(aParameterVector)
    renderData = substitute(data, "zSize", y) #This function could be provided by a template module like jinja for example. 
    simulationFrame = renderVRML(renderData)
    saveImage(simulationFrame)

(Please note: More information about Jinja can be found here - comment inline link not rendering properly above.)

Going back into the file and binding different elements to different quantities (for example, changing the transform which can rotate, scale, translate the box, or changing the appearance of the box by assigning a different colour), you can create any sort of "dashboard" output for your simulation....including irregular hexagons.

This technique is a direct application of Data Driven Documents but to a different substrate (than HTML or SVG) here.

Hope this helps.

Answer 3

For example you could take a look of example done by Alexandra Baumgart and Hazuki Okuda using Mathematica. This is implemented with Manipulate effectively creating a basic UI.

Code:

Manipulate[
Grid[{{Show[
ParametricPlot3D[{1.25Cos[t], 1.25 Sin[t],s+2-2w},{s,0,.25},{t,0,2Pi},PlotStyle->Directive[Opacity[1],Gray],Mesh->None, Lighting->"Neutral"], 
ParametricPlot3D[{r Cos[t], r Sin[t],2.25-2w},{r,0,1.25},{t,0,2Pi},PlotStyle->Directive[Opacity[1],Gray],Mesh->None, Lighting->"Neutral"],
ParametricPlot3D[{r Cos[t], r Sin[t],2-2w},{r,0,1.25},{t,0,2Pi},PlotStyle->Directive[Opacity[1],Gray],Mesh->None, Lighting->"Neutral"],
ParametricPlot3D[{1.25Cos[t], 1.25 Sin[t],s-2.25+2w},{s,0,.25},{t,0,2Pi},PlotStyle->Directive[Opacity[1],Gray],Mesh->None, Lighting->"Neutral"],
ParametricPlot3D[{r Cos[t], r Sin[t],-2.25+2w},{r,0,1.25},{t,0,2Pi},PlotStyle->Directive[Opacity[1],Gray],Mesh->None],
ParametricPlot3D[{r Cos[t], r Sin[t],-2+2w},{r,0,1.25},{t,0,2Pi},PlotStyle->Directive[Opacity[1],Gray],Mesh->None, Lighting->"Neutral"],  
ParametricPlot3D[{(1-s/2) Cos[t],(1-s/2) Sin[t],Max[0,-s+2 ]},{s,Min[2-(2^(2/3) (2 \[Pi]-Min[2Pi,3 V((w^2)/(.04))])^(1/3))/\[Pi]^(1/3)-w,1.99],2},{t,0, 2 Pi},PlotStyle->Directive[Opacity[1],Hue[a]],Mesh->None, Lighting->"Neutral"],
ParametricPlot3D[{r Cos[t], r Sin[t],Max[0,(2^(2/3) (2 \[Pi]-Min[2Pi,3 V((w^2)/(.04))])^(1/3))/\[Pi]^(1/3)-w]},{r,0, .000000000001+w+((2^(2/3) (2 \[Pi]-Min[2Pi,3 V((w^2)/(.04))])^(1/3))/\[Pi]^(1/3)-w)/2},{t,0,2Pi},PlotStyle->Directive[Opacity[1], Hue[a]],Mesh->None, Lighting->"Neutral"],
ParametricPlot3D[{r Cos[t], r Sin[t], Min[0,-(2^(2/3) (2 \[Pi]-Min[2Pi,3 V((w^2)/(.04))])^(1/3))/\[Pi]^(1/3)]},{r, 0, .00000000001+w+((2^(2/3) (2 \[Pi]-Min[2Pi,3 V((w^2)/(.04))])^(1/3))/\[Pi]^(1/3)-w)/2},{t,0,2Pi},PlotStyle->Directive[Opacity[1],Hue[a]],Mesh->None, Lighting->"Neutral"],
ParametricPlot3D[{(1-s/2) Cos[t],(1-s/2) Sin[t],Min[0,s -2 ]},{s,w,2},{t,0, 2 Pi},PlotStyle->Directive[Opacity[.2],Gray],Mesh->None, Lighting->"Neutral"],ParametricPlot3D[{(1-s/2) Cos[t],(1-s/2) Sin[t],Max[0,-s + 2 - w]},{s,w,2},{t,0, 2 Pi},PlotStyle->Directive[Opacity[.2],Gray],Mesh->None, Lighting->"Neutral"], 
ParametricPlot3D[{(1-s/2) Cos[t],(1-s/2) Sin[t],Min[0,s-2+ w]},{s,w,Min[2-(2^(2/3) (2 \[Pi]-Min[2Pi,3 V((w^2)/(.04))])^(1/3))/\[Pi]^(1/3)-w,2]},{t,0, 2 Pi},Mesh->None, PlotStyle->Directive[Opacity[1],Hue[a]], Lighting->"Neutral"],  ParametricPlot3D[{(w/2) Cos[t],(w/2) Sin[t], b},{t,0,2Pi}, {b, -2 + w, 0}, PlotStyle->Directive[Opacity[1], Hue[a]],Mesh->None, Lighting->"Neutral"],PlotRange->All,ImageSize->{300,300}, SphericalRegion-> True]},{Row[{Text["time to empty = "], Text[2Pi (.04)/(3w^2)],Text[" seconds"]}]}}],{start,ControlType->None},{end,ControlType->None},
{{V,.01,"time (seconds)"},0.01,34,.01,ControlType->Animator,AnimationRate->1,AnimationRunning->False,ImageSize->Small}, 
{{w,.05,"neck width (millimeters)"}, .05, .3,.01,Appearance->"Labeled"}, 
{{a,0,"color of sand"}, 0, 1,Appearance->"Labeled"}]

Source: http://demonstrations.wolfram.com/FlowTimeInAnHourglass/

Answer 4

Trying to set a standard mesh library, I believe we should start by determining the degrees of freedom.

For concave bissymetric hexagons with 2 mouths it could be:

1. mouth ax width (W_m)
2. relative top width (w_t = W_t / W_m)
3. relative bottom width (w_b = W_b / W_m)
4. relative top height
5. relative bottom height

The answers from Margus and A_A will be helpful for the realization.

Answer 5

If your goal is to visualize such hexagons in matlab, then fill and fill3 should do the trick. Here is an example code, which assumes that your hexagons are parameterized by two widths w1 and w2 and the d1 and d2 parameters which are lengths of the sides:

function draw_hexagon(w1, w2, d1, d2)
    a=(w1-w2)/2;
    b1=sqrt(d1^2 - a^2);
    b2=sqrt(d2^2 - a^2);

    xs=[-w1/2, w1/2, w2/2, w1/2, -w1/2, -w2/2];
    ys=[b1, b1, 0, -b2, -b2, 0];

    fill(xs, ys, 'b')
    axis square 
    grid on
end

It will produce the following for w1=4, w2=2, d1=2, d2=3:

And similarly for 3D:

function draw_hexagon_3d(w1, w2, d1, d2)
    a=(w1-w2)/2;
    b1=sqrt(d1^2 - a^2);
    b2=sqrt(d2^2 - a^2);

    xs=[-w1/2, w1/2, w2/2, w1/2, -w1/2, -w2/2];
    ys=[0, 0, 0, 0, 0, 0];
    zs=[b1, b1, 0, -b2, -b2, 0];

    fill3(xs, ys, zs, 'b')
    grid on
    axis square 
end

we get: