# how to create a FEM mesh with "RegionPlot" where the sample points are a function of the position

I need to resolve the diffusion equation in a domain with circular sources (holes).

I've tried defining the whole region with holes but the result is very mesh-dependent and not symmetric, so at the moment I'm imposing the symmetry myself by only solving 1/4th of the system.

cords = Table[{i, 0}, {i, {0, 5}}];
Ω = Apply[And, Norm[{x, y} - #]^2 > 1 & /@ cords]

mesh2 = ToElementMesh[
ImplicitRegion[
And[x > 0 &&
y > 0 && ( x - 2)^2 + y^2 <= 900, Ω], {x, y}],
"MaxBoundaryCellMeasure" -> .5, "ImproveBoundaryPosition" -> False,
"MaxCellMeasure" -> 10,
"BoundaryMeshGenerator" -> {"RegionPlot",
"SamplePoints" -> 300}];


This looks like this:

Having so many sample points makes the mesh creation really slow, so I would like to have a way to make the resolution position-dependent, having more resolution on the region plot near the centre and less in the outer border.

EDIT:

The problem I have is that I need to go to relatively large regions, and if I use the mesh refining option as user21 recommends, I end up with something like this:

A coarse approximation to the solution is two logarithmic decays from the centre of each particle, so I'd like the resolution of the boundary to go like that.

EDIT2:

somebody should have told me, "don't try FEM in V10.0", most of the other problems I was having disappeared after updating.

Not only the MeshRefinementFunction works as expected in 10.3 but also the solution is much more symmetric and continuous for the cases I'm looking at. I don't think it's only because the mesh is better made in the new version, seems like a more robust solver, at least from the user side of it.

Despite "MaxBoundaryCellMeasure" not accepting a function of position, the MeshRefinementFunction can be used to refine wherever one wants to. In my case I'm looking at something like this:

MeshRefinementFunction ->
Function[{vertices, area},
area > 0.0125 (0.1 +
If[Norm[Mean[vertices] - {2.5, 0}] < 5, 4,
4 + Norm[Mean[vertices] - {2.5, 0}]^2])


Which is pretty much exactly what I wanted to achieve when i posed the question. Thanks user21!

Update

The problem seems to persist in mma 11, at least in my mac.

Even thought the refinement function is quite fine:

• You might try constructing your own boundary mesh. Code will look like ToElementMesh[ ToBoundaryMesh["Coordinates" -> ..., "BoundaryElements" -> { ...}]] Jan 5, 2016 at 16:12
• There are a few missing pieces of information - what is Ω3, and what is Ω needed for? Why does the boundary mesh generator need to be region plot? Jan 5, 2016 at 16:26
• Sorry, omega is the region where do the meshing, omega 3 was a typo. Jan 5, 2016 at 17:45

You could use:

Needs["NDSolveFEM"]
ToElementMesh[
RegionDifference[
RegionDifference[Disk[{0, 0}, 1, {0, \[Pi]/2}],
Disk[{0, 0}, 1/25]], Disk[{3/10, 0}, 1/25]],
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.0005 (0.1 + 2 Norm[Mean[vertices]])]]["Wireframe"]


For a larger domain:

Needs["NDSolveFEM"]
mesh = ToElementMesh[
RegionDifference[
RegionDifference[Disk[{0, 0}, 5, {0, \[Pi]/2}],
Disk[{0, 0}, 1/25]], Disk[{3/10, 0}, 1/25]],
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.0005 (0.1 + Norm[Mean[vertices]])]];
mesh["Wireframe"[PlotRange -> {{-0.1, 1}, {-0.1, 1}}]]


• The problem with MeshRefinement is that it does precisely what it say, it refines an already existing mesh. The problem is that as I make the outer circle larger, the resolution of the smaller half circles diminish, so refining the mesh resulted in a very refined square in the worst case I tried, let me see if I can recover the problem. I guess what I'm looking for is something like a "MaxBoundaryCellMeasure" that accepts a function of the position instead of a single value. Jan 5, 2016 at 17:54
• for example, if I run your code with a outer radius of 2 the central half disk becomes a triangle, very resolved. Jan 5, 2016 at 17:58
• @tsuresuregusa, What do you mean with "the central half disk becomes a triangle"? The MeshRefinementFunction does the refinement during the mesh creation, there is no initial mesh that is refined. There is no way to have "MaxBoundaryCellMeasure" accept a function of position. Jan 5, 2016 at 18:05
• I hope the image explains what I mean by the disk becoming a triangle. Jan 5, 2016 at 18:08
• @tsuresuregusa, that's because you use the region plot boundary mesh generator. Use the default one, that works better. Possibly you'd need to get a newer version of Mathematica? Jan 5, 2016 at 18:13