Update: 12.1
Another way to generate the mesh is to make use of the OpenCascadeLink. For this we generate a flat cross section of the glass in 3D.
polygon =
Polygon[{{0, 0, 0}, {r2 + del, 0, 0}, {r2 + del, 0, L1}, {r1 + del,
0, L}, {r1, 0, L}, {r2, 0, L1}, {0, 0, L1}}];
Graphics3D[{FaceForm[], EdgeForm[Black], polygon}, Boxed -> False]
We load the link
Needs["OpenCascadeLink`"]
and convert the polygon into an OCCT shape:
shape = OpenCascadeShape[polygon];
We set up an axis of revolution and sweep the polygon.
axis = {{0, 0, 0}, {0, 0, 3/2 L}}; sweep = OpenCascadeShapeRotationalSweep[shape, axis, 2 [Pi]];
Here is a visual of the result:
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[sweep,
"ShapeSurfaceMeshOptions" -> {"LinearDeflection" -> 0.00125}];
Show[Graphics3D[{{Red, polygon}, {Blue, Thick, Arrow[axis]}}],
bmesh["Wireframe"], Boxed -> False]
You see the original polygon in red and the blue arrow is the rotational axis. From here we can generate the mesh in the same way:
mesh = ToElementMesh[bmesh, "MeshOrder" -> 1(*,
"MaxCellMeasure"\[Rule]0.000000005*)]
mesh["Wireframe"[
"MeshElementStyle" ->
Directive[Opacity[0.2], Specularity[White, 17], FaceForm[White],
EdgeForm[]]]]
This is a much better approximation of the geometry. Nevertheless finding the eigenvalues remains challenging as there is a strong dependency of the eigenvalues on the mesh.
Update: 12.1
Another way to generate the mesh is to make use of the OpenCascadeLink. For this we generate a flat cross section of the glass in 3D.
polygon =
Polygon[{{0, 0, 0}, {r2 + del, 0, 0}, {r2 + del, 0, L1}, {r1 + del,
0, L}, {r1, 0, L}, {r2, 0, L1}, {0, 0, L1}}];
Graphics3D[{FaceForm[], EdgeForm[Black], polygon}, Boxed -> False]
We load the link
Needs["OpenCascadeLink`"]
and convert the polygon into an OCCT shape:
shape = OpenCascadeShape[polygon];
We set up an axis of revolution and sweep the polygon.
axis = {{0, 0, 0}, {0, 0, 3/2 L}}; sweep = OpenCascadeShapeRotationalSweep[shape, axis, 2 [Pi]];
Here is a visual of the result:
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[sweep,
"ShapeSurfaceMeshOptions" -> {"LinearDeflection" -> 0.00125}];
Show[Graphics3D[{{Red, polygon}, {Blue, Thick, Arrow[axis]}}],
bmesh["Wireframe"], Boxed -> False]
You see the original polygon in red and the blue arrow is the rotational axis. From here we can generate the mesh in the same way:
mesh = ToElementMesh[bmesh, "MeshOrder" -> 1(*,
"MaxCellMeasure"\[Rule]0.000000005*)]
mesh["Wireframe"[
"MeshElementStyle" ->
Directive[Opacity[0.2], Specularity[White, 17], FaceForm[White],
EdgeForm[]]]]
This is a much better approximation of the geometry. Nevertheless finding the eigenvalues remains challenging as there is a strong dependency of the eigenvalues on the mesh.
You get a better mesh with a different boundary mesh generator:
(mesh = ToElementMesh[reg,
"BoundaryMeshGenerator" -> \
{"BoundaryDiscretizeRegion",
Method -> {"MarchingCubes", PlotPoints -> 33}},
"MeshOrder" -> 1,
"MaxCellMeasure"\[Rule]0.000000005])["Wireframe"]
For that mesh I get
Abs[vals]/(2 Pi)
(*{0.000502385, 0.000502385, 0.00072869, 0.00072869, \
0.000733392, 0.000733392, 0.0010404, 0.0010404, 0.00150767, \
0.00150767, 0.00151325, 0.00151325, 0.308656, 2238.88, 2238.88}*)
And the 14th mode looks like:
MeshRegion[
ElementMeshDeformation[mesh, Re[Through[funs[[14]]["ValuesOnGrid"]]],
"ScalingFactor" -> 10^9]]
ToTwo other comments,: the fact that NDEigensystem gives messages suggests to me that this mesh is still not good enough; as you see I also used MeshOrder->1 as I did not want to wait for a second order mesh to finish. But you might want to try that and a finer mesh. Probably by using more plot points. Perhaps generate the boundary mesh manually?
A second thing that come to mind is that I think you should have some rigid body modes because the glass stands on the table. Maybe experiment with
DirichletCondition[{u[t, x, y, z] == 0, v[t, x, y, z] == 0,
w[t, x, y, z] == 0}, x == 0]
Also, there is a nice Bell Acoustics customer example in the FEMAddOns. You can install that with
ResourceFunction["FEMAddOnsInstall"][]
and find it on the Applications guide page
FEMAddOns/guide/FEMApplications
or have a look at the cloud version of that notebook.
Hope this helps.
You get a better mesh with a different boundary mesh generator:
(mesh = ToElementMesh[reg,
"BoundaryMeshGenerator" -> \
{"BoundaryDiscretizeRegion",
Method -> {"MarchingCubes", PlotPoints -> 33}},
"MeshOrder" -> 1,
"MaxCellMeasure"\[Rule]0.000000005])["Wireframe"]
For that mesh I get
Abs[vals]/(2 Pi)
(*{0.000502385, 0.000502385, 0.00072869, 0.00072869, \
0.000733392, 0.000733392, 0.0010404, 0.0010404, 0.00150767, \
0.00150767, 0.00151325, 0.00151325, 0.308656, 2238.88, 2238.88}*)
And the 14th mode looks like:
MeshRegion[
ElementMeshDeformation[mesh, Re[Through[funs[[14]]["ValuesOnGrid"]]],
"ScalingFactor" -> 10^9]]
To other comments, the fact that NDEigensystem gives messages suggests to me that this mesh is still not good enough; as you see I also used MeshOrder->1 as I did not want to wait for a second order mesh to finish. But you might want to try that and a finer mesh. Probably by using more plot points. Perhaps generate the boundary mesh manually?
A second thing that come to mind is that I think you should have some rigid body modes because the glass stands on the table. Maybe experiment with
DirichletCondition[{u[t, x, y, z] == 0, v[t, x, y, z] == 0,
w[t, x, y, z] == 0}, x == 0]
Also, there is a nice Bell Acoustics customer example in the FEMAddOns. You can install that with
ResourceFunction["FEMAddOnsInstall"][]
and find it on the Applications guide page
FEMAddOns/guide/FEMApplications
or have a look at the cloud version of that notebook.
Hope this helps.
You get a better mesh with a different boundary mesh generator:
(mesh = ToElementMesh[reg,
"BoundaryMeshGenerator" -> \
{"BoundaryDiscretizeRegion",
Method -> {"MarchingCubes", PlotPoints -> 33}},
"MeshOrder" -> 1,
"MaxCellMeasure"\[Rule]0.000000005])["Wireframe"]
For that mesh I get
Abs[vals]/(2 Pi)
(*{0.000502385, 0.000502385, 0.00072869, 0.00072869, \
0.000733392, 0.000733392, 0.0010404, 0.0010404, 0.00150767, \
0.00150767, 0.00151325, 0.00151325, 0.308656, 2238.88, 2238.88}*)
And the 14th mode looks like:
MeshRegion[
ElementMeshDeformation[mesh, Re[Through[funs[[14]]["ValuesOnGrid"]]],
"ScalingFactor" -> 10^9]]
Two other comments: the fact that NDEigensystem gives messages suggests to me that this mesh is still not good enough; as you see I also used MeshOrder->1 as I did not want to wait for a second order mesh to finish. But you might want to try that and a finer mesh. Probably by using more plot points. Perhaps generate the boundary mesh manually?
A second thing that come to mind is that I think you should have some rigid body modes because the glass stands on the table. Maybe experiment with
DirichletCondition[{u[t, x, y, z] == 0, v[t, x, y, z] == 0,
w[t, x, y, z] == 0}, x == 0]
Also, there is a nice Bell Acoustics customer example in the FEMAddOns. You can install that with
ResourceFunction["FEMAddOnsInstall"][]
and find it on the Applications guide page
FEMAddOns/guide/FEMApplications
or have a look at the cloud version of that notebook.
Hope this helps.