Zero solution by solving the Laplace problem on part of the mes

General Discussion
1

Hello everyone,
I wanted to solve the Laplace problem on part of the mesh domain. After defining the region that interests me, the problem does not resolve.
<<catch an erreur in solve => set sol = 0 !!! >>

load “msh3”
load “medit”
load “UMFPACK64”

int gam0 = 95;
int gamin = 96;
int gamout = 97;
int Sminus = 98;
int Splus = 99 ;

real [int] A(2), B(2), C(2), D(2), E(2), a(2), b(2), c(2), d(2);
real [int] F(2), G(2), H(2), I(2), J(2);
real [int] K(2), L(2), M(2), Ng(2), AA(2) , BB(2), CC(2), DD(2);
A = [6,1.5]; B = [6, 2.5]; C = [12, 1.5]; D = [21,0];
F = [15, 2]; I = [12, 2.5]; J = [21, 4]; M = [21, 1]; Ng = [21, 3];
AA = [6,-1]; BB=[6,5]; CC=[21,-1]; DD=[21,5];
a=[8,1.8]; b=[10,1.8]; c=[10,2.2]; d=[8,2.2];

border c1(t=0,1){x=(1-t)A[0]+tC[0]; y=(1-t)A[1]+tC[1];label=gam0;}
border c3(t=0,1){x=(1-t)C[0]+tD[0]; y=(1-t)C[1]+tD[1];label=gam0;}

border c8(t=0,1){x=(1-t)D[0]+tM[0]; y=(1-t)D[1]+tM[1];label=gamout;}

border c9(t=0,1){x=(1-t)M[0]+tF[0]; y=(1-t)M[1]+tF[1];label=gam0;}
border c10(t=0,1){x=(1-t)F[0]+tNg[0]; y=(1-t)F[1]+tNg[1];label=gam0;}

border c11(t=0,1){x=(1-t)Ng[0]+tJ[0]; y=(1-t)Ng[1]+tJ[1];label=gamout;}

border c12(t=0,1){x=(1-t)J[0]+tI[0]; y=(1-t)J[1]+tI[1];label=gam0;}
border c13(t=0,1){x=(1-t)I[0]+tB[0]; y=(1-t)I[1]+tB[1];label=gam0;}

border c14(t=0,1){x=(1-t)B[0]+tA[0]; y=(1-t)B[1]+tA[1];label=gamin;}

border c001(t=0,1){x=(1-t)a[0]+td[0]; y=(1-t)a[1]+td[1];gam0;}
border c002(t=0,1){x=(1-t)d[0]+tc[0]; y=(1-t)d[1]+tc[1];gam0;}
border c003(t=0,1){x=(1-t)c[0]+tb[0]; y=(1-t)c[1]+tb[1];gam0;}
border c004(t=0,1){x=(1-t)b[0]+ta[0]; y=(1-t)b[1]+ta[1];gam0;}

border c01(t=0,1){x=(1-t)AA[0]+tCC[0]; y=(1-t)AA[1]+tCC[1];Splus;}
border c01b(t=0,1){x=(1-t)CC[0]+tD[0]; y=(1-t)CC[1]+tD[1];Splus;}
border c02(t=0,1){x=(1-t)M[0]+tNg[0]; y=(1-t)M[1]+tNg[1];Splus;}
border c02b(t=0,1){x=(1-t)J[0]+tDD[0]; y=(1-t)J[1]+tDD[1];Splus;}
border c03(t=0,1){x=(1-t)DD[0]+tBB[0]; y=(1-t)DD[1]+tBB[1];Sminus;}
border c04(t=0,1){x=(1-t)BB[0]+tB[0]; y=(1-t)BB[1]+tB[1]; Splus;}
border c05(t=0,1){x=(1-t)A[0]+tAA[0]; y=(1-t)A[1]+tAA[1];Splus;}

plot(c001(1) +c01b(1)+c002(1)+c02b(1)+c003(1)+c004(1) +c1(1)+c3(1)+c8(1)+c9(9)+c10(1)+c11(1)
+ c12(1)+c13(1)+c14(1)+c01(1)+c02(1)+c03(1)+c04(1)+c05(1)) ;

