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
verbosity=0;
// 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.