Skip to main content

finite element method - 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:


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



Answer



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


Comments