5
$\begingroup$

I want to create an ElementMesh for a disk.

For my specific scenario, some author advise to use a regular mesh, where the radial and tangential directions are divided in equally spaced intervals [1].

Moreover it's probably better to prevent the mesh element measure become too smalls near the center of the disk [2].

So, approaching the center of the disk, I sometimes double the angular step, and, as a first tentative, I built this mesh:

Mathematica graphics [as you can see radial step is nonconstant but this is not relevant for the rest]

I tried to use this mesh as domain for a simple Poisson equation solved by NDSolve. But when the resulting InterpolatingFunction is passed to Plot3D the kernel crash without giving any reason.

I suppose this is related to the fact my first tentative mesh doens't cover the disk whole domain and/or is non conformal.

What is a good-quality mesh with the previous requirements and supported by NDSolve, FEM and Plot functions?

UPDATE

[1] To be more precise the radial step has to be piecewise-constant, so that mesh elements doen't cross some internal, circular, boundary.

[2] The autor apparently advise to use a non constant angular step. I chosen to double the angular step approaching to the center for ease.

ABOUT THE KERNEL CRASH

The code to generate the previous mesh was a bit involved to take into account the specific needs for a piecewise constant radial step. I'm working on a code to generate a mesh without "gaps" so I don't have anymore the previous code. But, to reproduce the kernel crash:

mesh = Import["http://1drv.ms/1EJwm8D"] NDSolve[{Laplacian[u[x, y], {x, y}] == -1, DirichletCondition[u[x, y] == 0, Norm[{x, y}] >= 6.31436*10^6]}, u, {x, y} \[Element] mesh] Plot3D[u[x, y] /. First[%], {x, y} \[Element] mesh, PlotRange -> All]
Improve this question
$\endgroup$
10
  • $\begingroup$ Each row before the tangential density doubles, the mesh elements need to be pentagons, not trapezoids, in order to fill the small triangular space left at the top. $\endgroup$ Commented Apr 1, 2015 at 18:36
  • $\begingroup$ Does this mesh have hanging nodes (T-junctions)? Then it depends how the mesh is build. In any case if you see a crash it were good if you could send the data $\endgroup$ Commented Apr 1, 2015 at 18:59
  • $\begingroup$ @user21 Indeed I don't know how to "declare" hanging nodes... But at present I don't even have an idea on how to design a good mesh, with my requirements, and how to mesh the area where the angular density doubles. $\endgroup$ Commented Apr 1, 2015 at 19:07
  • $\begingroup$ @unlikely, I am not sure I understand do you have the above mesh or are you trying to create that (and in that case, what is crashing)? $\endgroup$ Commented Apr 1, 2015 at 19:17
  • 1
    $\begingroup$ @2012rcampion Pentagons are not available in an ElementMesh so I'm working on using 3 TriangleElement instead of 1 "trapezoid"... $\endgroup$ Commented Apr 4, 2015 at 11:29

1 Answer 1

Reset to default
7
$\begingroup$

Here is a way:

Needs["NDSolve`FEM`"] start = E^(1/4) // N; stop = E^2 // N; nx = 40; ny = 20; coordinates = Flatten[ Table[{r Cos[\[Theta]], r Sin[\[Theta]]}, {r, Table[Log[x], {x, start, stop, (stop - start)/(ny - 1)}]}, {\[Theta], 0., 2. Pi, (2. Pi - 0.)/(nx - 1)}], 1]; incidents = Flatten[Table[{j*nx + i, j*nx + i + 1, (j - 1)*nx + i + 1, (j - 1)*nx + i}, {i, 1, nx - 1}, {j, 1, ny - 1}], 1]; lc = Max[incidents] + 1; mesh = ToElementMesh["Coordinates" -> Join[coordinates, {{0, 0}}], "MeshElements" -> {QuadElement[incidents], TriangleElement[ Join[{lc}, #] & /@ Partition[Range[nx - 1], 2, 1, 1]]}];

This will create a mostly quad mesh. The center to the first row of quads is done with triangle elements. The power of start will set the inner circle. For 

start = E^(1/2) // N; stop = E^2 // N; nx = 10; ny = 3;

you get

enter image description here

For

start = E^(1/4) // N; stop = E^2 // N; nx = 40; ny = 20;

You get

enter image description here

You could inspect mesh["Quality"] and use MeshOrderAlteration[mesh, 2] to get a second order mesh. But then it might make sense to move the mid side nodes of the outer quad elements onto the circle.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Although some reference apparently recommend tu use a piecewise constant radial step, your coice is interesting: can you add few words about it? $\endgroup$ Commented Apr 2, 2015 at 10:33
  • $\begingroup$ As of my question the same reference apparently recommend to use a non constant tangential density (instead of a non constant radial density) to control mesh element measure and this is what I was trying to do. My problem is how to handle the transition from one tangential density to the double tangential density. I'm working on this. $\endgroup$ Commented Apr 2, 2015 at 10:36
  • $\begingroup$ The use of MeshOrderAlteration[mesh, 2]["MeshOrder"] (documented somewhere?) to change the mesh order is definitely interesting. I'll experiment with it once done with the first-order mesh. $\endgroup$ Commented Apr 2, 2015 at 10:38
  • $\begingroup$ @unlikely, In principal a mesh with hanging nodes is what you want. In this version (V10.1) a hanging node mesh is most likely only going to work for stationary problems, if at all. I think the better way to proceed here is to use triangle elements to transition from a row of quad elements to a next row of quad elements, but that's not super easy. As to my approach, it looks nice, but if you check the mesh quality, you'll see that the mean is not too good. I think with the approach you have in mind you might get a better quality mesh. If you come up with a solution. I'd be interested to see it $\endgroup$ Commented Apr 2, 2015 at 11:03
  • $\begingroup$ @unlikely, no MeshOrderAlteration is not documented. But it's easy. first the mesh, then the order you want. Possible choices 1 or 2. $\endgroup$ Commented Apr 2, 2015 at 11:04

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.