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[ 
    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:

enter image description here

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.


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:

zoom in for the region user21 proposes but with a 2 outer radius

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.


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!


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

enter image description here

Even thought the refinement function is quite fine:

enter image description here

  • $\begingroup$ You might try constructing your own boundary mesh. Code will look like ToElementMesh[ ToBoundaryMesh["Coordinates" -> ..., "BoundaryElements" -> { ...}]] $\endgroup$
    – george2079
    Commented Jan 5, 2016 at 16:12
  • $\begingroup$ 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? $\endgroup$
    – user21
    Commented Jan 5, 2016 at 16:26
  • $\begingroup$ Sorry, omega is the region where do the meshing, omega 3 was a typo. $\endgroup$ Commented Jan 5, 2016 at 17:45

1 Answer 1


You could use:

   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"]

enter image description here

For a larger domain:

mesh = ToElementMesh[
    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}}]]

enter image description here

  • $\begingroup$ 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. $\endgroup$ Commented Jan 5, 2016 at 17:54
  • $\begingroup$ for example, if I run your code with a outer radius of 2 the central half disk becomes a triangle, very resolved. $\endgroup$ Commented Jan 5, 2016 at 17:58
  • $\begingroup$ @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. $\endgroup$
    – user21
    Commented Jan 5, 2016 at 18:05
  • $\begingroup$ I hope the image explains what I mean by the disk becoming a triangle. $\endgroup$ Commented Jan 5, 2016 at 18:08
  • $\begingroup$ @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? $\endgroup$
    – user21
    Commented Jan 5, 2016 at 18:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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