1 /*******************************************************
 2  *  This file is part of the FreeFEM tutorial made by  *
 3  *      Pierre Jolivet <pierre@joliv.et>               *
 4  *                                                     *
 5  *  See https://joliv.et/FreeFem-tutorial for more     *
 6  *                                                     *
 7  *   Description: mesh coloring                        *
 8  *******************************************************/
 9 load "example1"
10 border a(t=0, pi) { x=cos(t); y=sin(t); label=1; }
11 border b(t=pi, 2*pi) { x=cos(t); y=sin(t); label=2; }
12 mesh Th = buildmesh(a(10) + b(10));
13 fespace Ph(Th, P0);
14 Ph color;
15 coloring(color[], Th);
16 int ncols = color[].max + 1;
17 {
18     int[int] icolor(color[].n);
19     for(int i = 0; i < color[].n; ++i)
20         icolor[i] = color[][i];
21     Th = change(Th, fregion = icolor[nuTriangle]);
22 }
 
24 int myRegion;
25 varf vPbNeumann(u, v) = int2d(Th, myRegion)(dx(u) * dx(v) + dy(u) * dy(v));
26 fespace Vh(Th, P1);
 
28 for(int i = 0; i < ncols; ++i) {
29     mesh ThCol = trunc(Th, abs(color - i) < 0.1, label = 10);
30     plot(ThCol, wait = 1);
31     Ph ind = abs(color - i) < 0.1 ? 1 : 0;
32     plot(ind, wait = 1, fill = 1);
33     myRegion = i;
34     matrix A = vPbNeumann(Vh, Vh);
35     cout << "Assembled matrix for color #" << (i + 1) << "\n" << A << endl;
36 }
 
38 int[int][int] colors(0);
39 Th = change(Th, flabel = nuTriangle * 3 + nuEdge);
40 plot(Th, wait = 1);
41 coloringBoundary(colors, Th);
42 for(int i = 0; i < colors.n; ++i) {
43     varf onL(u, v) = on(colors[i], u = 1);
44     fespace Wh(Th, P0edge);
45     Wh lab = onL(0, Vh, tgv = -1);
46     plot(lab, wait = 1, fill = 1);
47 }