Solving diffusion equation on mesh with multiple sub-regions

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:

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:

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"


int Nelements = 150;

// label references:
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) =
    c * v / dt +
    D * (
      dx(c) * dx(v) +
      dy(c) * dy(v)
  ) -
    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;

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 :frowning:

Posting the broken on-purpose image here to save you some copy and pasting

Its only natural when the mesh ends :slight_smile: As you are trying to get J.N to zero you
could also set D to zero [edit lol ] there. Preventing flow across the border just applies to
normal component.
See if this helps, I was going to test it first but apparently my FF doesn’t support
mmg2d and I didn’t have time to change it :slight_smile:

Totally went over my head that it isn’t as direct when the mesh is connected there, seems obvious in hindsight.

Thanks a lot for pointing it out and the reference! I’ll play around with the ideas you’ve given me :slight_smile:

@marchywka the penalty method seems to work, but leads to some funny instabilities to this rather simple diffusion problem.

To this end, I decided to try and define a discontinuous diffusion coefficient, but I got stuck on indexing/modifying dofs.

So say I have:

Vh dInterp = D;  // interpolated diffusion coefficient in the FE-space

Is there a way to select all indices of the dInterp[] array that correspond to elements with a specific label and set them to 0?

I also tried to define the following function:

func real diffCoeff(real xi, real eta, mesh M) {
  if (M(xi, eta).label == LABEL) {
    return 0;
  } else {
    return D;

and interpolate that into the fespace, but I get the following strange compile-timetype error:

   44 :   if (M(xi,eta) error operator   <N5Fem2D4MeshE>, <d>, <d> 
 List of choices 
         (        <N12_GLOBAL__N_18lgVertexE> :   <N5Fem2D4MeshE>, <l> )

The code I use is something like this ( mutatis mutandi )
but there are ways to define regions. Note that in general you need to check when
your create a new spatial variable. In this case IIRC it drops from the weak form
( rather puzzling lol ) but you may want to put it into continuity eqn with Fick’s law etc.

func real  fTf(real xx,real yy)
//real f=.1;
real xz=xx/szx; // (xx-szx*.5)/(1.0*szx);
real yz=yy/szy; //(yy-szy*.5)/(1.0*szy);
if ( yz>f ) return 0; 
if ( yz< -f  ) return 0;
if ( xz>f ) return 0; 
if ( xz< -f  ) return 0;
return 1;
} // fTf 

//Vh Tf=Th(x,y).region; // fTf(x,y);
Vh Tf= fTf(x,y);

1 Like

I see, thanks a lot for the support!

Turns out, there are quite a few fine details in working with FF :slight_smile: