I need to generate a 2D mesh using triangular elements inside an hexagonal domain. I want my mesh to respect all the symmetries of the hexagon symmetry group. Is there a simple way to do that?
Technically, I suggest you generate the vertices yourself, and use the triangulate command (Mesh Generation) to generate the mesh. However, I have no idea how you should generate the vertices.
What do you hope to gain? There is always the issue of floating point round off
error in these things. Some code seems to take great care of this,
fixed or indefinite precision math, while I’m not sure what others do for tolerances.
If it really is symmetric whatever you are running maybe should know that.
If you are going to solve in all 6 sections don’t you want your result to be fairly
robust against mesh details?
Thank you for your answer. I am solving the Schrodinger equation on hexagonal domains. If the system has inversion symmetry, each of the resulting eigenvalues is doubly degenerate. Nevertheless, if I solve the eignevalue problem on a mesh that is not perfectly symmetric (of course, neglecting the floating point round off errors you mentioned), the degeneracy of the energy eigenvalues is slightly broken.
Thank you for your quick response. Let’s assume I am able to generate the vertices manually as you suggested. I would also need to specify edge labels for the outer boundary and to give different labels to triangles in different regions of the mesh. This is because I want to simulate a core-shell hexagonal structure, with different material properties for the core and the shell region respectively. Would it be possible to add these properties to the mesh object after having generated it with the triangulate command?
I think you can achieve that with the change command. However, I have never used it myself. Consult the documentation and the forum on the change function for more info.
There is a 2D harmonic oscillator example. I just ran that with whatever mesh it makes for
a square. The mesh AFAICT is a perfect grid but I guess if you wanted to play with that you
could run it with a distroted mesh or maybe move a few verticieis and see what happens
to the energies. The first few are below and good to a few sig figs. You could
also play with the element types I guess. Are you looking for something subtle like
a small Jahn-Teller ? In any case you probably want to make sure it is not a function
of numerical parameters 
Degeneracy
of course may be a physical and math problem
Eigenvalue #0 = 1.00104
Eigenvalue #1 = 2.00187
Eigenvalue #2 = 2.00437
Eigenvalue #3 = 3.00332
Eigenvalue #4 = 3.00561
If I had to do that, I would create the elementary triangle (of the “coarse” mesh if you want), then rotate it, and glue it 6 times… This surely won’t have the D_6 symmetry, but I guess if you do it now with equilateral triangles 12 times your mesh should be invariant under the D_6 group.
Otherwise did you look at this thread?
Hope this helps
What about something like this:
real L = 1.0;
border C1(t=0, 1){x=L*t; y=0.0;}
border C2(t=0, 1){x=L*(1.0-t*0.5); y=L*(0.5*t*sqrt(3.0));}
border C3(t=0, 1){x=L*0.5*(1.0-t); y=L*0.5*(1.0-t)*sqrt(3.0);}
int nn = 20;
mesh Th1 = buildmesh(C1(nn) + C2(nn) + C3(nn));
mesh Th2 = movemesh(Th1,[L+(x*cos(pi/3.0)-y*sin(pi/3.0)),x*sin(pi/3.0)+y*cos(pi/3.0)]);
mesh Th3 = movemesh(Th1,[-(x*cos(pi/3.0)-y*sin(pi/3.0)),(x*sin(pi/3.0)+y*cos(pi/3.0))]);
mesh Thhalf1 = Th1+Th2+Th3;
mesh Thhalf2 = movemesh(Thhalf1,[x,L*sqrt(3.0)-y]);
mesh Thwhole = Thhalf1+Thhalf2;
plot(Thwhole);
Hi, thank you for your interest in the topic. Yeah, I agree with you. The most simple way to do that should be to create the “elementary triangle” and then use a combination of the symmetry operations of the group to get the desired mesh invariance. I did not see the thread yuo mentioned, thanks a lot.
Hi, thank you for your answer. Your proposal is inspiring. Probably the correct “elementay triangle” is half the starting equilateral triangle. Otherwise, for some values of nn, the number of edges for each border, we could lose the invariance under reflection with respect some planes, for instance x=0 and y=0 planes. But I’m going to develop your idea, since I think it’s the proper starting point.
Good point. Then simply do this:
real L = 1.0;
border C1(t=0, 1){x=0.5*L*t; y=0.0;}
border C2(t=0, 1){x=0.5*L; y=L*(0.5*t*sqrt(3.0));}
border C3(t=0, 1){x=L*0.5*(1.0-t); y=L*0.5*(1.0-t)*sqrt(3.0);}
int nn = 20;
mesh Th01 = buildmesh(C1(0.5*L*nn) + C2(0.5*L*sqrt(3.0)*nn) + C3(L*nn));
mesh Th02 = movemesh(Th01,[L-x,y]);
mesh Th1 = Th01 + Th02;
mesh Th2 = movemesh(Th1,[L+(x*cos(pi/3.0)-y*sin(pi/3.0)),x*sin(pi/3.0)+y*cos(pi/3.0)]);
mesh Th3 = movemesh(Th1,[-(x*cos(pi/3.0)-y*sin(pi/3.0)),(x*sin(pi/3.0)+y*cos(pi/3.0))]);
mesh Thhalf1 = Th1 + Th2 + Th3;
mesh Thhalf2 = movemesh(Thhalf1,[x,L*sqrt(3.0)-y]);
mesh Thwhole = Thhalf1+Thhalf2;
plot(Thwhole);
Thank you, it seems to work properly. That’s nice! Now I think that I can use the change command to redefine border lables and region labels the way I need.