mesh Th = buildmesh(c1(50)+c3(60)+ c8(10) +c9(50)+ c10(50)+c11(10)+c12(60)+c13(50)+c14(10)
+ c01(110)+c01b(10)+c02(15)+c02b(10)+c03(110)+c04(20)+c05(20)
+ c001(-7)+c002(-30)+c003(-7)+c004(-30)) ;

plot(Th, wait=true) ;

fespace Vh(Th, P2) ;

int road = Th(12., 2.).region ;
int air1 = Th(9., 2.).region ;
int air2 = Th(9., 1.).region ;
int air3 = Th(9., 4.).region ;
int air4 = Th(19., 2.).region ;

cout << "--------- Label number for road region = " << road << endl;
cout << "--------- Label number for air1 region = " << air1 << endl;
cout << "--------- Label number for air2 region = " << air2 << endl;
cout << "--------- Label number for air3 region = " << air3 << endl;
cout << "--------- Label number for air4 region = " << air4 << endl;

real U0 = 10. ;
real U00 = 0. ;

Vh U = U0*(region==road ? 1 : 0) + U00*(region==air1 ? 1 : 0) + U00*(region==air2 ? 1 : 0) + U00*(region==air3 ? 1 : 0) + U00*(region==air4 ? 1 : 0) ;

plot(U, fill=true, value = true);

Vh vv, vvh ;

problem Laplace(vv, vvh, solver=UMFPACK) = int2d(Th)(U*(dx(vv)*dx(vvh) + dy(vv)*dy(vvh)))
+ int1d(Th, gamin)(vvh)
+ on(gamout, vv=0)
;

Laplace ;
plot(vv, nbiso=50, fill=1, value=1);

Could you help me please?
Thanks

2

Are you sure that your problem is well-posed? Probably your matrix is not invertible because the factor U vanishes on some part of the domain (U00=0).
If you want to solve the Laplace problem on the road region you need to define a mesh corresponding to this region and set the unknown vv to have degrees of freedom only in this road region.

3

Hello everyone,
I wanted to solve the Laplace problem on part of the mesh domain. After defining the region that interests me, the problem does not resolve.
<<catch an erreur in solve => set sol = 0 !!! >>

load “msh3”
load “medit”
load “UMFPACK64”

int gam0 = 95;
int gamin = 96;
int gamout = 97;
int Sminus = 98;
int Splus = 99 ;

real [int] A(2), B(2), C(2), D(2), E(2), a(2), b(2), c(2), d(2);
real [int] F(2), G(2), H(2), I(2), J(2);
real [int] K(2), L(2), M(2), Ng(2), AA(2) , BB(2), CC(2), DD(2);
A = [6,1.5]; B = [6, 2.5]; C = [12, 1.5]; D = [21,0];
F = [15, 2]; I = [12, 2.5]; J = [21, 4]; M = [21, 1]; Ng = [21, 3];
AA = [6,-1]; BB=[6,5]; CC=[21,-1]; DD=[21,5];
a=[8,1.8]; b=[10,1.8]; c=[10,2.2]; d=[8,2.2];

border c1(t=0,1){x=(1-t)A[0]+tC[0]; y=(1-t)A[1]+tC[1];label=gam0;}
border c3(t=0,1){x=(1-t)C[0]+tD[0]; y=(1-t)C[1]+tD[1];label=gam0;}

border c8(t=0,1){x=(1-t)D[0]+tM[0]; y=(1-t)D[1]+tM[1];label=gamout;}

border c9(t=0,1){x=(1-t)M[0]+tF[0]; y=(1-t)M[1]+tF[1];label=gam0;}
border c10(t=0,1){x=(1-t)F[0]+tNg[0]; y=(1-t)F[1]+tNg[1];label=gam0;}

border c11(t=0,1){x=(1-t)Ng[0]+tJ[0]; y=(1-t)Ng[1]+tJ[1];label=gamout;}

border c12(t=0,1){x=(1-t)J[0]+tI[0]; y=(1-t)J[1]+tI[1];label=gam0;}
border c13(t=0,1){x=(1-t)I[0]+tB[0]; y=(1-t)I[1]+tB[1];label=gam0;}

border c14(t=0,1){x=(1-t)B[0]+tA[0]; y=(1-t)B[1]+tA[1];label=gamin;}

