# Diffusion localised on the map of France

Posted 6 years ago
9747 Views
|
7 Replies
|
13 Total Likes
|
 I have this program to simulate diffusion on the map of France ClearAll["Global*"] Needs["NDSolveFEM"] carto = DiscretizeGraphics[CountryData["France", {"Polygon", "Mercator"}]]  bmesh = ToBoundaryMesh[carto, "MaxBoundaryCellMeasure" -> 25, AccuracyGoal -> 1]; mesh = ToElementMesh[bmesh, MaxCellMeasure -> 5, "MaxBoundaryCellMeasure" -> 25]; mesh["Wireframe"]  op = -Laplacian[u[x, y], {x, y}] - 20; usol = NDSolveValue[{op == 1, DirichletCondition[u[x, y] == 0, True]},u, {x, y} \[Element] mesh]; Plot3D[usol[x, y], {x, y} \[Element] mesh, PlotTheme -> "Detailed", ColorFunction -> "Rainbow", PlotPoints -> 50]  Show[ContourPlot[usol[x, y], {x, y} \[Element] mesh, ColorFunction -> "Temperature"], bmesh["Wireframe"]] I obtained an image of the diffusion. But I want that the diffusion process is centered upon Paris. Do you have a solution ?Thanks! ~ André Dauphiné
7 Replies
Sort By:
Posted 6 years ago

I don't know if this solution is what you are going for - it is more of a kludge, and I can't say if it correctly models diffusion from a center to a regions boundaries.

## Method 1 - add an interior boundary and second DirichletCondition

But this is a quick method. First you use the "IncludePoints" option of ToBoundaryMesh to include a set of points describing a circular region around Paris. These points are then considered part of the boundary for NDSolve. Now you add a second DirichletCondition setting the value of u to some high constant at that interior boundary.

The two adjustable parameters here are the radius around Paris that you consider, and the value you set the function to at that boundary.

ClearAll["Global*"]
Needs["NDSolveFEM"]
parisradius = 0.3;
parisdiffusionvalue = 150;
carto = DiscretizeGraphics[
CountryData["France", {"Polygon", "Mercator"}]];
paris = First@
GeoGridPosition[
GeoPosition[
CityData[Entity["City", {"Paris", "IleDeFrance", "France"}]][
"Coordinates"]], "Mercator"];
bmesh = ToBoundaryMesh[carto, AccuracyGoal -> 1,
"IncludePoints" -> CirclePoints[paris, parisradius, 50]];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 5];
mesh["Wireframe"];
op = -Laplacian[u[x, y], {x, y}] - 20;
usol = NDSolveValue[{op == 1,
DirichletCondition[u[x, y] == 0, Norm[{x, y} - paris] > .6],
DirichletCondition[u[x, y] == parisdiffusionvalue,
Norm[{x, y} - paris] < .6]}, u, {x, y} \[Element] mesh];
Show[ContourPlot[usol[x, y], {x, y} \[Element] mesh,
PlotPoints -> 100, ColorFunction -> "Temperature"],
bmesh["Wireframe"]] // Quiet
Plot3D[usol[x, y], {x, y} \[Element] mesh, PlotTheme -> "Detailed",
ColorFunction -> "Rainbow", PlotPoints -> 50]


Thanks to kglr for showing how to convert city coordinates to the Mercator projection.

## Method 2 - Modify the differential equations

This seems to give a better result, but again it isn't based on any transport theory. All I'm doing is setting the Laplacian equal to a Gaussian function centered around Paris. The only important adjustable parameter is the width of the Gaussian,

ClearAll["Global*"]
Needs["NDSolveFEM"]
parisradius = 2;
parisdiffusionvalue = 150;
carto = DiscretizeGraphics[
CountryData["France", {"Polygon", "Mercator"}]];
paris = First@
GeoGridPosition[
GeoPosition[
CityData[Entity["City", {"Paris", "IleDeFrance", "France"}]][
"Coordinates"]], "Mercator"];
bmesh = ToBoundaryMesh[carto, AccuracyGoal -> 1];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 5];
mesh["Wireframe"];
op = -Laplacian[u[x, y], {x, y}];
usol = NDSolveValue[{op ==
150 Exp[-Norm[{x, y} - paris]^2/parisradius],
DirichletCondition[u[x, y] == 0, True]},
u, {x, y} \[Element] mesh];
Show[ContourPlot[usol[x, y], {x, y} \[Element] mesh,
PlotPoints -> 100, ColorFunction -> "Temperature"],
bmesh["Wireframe"]] // Quiet
Plot3D[usol[x, y], {x, y} \[Element] mesh, PlotTheme -> "Detailed",
ColorFunction -> "Rainbow", PlotPoints -> 50]


Posted 6 years ago
 Thanks very much A. Dauphiné
Posted 6 years ago
 I hope it helps, I am new to using finite element mesh. Is there a source I could look at for the diffusion equation you are modeling? Probably there is a better way to model a single source in Paris than the Gaussian used in the second plot.
Posted 6 years ago
 Yes, I would also be interested in what kind of general problem @André Dauphiné is trying to approach. André, could you please share?
Posted 5 years ago
 @Jason which version of Mathematica did you use to make these. I tried 10.4 and 11.1 with fresh kernel and i am not able to run your code. it is giving an error. Can you please clarify?
Posted 5 years ago
 Ali - I don't know what happened here, I was most certainly using version 10.xx - not sure which, but when I evaluate this now it doesn't work. The problem is the one mentioned by the OP here: http://community.wolfram.com/groups/-/m/t/1023182There is a bug with DiscretizeGraphics where it isn't working on multi-polygons. Make this change in the definition of carto in either of the code blocks above and it should work carto = DiscretizeGraphics[ CountryData["France", {"Polygon", "Mercator"}] /. Polygon[x : {{{_, _} ..} ..}] :> Polygon /@ x]; `
Posted 5 years ago
 Thanks this works ! On a different note, I do also wish that people at Wolfram can incorporate some ways to repair broken meshes/or self intersections in mesh
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments