# Solving diffusion equation on mesh with multiple sub-regions

**URL:** <https://community.freefem.org/t/solving-diffusion-equation-on-mesh-with-multiple-sub-regions/2954>\
**Category:** General Discussion\
**Created:** [February 9, 2024, 11:24pm UTC](https://community.freefem.org/t/solving-diffusion-equation-on-mesh-with-multiple-sub-regions/2954 "2024-02-09T23:24:04Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![vasilvas](https://avatars.discourse-cdn.com/v4/letter/v/3d9bf3/32.png) [@vasilvas](https://community.freefem.org/u/vasilvas)\
**Post date:** [February 9, 2024, 11:24pm UTC](https://community.freefem.org/t/solving-diffusion-equation-on-mesh-with-multiple-sub-regions/2954/1 "2024-02-09T23:24:04Z")

</div>

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:

 ![mesh level set](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/9/9ad2ae80ff33b3174e3275b62623815f6dc181cc.png)

Next, I set-up the diffusion problem initial condition, which looks something like this:

[![](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/5/5a910cedad4b43a58e5207295c4426f54d2ea672.png) ](https://i.imgur.com/u7SfKTF.png)

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:

```auto
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:

```auto
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 ☹

---

<div class="post-metadata">

**Author:** ![vasilvas](https://avatars.discourse-cdn.com/v4/letter/v/3d9bf3/32.png) [@vasilvas](https://community.freefem.org/u/vasilvas)\
**Post date:** [February 9, 2024, 11:26pm UTC](https://community.freefem.org/t/solving-diffusion-equation-on-mesh-with-multiple-sub-regions/2954/2 "2024-02-09T23:26:27Z")

</div>

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

 ![isolines](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/d/dd79f2b527e9ed8cc1a582c78cae3885c3bc744d.png)

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [February 10, 2024, 11:36am UTC](https://community.freefem.org/t/solving-diffusion-equation-on-mesh-with-multiple-sub-regions/2954/3 "2024-02-10T11:36:05Z")

</div>

Its only natural when the mesh ends 🙂 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 🙂

[https://community.freefem.org/t/normal-boundary-condition/341](https://community.freefem.org/t/normal-boundary-condition/341)

---

<div class="post-metadata">

**Author:** ![vasilvas](https://avatars.discourse-cdn.com/v4/letter/v/3d9bf3/32.png) [@vasilvas](https://community.freefem.org/u/vasilvas)\
**Post date:** [February 10, 2024, 12:50pm UTC](https://community.freefem.org/t/solving-diffusion-equation-on-mesh-with-multiple-sub-regions/2954/4 "2024-02-10T12:50:53Z")

</div>

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 🙂

---

<div class="post-metadata">

**Author:** ![vasilvas](https://avatars.discourse-cdn.com/v4/letter/v/3d9bf3/32.png) [@vasilvas](https://community.freefem.org/u/vasilvas)\
**Post date:** [February 10, 2024, 5:08pm UTC](https://community.freefem.org/t/solving-diffusion-equation-on-mesh-with-multiple-sub-regions/2954/5 "2024-02-10T17:08:46Z")

</div>

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

```auto
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:

```auto
func real diffCoeff(real xi, real eta, mesh M) {
// https://doc.freefem.org/references/global-variables.html#label
  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:

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

```

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [February 10, 2024, 6:11pm UTC](https://community.freefem.org/t/solving-diffusion-equation-on-mesh-with-multiple-sub-regions/2954/6 "2024-02-10T18:11:41Z")

</div>

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.

```auto
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);

```auto

```

---

<div class="post-metadata">

**Author:** ![vasilvas](https://avatars.discourse-cdn.com/v4/letter/v/3d9bf3/32.png) [@vasilvas](https://community.freefem.org/u/vasilvas)\
**Post date:** [February 10, 2024, 10:18pm UTC](https://community.freefem.org/t/solving-diffusion-equation-on-mesh-with-multiple-sub-regions/2954/7 "2024-02-10T22:18:08Z")

</div>

I see, thanks a lot for the support!

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