border c001(t=0,1){x=(1-t)a[0]+td[0]; y=(1-t)a[1]+td[1];gam0;}
border c002(t=0,1){x=(1-t)d[0]+tc[0]; y=(1-t)d[1]+tc[1];gam0;}
border c003(t=0,1){x=(1-t)c[0]+tb[0]; y=(1-t)c[1]+tb[1];gam0;}
border c004(t=0,1){x=(1-t)b[0]+ta[0]; y=(1-t)b[1]+ta[1];gam0;}

border c01(t=0,1){x=(1-t)AA[0]+tCC[0]; y=(1-t)AA[1]+tCC[1];Splus;}
border c01b(t=0,1){x=(1-t)CC[0]+tD[0]; y=(1-t)CC[1]+tD[1];Splus;}
border c02(t=0,1){x=(1-t)M[0]+tNg[0]; y=(1-t)M[1]+tNg[1];Splus;}
border c02b(t=0,1){x=(1-t)J[0]+tDD[0]; y=(1-t)J[1]+tDD[1];Splus;}
border c03(t=0,1){x=(1-t)DD[0]+tBB[0]; y=(1-t)DD[1]+tBB[1];Sminus;}
border c04(t=0,1){x=(1-t)BB[0]+tB[0]; y=(1-t)BB[1]+tB[1]; Splus;}
border c05(t=0,1){x=(1-t)A[0]+tAA[0]; y=(1-t)A[1]+tAA[1];Splus;}

plot(c001(1) +c01b(1)+c002(1)+c02b(1)+c003(1)+c004(1) +c1(1)+c3(1)+c8(1)+c9(9)+c10(1)+c11(1)

  • c12(1)+c13(1)+c14(1)+c01(1)+c02(1)+c03(1)+c04(1)+c05(1)) ;

mesh Th = buildmesh(c1(50)+c3(60)+ c8(10) +c9(50)+ c10(50)+c11(10)+c12(60)+c13(50)+c14(10)

  • c01(110)+c01b(10)+c02(15)+c02b(10)+c03(110)+c04(20)+c05(20)
  • c001(-7)+c002(-30)+c003(-7)+c004(-30)) ;

plot(Th, wait=true) ;

fespace Vh(Th, P2) ;

int road = Th(12., 2.).region ;
int air1 = Th(9., 2.).region ;
int air2 = Th(9., 1.).region ;
int air3 = Th(9., 4.).region ;
int air4 = Th(19., 2.).region ;

cout << "--------- Label number for road region = " << road << endl;
cout << "--------- Label number for air1 region = " << air1 << endl;
cout << "--------- Label number for air2 region = " << air2 << endl;
cout << "--------- Label number for air3 region = " << air3 << endl;
cout << "--------- Label number for air4 region = " << air4 << endl;

real U0 = 10. ;
real U00 = 0. ;

Vh U = U0*(region==road ? 1 : 0) + U00*(region==air1 ? 1 : 0) + U00*(region==air2 ? 1 : 0) + U00*(region==air3 ? 1 : 0) + U00*(region==air4 ? 1 : 0) ;

plot(U, fill=true, value = true);

Vh vv, vvh ;

problem Laplace(vv, vvh, solver=UMFPACK) = int2d(Th)(U*(dx(vv)*dx(vvh) + dy(vv)*dy(vvh)))

  • int1d(Th, gamin)(vvh)
  • on(gamout, vv=0)
    ;

Laplace ;
plot(vv, nbiso=50, fill=1, value=1);

Could you help me please?
Thanks

4

You could probably use the CG solver to get an answer or, IIRC I’ve played with this, use an “indcator”
fespace to enforce the laplace or identity ( u=1) in the various regions to make
it invertible.

5

The problem is well posed and can be easily solved (only) on the road region. However, my goal is to immerse this road region in a larger region to subsequently model the pollution. So I thought of defining two domains (road + air). On the first one, I wanted to start by solving the Laplace problem and then the problem of pollution on the whole domain (air + road).

6

can you upload the thing? my copy/paste messes it up.
thanks.

7

Transport-mod-bis (5).edp (3.5 KB)
Here is the code to download.
Using CG solver did not solve the problem.

8

Is this what you wanted?

Transport-mod.edp (3.6 KB)

diff Transpo*
80,84c80,81
< // problem Laplace(vv, vvh, solver=CG) = int2d(Th)(U*(dx(vv)dx(vvh) + dy(vv)dy(vvh)))
< problem Laplace(vv, vvh) = int2d(Th)(U(dx(vv)dx(vvh) + dy(vv)dy(vvh)))
< +int2d(Th)((1-U)((vv)(vvh) + (vv)(vvh)))
<
< + int1d(Th, gamin)(vvh)

