Greetings FreeFem community,
I am quite new to FreeFem/mmg tools and I’m trying to set-up a diffusion equation on a mesh with multiple regions (generated via supplying a level-set function to mmg2d)
What I want to do is “carve out” an initial sub-region via mmg and force the solution to be constant on one side of the sub-region and be integrated “normally” on the other side, with a Neumann boundary separating them both.
So as a first step, I managed to get mmg to remesh based on an implicit curve (circle) and successfully obtained a mesh with two regions:
Next, I set-up the diffusion problem initial condition, which looks something like this:
And then I set up the problem itself and integrate it with steps dt until Tmax is reached. What I’d expect with a zero-flux boundary is for the outer region to “diffuse” and the inner region to stay at the same constant value.
However, this is not what I observe:
i.imgur.com/W9WQlqQ.png
Here, obviously the isolines are passing through the circular boundary, i.e. the zero-flux condition was not respected.
For the diffusion equation Neuman b.c.s are natural, so there isn’t an operator like on() to set them explicitly, so the problem is elsewhere.
Maybe I am misunderstanding mmg2d’s resulting mesh, and elements with label = 10 is not understood by FreeFem as a boundary? ref: https://forum.mmgtools.org/t/mmg3d-how-to-identify-the-boundary-label/380
So this is where I need help - how do I set-up this problem properly?
Here’s a simplified script that illustrates my current attempts to approach this:
load "distance"
load "mmg"
load "ffrandom"
srandomdev();
int Nelements = 150;
// label references: https://forum.mmgtools.org/t/mmg3d-how-to-identify-the-boundary-label/380
int levelSetBoundaryLabel = 10;
int insideLevelSetRegionLabel = 3;
int outsideLevelSetRegionLabel = 2;
real dt = 0.1;
real Tmax = 5;
real D = 0.1;
int cIni = 5;
real cEq = 0.1;
real xA, yA, R, hmin;
xA = 0.5;
yA = 0.5;
R = 0.1;
hmin = 0.001;
func real ic(real xi, real eta) {
if ((xi - xA) ^ 2 + (eta - yA) ^ 2 <= R ^ 2) {
return cEq; // constant inside circle
} else {
return cIni * randreal1(); // zero outside
}
}
func c0 = ic(x, y);
mesh Th = square(Nelements, Nelements);
fespace Vh(Th, P1);
Vh levelSet, signedDistance;
Vh cOld, v, c = c0;
levelSet = (x - xA) ^ 2 + (y - yA) ^ 2 - R ^ 2;
distance(Th, levelSet, signedDistance[]);
Th = mmg2d(Th, iso = 1, ls = 0.0, metric = signedDistance[], hmin = hmin);
problem diffusion(c, v) =
int2d(Th)(
c * v / dt +
D * (
dx(c) * dx(v) +
dy(c) * dy(v)
)
) -
int2d(Th)(
cOld * v / dt
);
plot(Th, wait = true); // The generated mesh from the level set function
plot(c, fill = true, value = true, wait = true); // the initial scalar field
for (real t = 0; t < Tmax; t += dt) {
cOld = c;
diffusion;
}
plot(c, value = true, wait = 1); // The result after integration
P.S.: I understand that it would be much easier if I just generated the circle with a parametrized border, but this is part of a level-set problem I’m trying to set-up.
P.S.2: It seems that as a new user I’m not allowed to have more than one image directly in the post and two links so I had to break some of the links on purpose. Sorry about that
As you are trying to get J.N to zero you