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]
$\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$ Apr 1 '15 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$
    – user21
    Apr 1 '15 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$
    – unlikely
    Apr 1 '15 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$
    – user21
    Apr 1 '15 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$
    – unlikely
    Apr 4 '15 at 11:29
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.

$\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$
    – unlikely
    Apr 2 '15 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$
    – unlikely
    Apr 2 '15 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$
    – unlikely
    Apr 2 '15 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$
    – user21
    Apr 2 '15 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$
    – user21
    Apr 2 '15 at 11:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.