# How to create union of domains and define functions on subdomains?

**URL:** <https://community.freefem.org/t/how-to-create-union-of-domains-and-define-functions-on-subdomains/1766>\
**Category:** General Discussion\
**Created:** [May 22, 2022, 6:18am UTC](https://community.freefem.org/t/how-to-create-union-of-domains-and-define-functions-on-subdomains/1766 "2022-05-22T06:18:16Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jade](https://avatars.discourse-cdn.com/v4/letter/j/a8b319/32.png) [@Jade](https://community.freefem.org/u/Jade)\
**Post date:** [May 22, 2022, 6:18am UTC](https://community.freefem.org/t/how-to-create-union-of-domains-and-define-functions-on-subdomains/1766/1 "2022-05-22T06:18:16Z")

</div>

I’d like to solve an PDE in a square with a source term g defined piecewisely. The function g(x,y) has one value in an union of 3 disks inside the square and and another value outside the union. I can find the all the intersection points of the circles and define multiple borders, but it will be quite awkward to do so. Besides, I am going to solve the problem with many such problems with 3 disks of random centers and radii. So the number of disjoint boundaries may vary from case to case. I have the following questions.

1. Can we define union of subdomains in FreeFem? If not, I guess I would use external programs such as gmsh to create the mesh, but then I am not sure if it’s easy to identify the internal boundaries since the number of disjoint boundaries vary from one problem to another.

2. Can we define a function that’s only supported in a subdomain? For the piecewise function g, we can define it on the whole domain using logical operators or if-else statements. But it would nicer if we can define separate functions on subdomains.

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [May 22, 2022, 11:49am UTC](https://community.freefem.org/t/how-to-create-union-of-domains-and-define-functions-on-subdomains/1766/2 "2022-05-22T11:49:31Z")

</div>

1. You can use the `+` operator with meshes, read the documentation.
2. You can define multiple `fespace` on different meshes.

---

<div class="post-metadata">

**Author:** ![Jade](https://avatars.discourse-cdn.com/v4/letter/j/a8b319/32.png) [@Jade](https://community.freefem.org/u/Jade)\
**Post date:** [May 22, 2022, 1:17pm UTC](https://community.freefem.org/t/how-to-create-union-of-domains-and-define-functions-on-subdomains/1766/3 "2022-05-22T13:17:32Z")

</div>

I went over the Mesh generation part of the doc. Yes we can add meshes. I tried the following code.

```auto
  border a(t=0, 2*pi){x=cos(t); y=sin(t); label=1;}
  border b(t=0, 2*pi){x=0.5+cos(t); y=sin(t); label=2;}
  mesh Th1 = buildmesh(a(50));
  mesh Th2 = buildmesh(b(50));
  mesh Th = Th1+Th2;
  plot(Th);

```

Then the mesh looks like the following

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/7/79deafdc903d59b6bb57ec8a299a521d2216d931.png)

The meshes overlap! Maybe I missed some other parts of the doc?

If there are only two disks, I can easily define the whole boundary explicitly and then mesh it. However, the situation becomes very complicated if there 3 disks, which may or may not overlap.

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [May 22, 2022, 1:25pm UTC](https://community.freefem.org/t/how-to-create-union-of-domains-and-define-functions-on-subdomains/1766/4 "2022-05-22T13:25:18Z")

</div>

Ah, if the meshes overlap in a non-confirming way, then unless you can define the border explicitly, you cannot use `buildmesh` and the operator `+`.

---

<div class="post-metadata">

**Author:** ![Jade](https://avatars.discourse-cdn.com/v4/letter/j/a8b319/32.png) [@Jade](https://community.freefem.org/u/Jade)\
**Post date:** [May 22, 2022, 1:37pm UTC](https://community.freefem.org/t/how-to-create-union-of-domains-and-define-functions-on-subdomains/1766/5 "2022-05-22T13:37:05Z")

</div>

Thanks! It helps a lot! I guess I will resolve to gmsh then.

---

<div class="post-metadata">

**Author:** ![Jade](https://avatars.discourse-cdn.com/v4/letter/j/a8b319/32.png) [@Jade](https://community.freefem.org/u/Jade)\
**Post date:** [May 24, 2022, 1:02pm UTC](https://community.freefem.org/t/how-to-create-union-of-domains-and-define-functions-on-subdomains/1766/6 "2022-05-24T13:02:37Z")

</div>

I managed to generate the requried mesh using gmsh as in the following code:

```auto
SetFactory("OpenCASCADE");
Rectangle (1) = {-3, -3, 0, 6, 6};
DBC() = CombinedBoundary{Surface{1};}; // Dirichlet BC
Physical Curve(1) = DBC(); // whole boundary, can use index to access different parts
Disk (2) = {0, 0, 0, 0.5};
Disk (3) = {0.8, 0, 0, 0.4};
Disk (4) = {0.3, 0.6, 0, 0.3};
v() = BooleanUnion{ Surface{2}; Delete; }{ Surface{3, 4}; Delete; }; // combine all together with conformal boundaries
w() = BooleanFragments{Surface{1}; Delete; }{Surface{v()}; Delete; };
Physical Surface(1) = w(); // whole domain, can use index to access different parts

Mesh.Algorithm = 6;
Mesh.CharacteristicLengthMin = 0.1;
Mesh.CharacteristicLengthMax = 0.1;

```

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/f/fa28516af59a4182c890b7851d5a1d03ade1e1dc.jpeg)

It has only one physical boundary, which is the boundary of the square. Then I follow this [example](http://jolivet.perso.enseeiht.fr/FreeFem-tutorial/section_8/example14.edp.html) and wrote the following code (loading .mesh or .msh file using `gmshload` failed, so I followed your suggestion to use `DM`):

```auto
load "gmsh"
load "PETSc-complex"
macro dimension()2// EOM
include "macro_ddm.idp"

// import mesh from gmsh
DM dm("3disks.msh");
if(!mpirank)
    mesh Th = dm;
broadcast(processor(0), Th);
plot(Th);

Mat<complex> A;
func Pk = P2;
createMat(Th, A, Pk);

// define the bilinear form
real k = 2.0;
macro grad(u)[dx(u), dy(u)]// EOM
varf vPb(u, v) = int2d(Th)(- (grad(u)' * grad(v)) + k^2 * u * v) + on(1, u=0);
fespace Vh(Th, Pk);
A = vPb(Vh, Vh, tgv=-1);

// define the linear form
int nRhs = 1;
complex[int, int] rhs(Vh.ndof, nRhs), sol(Vh.ndof, nRhs);
for(int i = 0; i < nRhs; ++i) {
    func g = (x^2+y^2<0.5^2) | ((x-0.8)^2+y^2<0.4^2) | ((x-0.3)^2+(y-0.6)^2<0.3^2); // will use multiple rhs later
    varf vRhs(u, v) = int2d(Th)(g * v) + on(1, u=0);
    rhs(:, i) = vRhs(0, Vh, tgv=-1);
}

// solve the system
complex[int, int] B(A.n, rhs.m), X(A.n, rhs.m);
ChangeNumbering(A, rhs, B);
set(A, sparams = "-ksp_type preonly -pc_type lu");
KSPSolve(A, B, X);
ChangeNumbering(A, sol, X, inverse = true, exchange = true);
for(int i = 0; i < rhs.m; ++i) {
    Vh u;
    u[] = sol(:, i).re;
    plotD(Th, u, cmm="u_"+i);
}

```

The code runs well if i use `-np 1`. However, when it goes to multiple processes, the the plot only show parts of the solution. So I am afraid I can’t not output the global solution.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/6/636d153fa42907717ed19dcca6c3fe1b27c2ab23.png)

It looks like a mixing of some of the decomposed domains.

If I use the the built-in mesh such as

```auto
border a(t=0, 2*pi) { x=cos(t); y=sin(t); label=1; }
mesh Th;
if(!mpirank)
    Th = buildmesh(a(300));
broadcast(processor(0), Th);

```

Then it’s ok to use multiple processes. Is there something we should take care for imported mesh when using PETSc?

---

<div class="post-metadata">

**Author:** ![Jade](https://avatars.discourse-cdn.com/v4/letter/j/a8b319/32.png) [@Jade](https://community.freefem.org/u/Jade)\
**Post date:** [May 24, 2022, 1:22pm UTC](https://community.freefem.org/t/how-to-create-union-of-domains-and-define-functions-on-subdomains/1766/7 "2022-05-24T13:22:05Z")

</div>

Problem solved if I use the `msh2` format and use `gmshload`:

```auto
gmsh 3disks.geo -2 -o 3disks.msh -format msh2

```

But the .mesh format always fail to be imported:

```auto
'/Applications/FreeFem++.app/Contents/ff-4.11/ff-petsc/r/bin/mpiexec' -np 8 /Applications/FreeFem++.app/Contents/ff-4.11/bin/FreeFem++-mpi -nw 'randomSource.edp' -ne -wg
-- FreeFem++ v4.11 (Lun 11 avr 2022 22:25:22 CEST - git v4.11)
   file : randomSource.edp
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue parallelempi 
 load: init metis (v 5 )
 sizestack + 1024 =10352 ( 9328 )

mshptg bug in number of points 0 > 0 == max nb points 
 Sorry Error build delaunay triangle error = 1
  current line = 10 mpirank 0 / 8
Exec error : Error mshptg8_
   -- number :1
Exec error : Error mshptg8_
   -- number :1
 err code 8 , mpirank 0

```

In the YouTube video I remembered you suggest use the .mesh rather than the .msh format.

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [May 24, 2022, 1:34pm UTC](https://community.freefem.org/t/how-to-create-union-of-domains-and-define-functions-on-subdomains/1766/8 "2022-05-24T13:34:13Z")

</div>

This is wrong:

```auto
DM dm("3disks.msh");
if(!mpirank)
    mesh Th = dm;
broadcast(processor(0), Th);
plot(Th);

```

You need to do:

```auto
mesh Th;
if(!mpirank) {
    DM dm("3disks.msh", communicator = mpiCommSelf);
    Th = dm;
}
broadcast(processor(0), Th);
plot(Th);

```

---

<div class="post-metadata">

**Author:** ![Jade](https://avatars.discourse-cdn.com/v4/letter/j/a8b319/32.png) [@Jade](https://community.freefem.org/u/Jade)\
**Post date:** [May 24, 2022, 1:46pm UTC](https://community.freefem.org/t/how-to-create-union-of-domains-and-define-functions-on-subdomains/1766/9 "2022-05-24T13:46:45Z")

</div>

This works! Thank you for your prompt reply!

I get a warning:

```auto
Warning: -- Your set of boundary condition is incompatible with the mesh label.

```

However, from my gmsh file I can see there is only one boundary which has label 1. I used your code below and doubled checked the label is 1:

```auto
int[int] lab = labels(Th);
varf onG(u, v) = on(lab, u = label);
fespace Vh(Th, P1);
Vh u;
u[] = onG(0, Vh, tgv = -1);
plot(u, wait=1);

```

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/f/fedd8ad2ddaa90260e5a334a3b49939ecdee24f9.png)

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [May 24, 2022, 1:55pm UTC](https://community.freefem.org/t/how-to-create-union-of-domains-and-define-functions-on-subdomains/1766/10 "2022-05-24T13:55:21Z")

</div>

This warning appears when you have floating subdomains, i.e., subdomains that do not touch the boundary. It’s benign.