Hello,
I have a problem with two meshes, say Th1 and Th2, with Th1 extracted from Th2 (Th1 subset of Th2). I have to integrate in Th1 a product u1*u2, with u1 in Th1, u2 in Th2, and u1 and u2 say of class P1.
I prefer to make the integration “by hand” (not using int2d(Th1)), and for this I need the values of u2 on the nodes of Th1.
I thought I could call u2 as u2[Th1[nT1][jj]], with nT1 triangle in Th1, but it seems that this call refers to different values of u2.
So, it seems that the numbering of triangles and nodes/triangle of Th2 and Th1 is not the same in the intersection. I have attached a simple code which demonstrate this.
My question is how one can call the right values of u2 in Th1, knowing the triangle index nT1 and jj in Th1 .
Thank you!
load “Curvature”
int lb11=11,lb12=12,lb21=21;
// borders
border b11(t=0,2) {
if (t<=1) {x=1; y=t;}
else {x=2-t; y=1.;}
label=lb11;}
border b12(t=0,2) {
if (t<=1) {x=0; y=1.-t;}
else {x=t-1.; y=0;}
label=lb12;}
border b21(t=4,10) {
if (t<=5) {x=t-3; y=0;}
else if (t<=7) {x=2.; y=t-5.;}
else if (t<=9) {x=9-t; y=2.;}
else {x=0;y=11-t;}
label=lb21;}
//plot(b11(2)+b12(2),wait=1,dim=2);
//meshes
//plot(b21(6),dim=2,wait=1);
mesh Th2 = buildmesh(b11(2)+b12(2)+b21(6));
//plot(Th2,wait=1);
real[int,int] XYb11(3,1);
real lengthb11=extractborder(Th2,lb11,XYb11);
border newb11(t=0,XYb11.m-1){ P.x=XYb11(0,t); P.y=XYb11(1,t);
label=lb11;}
Th2=buildmesh(newb11((XYb11.m-1))+b12(2)+b21(6),fixeborder=1);
plot(Th2,wait=1,cmm=“Th2”);
int reg1=Th2(0.1,0.1).region,
reg2=Th2(1.9,1.9).region;
mesh Th1=trunc(Th2,region==reg1);
//plot(Th1,wait=1,cmm=“Th1/reg1”);
mesh th1=Th1;
plot(Th1,Th2,wait=1,cmm=“Th1, Th2”);
plot(th1,Th2,wait=1,cmm=“th1, Th2”);
fespace V2(Th2,P1); V2 u2=x+4y;
fespace V1(Th1,P1); V1 u1=x+4y;
for (int i=0; i < Th2.nt; i++) {
cout << "i2/Th2 = " << i << "; "
<< Th2[i][0] << " " << Th2[i][1] << " " << Th2[i][2] << endl;
if (i<Th1.nt) {
cout << "i1/Th1 = " << i << "; "
<< Th1[i][0] << " " << Th1[i][1] << " " << Th1[i][2] << endl;
cout << "i1/th1 = " << i << "; "
<< th1[i][0] << " " << th1[i][1] << " " << th1[i][2] << endl;
for (int j=0; j<3; j++) {
if (Th2[i][j]!=Th1[i][j])
cout << "nodes are different in triangle i = " << i << endl;
if (u2[Th2[i][j]] != u1[Th1[i][j]])
cout << "u2!=u1 in triangle i = " << i << endl;
} } }