How to Access Appropriate Triangles in the Mesh

Hi everyone!

I am not sure if I am implementing the correct code for this problem:

For my case, denote g1=\gamma, gdiv =\gamma. The domain \Omega lies inside the rectangular domain O, where \Gamma:=\partial\Omega. Denote T be the triangles inside the mesh and \mathcal{T}_h^O, be a quasi uniform mesh consisting of triangles of max diameter h on O and consider also the following :

  • \mathcal{T}_h = \{T\in\mathcal{T}_h^O: T\cap \Omega\neq \emptyset \}
  • \Omega_h = (\cup_{T\in\mathcal{T}_h}T), where its boundary is $Gamma_h$
  • \mathcal{T}_h^\Gamma = \{T\in\mathcal{T}_h: T\cap \Gamma = \emptyset \}
  • \Omega_h^\Gamma = (\cup_{T\in\mathcal{T}_h^\Gamma}T), where its internal boundary is \Gamma_h^i

My questions are:

  • Is it possible to access the triangles defined on \Omega_h, \Omega_h^\Gamma, \Gamma_h and \Gamma_h^i?
  • Which set satisfies the FreeFem command int1d(Th, levelset=phi)/int2d(Th, levelset=phi), where \Omega is the zero level set of \phi?

I really need to differentiate the integrals on these sets. I attempted to code the above variational formulation.

varf va1 (u, v) = int2d(Th,levelset=phi)(((dx(u) * dx(v)) + (dy(u)*dy(v))))
                    + int1d(Th,levelset=phi)(g1*((dx(u)*dx(v)) + (dy(u)*dy(v))))
                    + intalledges(Th)(delta*diam*(jump(dn(u))*jump(dn(v))));
  varf va2 (yh1, zh1) = int1d(Th, levelset=phi)(gdiv*dx(yh1)*dx(zh1))
                      + int1d(Th,levelset=phi)(g1*yh1*zh1);
  varf va3 (yh2, zh2) = int1d(Th, levelset=phi)(gdiv*dy(yh2)*dy(zh2))
                      + int1d(Th,levelset=phi)(g1*yh2*zh2);
  varf va4 (yh1, zh2) = int1d(Th, levelset=phi)(gdiv*dx(yh1)*dy(zh2));
  varf va5 (yh2, zh1) = int1d(Th, levelset=phi)(gdiv*dy(yh2)*dx(zh1));
  varf va6 (yh1, v) = int1d(Th, levelset=phi)(g1*yh1*dx(v));
  varf va7 (yh2, v) = int1d(Th, levelset=phi)(g1*yh2*dy(v));
  varf va8 (u, zh1) = int1d(Th, levelset=phi)(g1*dx(v)*zh1);
  varf va9 (u, zh2) = int1d(Th, levelset=phi)(g1*dy(v)*zh2);
  varf vL1 (u, v) = int2d(Th,levelset=phi)(f*v)
                   + int1d(Th,levelset=phi)(g*v);
  varf vL2 (yh1, zh1) = int1d(Th,levelset=phi)(gdiv*f*dx(zh1));
  varf vL3 (yh1, zh2) = int1d(Th,levelset=phi)(gdiv*f*dy(zh2));
  varf vb (u, v) = int1d(Th,levelset=phi)(1.*v);

  matrix A1 = va1(Vh, Vh);
  matrix A2 = va2(Vh, Vh);
  matrix A3 = va3(Vh, Vh);
  matrix A4 = va4(Vh, Vh);
  matrix A5 = va5(Vh, Vh);
  matrix A6 = va6(Vh, Vh);
  matrix A7 = va7(Vh, Vh);
  matrix A8 = va8(Vh, Vh);
  matrix A9 = va9(Vh, Vh);
  real[int] b1 = vL1(0, Vh);
  real[int] b2 = vL2(0, Vh);
  real[int] b3 = vL3(0, Vh);

  real[int] bb(3*p+1), xx(3*p+1), l(2*p+1);
  bb = [b1, b2, b3, 0.0];

  real[int] bd = vb(0, Vh);
  real[int] B(3*p);
  B(0:p-1)=bd;// rhs part for ut1
  B(p:3*p-1)=0.0;     // rhs part for pt
  // Block matrix
  matrix A = [[A1,A6,A7], [A8,A2,A5], [A9,A4,A3]];
  matrix<real> AA = [ [ A, B ], [ B', 0.0 ] ];
  set(AA, solver=CG);
  // Solve and set values
  xx = AA^-1 * bb;
  [u[],l] = xx;

I am not sure if this is the correct approximation of the solution u.
Thank you.

Hello, I don’t know if it will solve your problem, but to differenciate each triangle on my mesh, I do the following:

int globalH = 2;
mesh P = square(globalH,globalH);

int NbTriangles = P.nt; //number of triangles

fespace V0(P,P0);
V0 u0;

//---------------------------get triangle numbering
real[int] elemGlobal(NbTriangles);
for(int i = 0; i < NbTriangles; i++){
 u0[][i] = i;