7
$\begingroup$

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:

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.

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:

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.

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.

enter image description here

Even thought the refinement function is quite fine:

enter image description here

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

You could use:

Needs["NDSolve`FEM`"]
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"]

enter image description here

For a larger domain:

Needs["NDSolve`FEM`"]
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}}]]

enter image description here

$\endgroup$
  • $\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$ – tsuresuregusa Jan 5 '16 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$ – tsuresuregusa Jan 5 '16 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 Jan 5 '16 at 18:05
  • $\begingroup$ I hope the image explains what I mean by the disk becoming a triangle. $\endgroup$ – tsuresuregusa Jan 5 '16 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 Jan 5 '16 at 18:13

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.