# Diffusion localised on the map of France

Posted 2 years ago
3517 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 2 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"]
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}] - 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"]
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 ==
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 1 year 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 1 year 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 1 year 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
Posted 2 years ago
 Thanks very much A. Dauphiné