Effectively 1D flux boundary condition in 3D simulation

General Discussion
1

I want to model the following phenomena:
two bodies with thickness d made of the same material A with the same conductivity k are separated with a very thin (delta << d) plate from material B with a very high conductivity (k2 >> k). To reduce the computational burden of meshing the plate, it is treated as a 2-dimensional surface and is introduced to the weak form as int2d(plate) k2 delta nabla_t u nabla_t v, where nabla_t is 2D gradient in the plate’s plane (not to confuse with a directional derivative notation).

Everything works well, but I have a problem with Robin’s/Neumann’s boundary condition. If I introduce the BC in a standard way. (int2d(gamma_N) v*q), the flux will be imposed only on the surface of material A, and effectively there will be a no-flux condition at the boundary of material B (*). The flux at the B’s boundary can have a crucial effect on the phenomena, so I would like to control it.

Is it possible from the FEM and FreeFEM perspective to introduce the flux on such a boundary?

From the mathematical perspective on a rectangular domain, int2d f(x) dxdy = deltay int1d f(x)dx, so theoretically I could introduce delta int1d(plate) v q, but int1d is not accepted in 3D formulation by FreeFEM, so I am looking for other solutions.

The following MRE script is what I am working on currently:

load "gmsh"
load "msh3"
mesh3 Th=gmshload3("cube.msh");
fespace Vh(Th, P2);
Vh u, v;

real k = 1;
real kratio = 1000;
real k2 = kratio*k;
real delta = 0.01;
real q=100;
macro grad(u) [dx(u), dy(u), dz(u)]//EOM
solve Laplace(u, v)
= int3d(Th)(
     k*( grad(u)'*grad(v) )//' 
     )
+ int2d(Th,2)(  //thin plate effect
	k2*delta*(dx(u)*dx(v)+dy(u)*dy(v))  
	)
- int2d(Th,3)( //Neumann BC at inlet
	v*q
	)
+ on(5, u=0) //Dirichlet BC at outlet
;
plot(u, wait=true, nbiso=10, value=true);

//Flux calculations
meshL ThL = Sline(100);
real xmin = 0;
real xmax = 1;
real ymin = 0;
real ymax = ymin;
real zmin = 0.5;
real zmax = zmin;
meshL ThL1=movemesh(ThL, [x*(xmax-xmin)+xmin,x*(ymax-ymin)+ymin,x*(zmax-zmin)+zmin]);
meshL ThL2=movemesh(ThL, [x*(xmax-xmin)+xmin,x*(ymax-ymin)+ymin+1,x*(zmax-zmin)+zmin]);
fespace VhL1(ThL1, P2);
fespace VhL2(ThL2, P2);
VhL1 ul1;
VhL2 ul2;

real fluxIn = -int2d(Th,3)(k*dy(u)) ;//'
real fluxOut = -int2d(Th,5)(k*([dx(u), dy(u), dz(u)])'*[N.x,N.y, N.z]) ;//'
real fluxInL = -int1d(ThL1)(k2*delta*dy(u));
real fluxOutL = -int1d(ThL2)(k2*delta*dy(u));
cout<<"\nFlux in: "<<fluxIn<<", "<<fluxInL<<", sum: "<<fluxIn+fluxInL<<"\nFlux out: "<<fluxOut<<", "<<fluxOutL<<", sum: "<<fluxOut+fluxOutL<<"\n\n";

The mesh file:
https://mega.nz/file/SvZRgTBa#_RnDooT0YIIdpq8bIfncfaKIcN83vZue1TLCaHbQpDw
Labels: 3 - inlet, 5 - outlet, 100 - line at inlet, 101 - line at outlet, 2 - plate surface

(*) In the following script, there is a flux associated with the boundary, but it is correctly approaching 0 when increasing the mesh density

2

Hello,
You need to write a linear form on your (3d) fespace, representing the integral on an edge (or family of edges).
You can try to define directly the associated vector. For that you have to build the list of dof that are located on this edge, and then for each of these dof define the associated coefficient.

In the case the fespace is P1 the dofs are the vertices of the mesh and this can be done as

int[int] listdofedge(Th.nv);// list of dofs on edge, Th.nv is an upper bound on the number
int count=0;
for (int iv=0;iv<Th.nv;iv++){// loop on vertices
 real xx=Th(iv).x,yy=Th(iv).y,zz=Th(iv).z;// coordinates of vertex
 if ((abs(zz-0.5)<1.e-7)&&((abs(xx*(1.-xx))<1.e-7)||(abs(yy*(1.-yy))<1.e-7))) {// if vertex is on edge
  listdofedge(count)=iv;// select this vertex in list
  count++;
 }
}
listdofedge.resize(count);
cout << "Th.nv " << Th.nv << endl;
cout << "found " << count << " dof on edge"<< endl;

real[int] linedge(Th.nv);//linear form of the integral on edge
linedge=0.;
for (int ic=0;ic<count;ic++){//loop on dof on edge
 linedge(listdofedge(ic))=1./10.;// delta_x=1/10
}

Then you have to replace solve by varf, define your matrix and right-hand side,
and add delta*q*linedge to your right-hand side before solving the linear system.

3

I appreciate your help. I was able to implement your solution for P1 and it seems like it works:
3D_bd_flux_BC.edp (2.6 KB)
I don’t know why I had to multiply the flux by 100 (edge length is 1, so the specific flux should have the same value as the total flux), but that’s a minor concern.

I have a lot of questions, but let me ask the crucial ones:

  • we need to manually obtain vertices by checking their location. Is it possible to filter them by region? The mesh has lines labeled 100 and 101, they can be isolated in other software like Paraview. Also, how do we get “hidden” nodes when using P2?
  • we need to manually perform the integration. Here it is easy, as for constant q and known element length, we just need to integrate the associated shape function. Is there a way we can generalize this to non-uniform mesh, P2, and functional q? I think if we could get the line elements with associated nodes and their type (“hidden” and “normal” nodes), it would be possible.
4

I don’t think it is possible to select the vertices by edge region. There is a label attached to each vertex, Th(iv).label but it refers to the interior or 2d boundary, not to 1d boundary.
Indeed reading your mesh file, FreeFem says “edges in 3D mesh are not considered yet in freefem++, skeep data”.
About your plan to define the linear form in general, it is a nightmare. It is probably better to wait for future updates of FF that could manage with line elements in 3d!

About getting all the nodes for P2, this is possible.
Look at the documentation p.237 “finite elements connectivity”

For P2 in 2d, setting fespace Vh(Th,P2); the quantity Vh(k,j) is the index as dof of Vh, where k is the index of a triangle and j=0,…,5 the dof in triangle k. (j=0,1,2 is for the vertices, j=3,4,5 is for the middle of edges, 3 is the edge opposite to vertex 0, 4 is for the edge opposite to vertex 1, 5 is for the edge opposite to vertex 2)
In particular, setting

Vh u;
u[]=0;
u[](Vh(k,j))=1;

You get that u is the basis function of the dof of index j in triangle k.