problem Laplace(vv, vvh, solver=CG) = int2d(Th)(U*(dx(vv)*dx(vvh) + dy(vv)*dy(vvh)))
+ int1d(Th, gamin)(vvh)

9

Anis,
The answer depends on the boundary conditions you want. If you solve a Laplace problem for vv on the “road” region, you have to say what are the associated conditions you want on the boundary of the “road” domain.
Is it “vanishing normal derivative” everywhere, except on “gamout” where you set vv=0?
I am not sure that the proposition of Mike does that.
Or maybe “realistic” boundary conditions are more involved and depend on what happens on the “air” region?
Francois.

10

Not exactly.
I should have a solution that resembles the result of this code if I run the calculation on the road region only.
Transport-mod-bis.edp (2.0 KB)

11

yes, you need to fix the boundary conditions to get an exact match.
See “natural boundary conditions” to get them for the road-only example.
I took some of this and your mesh and ended up with this code
Transport-mod.edp (3.9 KB)

https://community.freefem.org/t/normal-boundary-condition/341

12

Since you have the right solution vv from your code on the road region only, what is your problem?
If you like you can extend this solution outside the road region by solving another Laplace problem for vvext, set in the air region only, with nonhomogeneous Dirichlet boundary condition vvext=vv at the interface between the two regions (and for example Neumann condition on the external boundary).

13

If I define and discretize only the road region, the Laplace problem can be solved on this region (no worries). However, if I define and discretize the entire domain (air + road) I can’t figure out how to calculate the solution of the problem only on the road region despite that I have defined the two regions correctly: Vh U = U0*(region==road? 1 : 0) + U00*(region==air1 ? 1 : 0) + U00*(region==air2 ? 1 : 0) + U00*(region==air3 ? 1 : 0) + U00*(region==air4 ? 1 : 0) ;
In other words, how could I integrate this definition into the weak formulation in order to have a solution on the road region only?
According to [MIKE MARCHYWKA], it is necessary to modify the conditions on gam0 ie at the interface between the two regions. Even taking vv=0 on gam0, it doesn’t work.
Thank you in advance for your help.
Anis

14

I think that gluing the two pieces I described above is equivalent to take your first version of the code and simply setting

real U00 = 1e-8 ;

in order to make it invertible.

15

Thanks for the suggestion. This allowed me to avoid a null solution, however I obtain get a solution on the whole domain and not only on the road.

16

You can make a cutoff by setting zero outside the road:

 fespace Vhr(Th, P2dc) ;
 Vhr vvr;
 vvr=vv*(region==road ? 1 : 0);
 plot(vvr, nbiso=50, fill=1, value=1, wait=true);

full file is attached.
laplace.edp (3.8 KB)

17

Zero is still a something. I think if you take the complete solution and also make your
road-only mesh you can just assign one solution to the other and ff will use
the points as needed.

In reality though the “road only” solution imposes a “natural” bndary condition that won’t
exist in the full problem unless there is something about it that makes it so.
Probably you have issues like sources and maybe diffusion rates that vary
without an actual boundary condition.

18

The boundary conditions used in the code posted by Anis as “Transport-mod-bis.edp” that solves the “road” alone is normal derivative vanishing (Neumann) on the interface between “road” and “air”.
The code for the full domain with U00 = 1e-8 provides the same solution with the same boundary condition in the “road” domain, which as I said is extended outside by nonhomogeneous Dirichlet for the “air” part.
If you consider another boundary condition in the “road” part, the code for the full domain will not get it (in the present version).

19

What was wrong with the code I posted earlier using tgv on rhs and the normal along the
road boundary? I guess if you need non-zero entries along A diagonal that may be easier
than adding new entries.

20

Thank you François for the code << laplace.edp >>. It allows you to calculate a solution only on the road after having calculated it over the entire domain (road + air).
Moreover and as you mentioned above, a Dirichlet condition on the boundary of the road “gam0=0” with a homogeneous Neumann condition on gamout change everything and give a solution that does not seem realistic to me. Could you tell me what to do to have the right code with these two conditions. It is clear that the problem comes from the condition “gam0=0”.
Here is the code with the new conditions :
laplace (1).edp (3.9 KB)
Thank you in advance