akiss
March 17, 2023, 3:59pm
1
Hello!
Given a mesh defined in FreeFem (in 2D for the moment), I would be interested by a distance map to a part of its boundary as a finite element field on it.
For instance, if I have a square mesh with label 1 for the bottom, label 2 for the right side, label 3 for the top and label 4 for the left side, I would like to have a distance map to the boundary part defined by labels 1 and 2.
What would be the way to have that?
Thanks a lot in advance!
You can use the plugin distante,
see :
load "distance"
mesh Th;
Th=square(100,100);
if(1)
{
mesh Th2;
mesh3 Th3;
int ltube=1,lup=3,lout=4,lin =2;
int l3d = 5;
real D3d=1;
{
int Nbtube= 4;
real Dx=1,Dy=1,La=3*Dx,Lh=5*Dy;
int Ndx=19,Ndy=19, Nh=(Lh-Dy)*Ndy/Dy,Na=(La-Dx)*Ndx/Dx ;
mesh Th0;
{
int[int] ll0=[1,0,0,1];
This file has been truncated. show original
load "distance"
mesh Th=buildmesh("g.gmesh",nbvx=10000);
fespace Vh(Th,P1);
real[int] bb(4);
boundingbox(Th,bb);
// b[0] xmin, xmax, ymin ymax
real xc1 = bb[0]*0.6 + bb[1]*0.4;
real xc2 = bb[0]*0.2 + bb[1]*0.8;
real yc = bb[2]*0.3 + bb[3]*0.7;
real r2 = ((bb[3]-bb[2])*.1)^2;
cout << xc1 << " "<< xc2 << " " << yc << "r2 = " << r2 << endl;
cout << bb << endl;
func f0 = ((x-xc1)^2 + (y-yc)^2- r2)*((x-xc2)^2 + (y-yc)^2 - r2);
Vh u = f0, v;
plot(u, dim=3, value=1,wait=1);
verbosity=1;
This file has been truncated. show original
to build the iso value of the border
load "distance"
int[int] lb=labels(Th);:
varf vb(u,v) = on(lb,u=1):
Vh ub= vb(0,Vh,tgv=1);// 1 on all arder and zero inside
ub[] = ub[] ? 0 : 1; // 0 on border and 1 inside
Vh dist;
distance(Th,ub,dist);
plot(ub,wait=1);