Simplify boundary conditions

zhukunw · February 1, 2025, 6:37am

Hi all, I am solving a convection diffusion equation on a domain made up of many short line segments. I am wondering if there’s a shortcut for setting the boundaries and conditions so that I do not have to manually do that each time. Please let me know if you have seen any way before. Thanks so much! Here’s the code: // Adjustable parameters
int meshsizer = 6010; // Mesh size, how many grid/little chambers are on one side
int meshsizezfull = 2010; // Mesh size, how many grid/little chambers are on one side
int meshsizez = 2*10;
// THIS SHOULD BE READ FROM A SPECIFIED DOCUMENT FROM PYTHON
real k = 0.134; // [mm^2/s] Thermal Diffusivity (at average temperature of water domain)
real Tleft = 273.15; // [K] Left boundary temperature
real dTdr = 0.025; // [K/mm] 0.025
real deltaT = 30.0; // [K] after meeting on 1/31/2025, 8. changing deltaT doesn’t change the ratio
real sigmax = 45.0; // [mm] after meeting on 1/31/2025, decreasing sigmax does increase the ratio closer to 1.

//Load boundary points’ x and y coordinates
ifstream bdryy(“C:/Users/zhuku/Desktop/Figure_34hr/boundary_y_coordinates.txt”); // function mask_ice_out
real[int] ycoors(75);
ifstream bdryx(“C:/Users/zhuku/Desktop/Figure_34hr/boundary_x_coordinates.txt”);
real[int] xcoors(75);
string boundary;
for (int i = 74; i >= 0; i–){
string boundary;
getline(bdryy, boundary);
ycoors(i) = atof(boundary);
getline(bdryx, boundary);
xcoors(i) = atof(boundary);
}
real imagelength = 150.5; //[mm]
real imageheight = 50.; //[mm]
real icesurfacesidelen = xcoors[0]; // [mm]
real icebottomsidelen = xcoors[xcoors.n-1]; // [mm]
real surfacesidelen = imagelength - icesurfacesidelen; // [mm] imagelength - xcoors[0]
real bottomsidelen = imagelength - icebottomsidelen; // [mm] imagelength - xcoors[36]
real waterheight = ycoors[0];// [mm] Don’t need to specify the index either
//func Tsurface = Tleft + dTdr * (x - icesurfacesidelen); //[K] water surface boundary temperature
func Tsurface = Tleft + deltaT * erf(x/sigmax); //[K] water surface boundary temperature

mesh Mesh;
border bottomside(t=0., 1.){
x=icebottomsidelen + bottomsidelen*t;
y=0.;
label=1;};

border rightside(t=0., 1.){
x=icebottomsidelen + bottomsidelen;
y=waterheight * t;
label=2;};

border watersurface(t=0., 1.){
x=icesurfacesidelen + surfacesidelen * (1.-t);
cout << x;
y=waterheight;
label=3;};

// 1/24/2025 Define a batch to replace the indexes with i, think of array which collects all the boundaries
// Try automize the border and see if it produces the right boundary
border leftside1(t=0., 1.) {x = (1 - t) * xcoors[0] + t * xcoors[1]; y = (1 - t) * ycoors[0] + t * ycoors[1]; label = 4;};
border leftside2(t=0., 1.) {x = (1 - t) * xcoors[1] + t * xcoors[2]; y = (1 - t) * ycoors[1] + t * ycoors[2]; label = 5;};
border leftside3(t=0., 1.) {x = (1 - t) * xcoors[2] + t * xcoors[3]; y = (1 - t) * ycoors[2] + t * ycoors[3]; label = 6;};
border leftside4(t=0., 1.) {x = (1 - t) * xcoors[3] + t * xcoors[4]; y = (1 - t) * ycoors[3] + t * ycoors[4]; label = 7;};
border leftside5(t=0., 1.) {x = (1 - t) * xcoors[4] + t * xcoors[5]; y = (1 - t) * ycoors[4] + t * ycoors[5]; label = 8;};
border leftside6(t=0., 1.) {x = (1 - t) * xcoors[5] + t * xcoors[6]; y = (1 - t) * ycoors[5] + t * ycoors[6]; label = 9;};
border leftside7(t=0., 1.) {x = (1 - t) * xcoors[6] + t * xcoors[7]; y = (1 - t) * ycoors[6] + t * ycoors[7]; label = 10;};
border leftside8(t=0., 1.) {x = (1 - t) * xcoors[7] + t * xcoors[8]; y = (1 - t) * ycoors[7] + t * ycoors[8]; label = 11;};
border leftside9(t=0., 1.) {x = (1 - t) * xcoors[8] + t * xcoors[9]; y = (1 - t) * ycoors[8] + t * ycoors[9]; label = 12;};

border leftside10(t=0., 1.) {x = (1 - t) * xcoors[9] + t * xcoors[10]; y = (1 - t) * ycoors[9] + t * ycoors[10]; label = 13;};
border leftside11(t=0., 1.) {x = (1 - t) * xcoors[10] + t * xcoors[11]; y = (1 - t) * ycoors[10] + t * ycoors[11]; label = 14;};
border leftside12(t=0., 1.) {x = (1 - t) * xcoors[11] + t * xcoors[12]; y = (1 - t) * ycoors[11] + t * ycoors[12]; label = 15;};
border leftside13(t=0., 1.) {x = (1 - t) * xcoors[12] + t * xcoors[13]; y = (1 - t) * ycoors[12] + t * ycoors[13]; label = 16;};
border leftside14(t=0., 1.) {x = (1 - t) * xcoors[13] + t * xcoors[14]; y = (1 - t) * ycoors[13] + t * ycoors[14]; label = 17;};
border leftside15(t=0., 1.) {x = (1 - t) * xcoors[14] + t * xcoors[15]; y = (1 - t) * ycoors[14] + t * ycoors[15]; label = 18;};
border leftside16(t=0., 1.) {x = (1 - t) * xcoors[15] + t * xcoors[16]; y = (1 - t) * ycoors[15] + t * ycoors[16]; label = 19;};
border leftside17(t=0., 1.) {x = (1 - t) * xcoors[16] + t * xcoors[17]; y = (1 - t) * ycoors[16] + t * ycoors[17]; label = 20;};
border leftside18(t=0., 1.) {x = (1 - t) * xcoors[17] + t * xcoors[18]; y = (1 - t) * ycoors[17] + t * ycoors[18]; label = 21;};

border leftside19(t=0., 1.) {x = (1 - t) * xcoors[18] + t * xcoors[19]; y = (1 - t) * ycoors[18] + t * ycoors[19]; label = 22;};
border leftside20(t=0., 1.) {x = (1 - t) * xcoors[19] + t * xcoors[20]; y = (1 - t) * ycoors[19] + t * ycoors[20]; label = 23;};
border leftside21(t=0., 1.) {x = (1 - t) * xcoors[20] + t * xcoors[21]; y = (1 - t) * ycoors[20] + t * ycoors[21]; label = 24;};
border leftside22(t=0., 1.) {x = (1 - t) * xcoors[21] + t * xcoors[22]; y = (1 - t) * ycoors[21] + t * ycoors[22]; label = 25;};
border leftside23(t=0., 1.) {x = (1 - t) * xcoors[22] + t * xcoors[23]; y = (1 - t) * ycoors[22] + t * ycoors[23]; label = 26;};
border leftside24(t=0., 1.) {x = (1 - t) * xcoors[23] + t * xcoors[24]; y = (1 - t) * ycoors[23] + t * ycoors[24]; label = 27;};
border leftside25(t=0., 1.) {x = (1 - t) * xcoors[24] + t * xcoors[25]; y = (1 - t) * ycoors[24] + t * ycoors[25]; label = 28;};
border leftside26(t=0., 1.) {x = (1 - t) * xcoors[25] + t * xcoors[26]; y = (1 - t) * ycoors[25] + t * ycoors[26]; label = 29;};
border leftside27(t=0., 1.) {x = (1 - t) * xcoors[26] + t * xcoors[27]; y = (1 - t) * ycoors[26] + t * ycoors[27]; label = 30;};

border leftside28(t=0., 1.) {x = (1 - t) * xcoors[27] + t * xcoors[28]; y = (1 - t) * ycoors[27] + t * ycoors[28]; label = 31;};
border leftside29(t=0., 1.) {x = (1 - t) * xcoors[28] + t * xcoors[29]; y = (1 - t) * ycoors[28] + t * ycoors[29]; label = 32;};
border leftside30(t=0., 1.) {x = (1 - t) * xcoors[29] + t * xcoors[30]; y = (1 - t) * ycoors[29] + t * ycoors[30]; label = 33;};
border leftside31(t=0., 1.) {x = (1 - t) * xcoors[30] + t * xcoors[31]; y = (1 - t) * ycoors[30] + t * ycoors[31]; label = 34;};
border leftside32(t=0., 1.) {x = (1 - t) * xcoors[31] + t * xcoors[32]; y = (1 - t) * ycoors[31] + t * ycoors[32]; label = 35;};
border leftside33(t=0., 1.) {x = (1 - t) * xcoors[32] + t * xcoors[33]; y = (1 - t) * ycoors[32] + t * ycoors[33]; label = 36;};
border leftside34(t=0., 1.) {x = (1 - t) * xcoors[33] + t * xcoors[34]; y = (1 - t) * ycoors[33] + t * ycoors[34]; label = 37;};
border leftside35(t=0., 1.) {x = (1 - t) * xcoors[34] + t * xcoors[35]; y = (1 - t) * ycoors[34] + t * ycoors[35]; label = 38;};
border leftside36(t=0., 1.) {x = (1 - t) * xcoors[35] + t * xcoors[36]; y = (1 - t) * ycoors[35] + t * ycoors[36]; label = 39;};

border leftside37(t=0., 1.) {x = (1 - t) * xcoors[36] + t * xcoors[37]; y = (1 - t) * ycoors[36] + t * ycoors[37]; label = 40;};
border leftside38(t=0., 1.) {x = (1 - t) * xcoors[37] + t * xcoors[38]; y = (1 - t) * ycoors[37] + t * ycoors[38]; label = 41;};
border leftside39(t=0., 1.) {x = (1 - t) * xcoors[38] + t * xcoors[39]; y = (1 - t) * ycoors[38] + t * ycoors[39]; label = 42;};
border leftside40(t=0., 1.) {x = (1 - t) * xcoors[39] + t * xcoors[40]; y = (1 - t) * ycoors[39] + t * ycoors[40]; label = 43;};
border leftside41(t=0., 1.) {x = (1 - t) * xcoors[40] + t * xcoors[41]; y = (1 - t) * ycoors[40] + t * ycoors[41]; label = 44;};
border leftside42(t=0., 1.) {x = (1 - t) * xcoors[41] + t * xcoors[42]; y = (1 - t) * ycoors[41] + t * ycoors[42]; label = 45;};
border leftside43(t=0., 1.) {x = (1 - t) * xcoors[42] + t * xcoors[43]; y = (1 - t) * ycoors[42] + t * ycoors[43]; label = 46;};
border leftside44(t=0., 1.) {x = (1 - t) * xcoors[43] + t * xcoors[44]; y = (1 - t) * ycoors[43] + t * ycoors[44]; label = 47;};
border leftside45(t=0., 1.) {x = (1 - t) * xcoors[44] + t * xcoors[45]; y = (1 - t) * ycoors[44] + t * ycoors[45]; label = 48;};
border leftside46(t=0., 1.) {x = (1 - t) * xcoors[45] + t * xcoors[46]; y = (1 - t) * ycoors[45] + t * ycoors[46]; label = 49;};
border leftside47(t=0., 1.) {x = (1 - t) * xcoors[46] + t * xcoors[47]; y = (1 - t) * ycoors[46] + t * ycoors[47]; label = 50;};
border leftside48(t=0., 1.) {x = (1 - t) * xcoors[47] + t * xcoors[48]; y = (1 - t) * ycoors[47] + t * ycoors[48]; label = 51;};
border leftside49(t=0., 1.) {x = (1 - t) * xcoors[48] + t * xcoors[49]; y = (1 - t) * ycoors[48] + t * ycoors[49]; label = 52;};
border leftside50(t=0., 1.) {x = (1 - t) * xcoors[49] + t * xcoors[50]; y = (1 - t) * ycoors[49] + t * ycoors[50]; label = 53;};
border leftside51(t=0., 1.) {x = (1 - t) * xcoors[50] + t * xcoors[51]; y = (1 - t) * ycoors[50] + t * ycoors[51]; label = 54;};
border leftside52(t=0., 1.) {x = (1 - t) * xcoors[51] + t * xcoors[52]; y = (1 - t) * ycoors[51] + t * ycoors[52]; label = 55;};
border leftside53(t=0., 1.) {x = (1 - t) * xcoors[52] + t * xcoors[53]; y = (1 - t) * ycoors[52] + t * ycoors[53]; label = 56;};
border leftside54(t=0., 1.) {x = (1 - t) * xcoors[53] + t * xcoors[54]; y = (1 - t) * ycoors[53] + t * ycoors[54]; label = 57;};
border leftside55(t=0., 1.) {x = (1 - t) * xcoors[54] + t * xcoors[55]; y = (1 - t) * ycoors[54] + t * ycoors[55]; label = 58;};
border leftside56(t=0., 1.) {x = (1 - t) * xcoors[55] + t * xcoors[56]; y = (1 - t) * ycoors[55] + t * ycoors[56]; label = 59;};
border leftside57(t=0., 1.) {x = (1 - t) * xcoors[56] + t * xcoors[57]; y = (1 - t) * ycoors[56] + t * ycoors[57]; label = 60;};
border leftside58(t=0., 1.) {x = (1 - t) * xcoors[57] + t * xcoors[58]; y = (1 - t) * ycoors[57] + t * ycoors[58]; label = 61;};
border leftside59(t=0., 1.) {x = (1 - t) * xcoors[58] + t * xcoors[59]; y = (1 - t) * ycoors[58] + t * ycoors[59]; label = 62;};
border leftside60(t=0., 1.) {x = (1 - t) * xcoors[59] + t * xcoors[60]; y = (1 - t) * ycoors[59] + t * ycoors[60]; label = 63;};
border leftside61(t=0., 1.) {x = (1 - t) * xcoors[60] + t * xcoors[61]; y = (1 - t) * ycoors[60] + t * ycoors[61]; label = 64;};
border leftside62(t=0., 1.) {x = (1 - t) * xcoors[61] + t * xcoors[62]; y = (1 - t) * ycoors[61] + t * ycoors[62]; label = 65;};
border leftside63(t=0., 1.) {x = (1 - t) * xcoors[62] + t * xcoors[63]; y = (1 - t) * ycoors[62] + t * ycoors[63]; label = 66;};
border leftside64(t=0., 1.) {x = (1 - t) * xcoors[63] + t * xcoors[64]; y = (1 - t) * ycoors[63] + t * ycoors[64]; label = 67;};
border leftside65(t=0., 1.) {x = (1 - t) * xcoors[64] + t * xcoors[65]; y = (1 - t) * ycoors[64] + t * ycoors[65]; label = 68;};
border leftside66(t=0., 1.) {x = (1 - t) * xcoors[65] + t * xcoors[66]; y = (1 - t) * ycoors[65] + t * ycoors[66]; label = 69;};
border leftside67(t=0., 1.) {x = (1 - t) * xcoors[66] + t * xcoors[67]; y = (1 - t) * ycoors[66] + t * ycoors[67]; label = 70;};
border leftside68(t=0., 1.) {x = (1 - t) * xcoors[67] + t * xcoors[68]; y = (1 - t) * ycoors[67] + t * ycoors[68]; label = 71;};
border leftside69(t=0., 1.) {x = (1 - t) * xcoors[68] + t * xcoors[69]; y = (1 - t) * ycoors[68] + t * ycoors[69]; label = 72;};
border leftside70(t=0., 1.) {x = (1 - t) * xcoors[69] + t * xcoors[70]; y = (1 - t) * ycoors[69] + t * ycoors[70]; label = 73;};
border leftside71(t=0., 1.) {x = (1 - t) * xcoors[70] + t * xcoors[71]; y = (1 - t) * ycoors[70] + t * ycoors[71]; label = 74;};
border leftside72(t=0., 1.) {x = (1 - t) * xcoors[71] + t * xcoors[72]; y = (1 - t) * ycoors[71] + t * ycoors[72]; label = 75;};
border leftside73(t=0., 1.) {x = (1 - t) * xcoors[72] + t * xcoors[73]; y = (1 - t) * ycoors[72] + t * ycoors[73]; label = 76;};
border leftside74(t=0., 1.) {x = (1 - t) * xcoors[73] + t * xcoors[74]; y = (1 - t) * ycoors[73] + t * ycoors[74]; label = 77;};

// Build the mesh the function is called the buildmesh
Mesh = buildmesh(
bottomside(meshsizer)
+ rightside(meshsizezfull)
+ watersurface(meshsizer)
+ leftside1(meshsizez)
+ leftside2(meshsizez)
+ leftside3(meshsizez)
+ leftside4(meshsizez)
+ leftside5(meshsizez)
+ leftside6(meshsizez)
+ leftside7(meshsizez)
+ leftside8(meshsizez)
+ leftside9(meshsizez)
+ leftside10(meshsizez)
+ leftside11(meshsizez)
+ leftside12(meshsizez)
+ leftside13(meshsizez)
+ leftside14(meshsizez)
+ leftside15(meshsizez)
+ leftside16(meshsizez)
+ leftside17(meshsizez)
+ leftside18(meshsizez)
+ leftside19(meshsizez)
+ leftside20(meshsizez)
+ leftside21(meshsizez)
+ leftside22(meshsizez)
+ leftside23(meshsizez)
+ leftside24(meshsizez)
+ leftside25(meshsizez)
+ leftside26(meshsizez)
+ leftside27(meshsizez)
+ leftside28(meshsizez)
+ leftside29(meshsizez)
+ leftside30(meshsizez)
+ leftside31(meshsizez)
+ leftside32(meshsizez)
+ leftside33(meshsizez)
+ leftside34(meshsizez)
+ leftside35(meshsizez)
+ leftside36(meshsizez)
+ leftside37(meshsizez)
+ leftside38(meshsizez)
+ leftside39(meshsizez)
+ leftside40(meshsizez)
+ leftside41(meshsizez)
+ leftside42(meshsizez)
+ leftside43(meshsizez)
+ leftside44(meshsizez)
+ leftside45(meshsizez)
+ leftside46(meshsizez)
+ leftside47(meshsizez)
+ leftside48(meshsizez)
+ leftside49(meshsizez)
+ leftside50(meshsizez)
+ leftside51(meshsizez)
+ leftside52(meshsizez)
+ leftside53(meshsizez)
+ leftside54(meshsizez)
+ leftside55(meshsizez)
+ leftside56(meshsizez)
+ leftside57(meshsizez)
+ leftside58(meshsizez)
+ leftside59(meshsizez)
+ leftside60(meshsizez)
+ leftside61(meshsizez)
+ leftside62(meshsizez)
+ leftside63(meshsizez)
+ leftside64(meshsizez)
+ leftside65(meshsizez)
+ leftside66(meshsizez)
+ leftside67(meshsizez)
+ leftside68(meshsizez)
+ leftside69(meshsizez)
+ leftside70(meshsizez)
+ leftside71(meshsizez)
+ leftside72(meshsizez)
+ leftside73(meshsizez)
+ leftside74(meshsizez)
);

//Load uvel and wvel values from a file
ifstream fuvel(“C:/Users/zhuku/Desktop/water_side_temp_dist/U_velocity.txt”); // this file only contains one column of data which is the U velocity
real[int] Uvel(Mesh.nv);
ifstream fwvel(“C:/Users/zhuku/Desktop/water_side_temp_dist/W_velocity.txt”);
real[int] Wvel(Mesh.nv);

string line;
for (int i =0; i < Mesh.nv; i++){
getline(fuvel, line);
Uvel(i) = atof(line);
getline(fwvel, line);
Wvel(i) = atof(line);
}

fespace Vh(Mesh, P1);
Vh T, v, uvel, wvel; //uvel and wvel also need to be a P1 field

plot(Mesh,wait=1);
// Think of saving the mesh first //Save mesh
savemesh(Mesh,“C:/Users/zhuku/Desktop/water_side_temp_dist/waterhshape.msh”);

fb77 · February 1, 2025, 2:06pm

Hello,
You have to define a multi-border as in Input Border Data from an External File in FreeFEM++ - #3 by fb77

The principle is to define
border a(t=0.,1.;i){ x=(1.-t)*bdrydatx(i)+t*bdrydatx(i+1); y=(1.-t)*bdrydaty(i)+t*bdrydaty(i+1);label=1;}

In this definition, i stands for an integer varying from 0 to n-1, where n is to be given when you plot or build the mesh by calling a().

Then you write
mesh Th=buildmesh(a(discr));
where discr is an array containing the number of segments that you want when you split each border element.
Then the above n is the length of discr.

Example (here n=3):

real[int] bdrydatx=[1,3,2,1];
real[int] bdrydaty=[0,1,4,0];
int[int] discr=[1,2,2];
border a(t=0.,1.;i){ x=(1.-t)*bdrydatx(i)+t*bdrydatx(i+1); y=(1.-t)*bdrydaty(i)+t*bdrydaty(i+1);label=1;}
mesh Th=buildmesh(a(discr));
plot(Th);

It gives

The lower border is split in 1 piece, the two other borders are split in 2 pieces because discr=[1,2,2].
You can use several multi-borders of this type in a command like
mesh Th=buildmesh(a1(discr1)+a2(discr2));

zhukunw · May 30, 2025, 6:28pm

Thank you so much! This has significantly simplified the code while producing the same result!

zhukunw · November 16, 2025, 8:42pm

Hello: I think I am having more questions regarding to the boundary that consists of multiple segments. The thing is I have build the mesh and solved the heat equation inside, but now I need to calculate heat flux along the boundary segment by segment. I really need to know the local heat flux. How can I break this down? I need to store the flux at each segment into an array. Here’s my code: //Solve the temperature distribution in ice THIS STAYS THE SAME AS NO PARAMETER TUNING IS INVOLVED

// Part 2 for ice

fespace Vhice(Thi, P1);

Vhice T, v; //uvel and wvel also need to be a P1 field

// Define the variables

real templeftice = 267.0; //266.55; // [K] Left boundary temperature all u change to T

real temprightice = 273.15;

real R = 20.5;

real kice = 1.15; //0.134; 0.58

solve thermic(T, v)

= int2d(Thi)(

    + (x+R) \* (kice\*(

          dx(T) \* dx(v)

        + dy(T) \* dy(v)

        )

    ))

// Neumann Conditions

+ int1d(Thi, bottomsideice)(0\*v)

+ int1d(Thi, icesurfaceice)(0\*v)

// Dirchlet Conditions

+ on(icefrontice, T=temprightice)

+ on(leftsideice, T=temprightice);

//real totalfluxright = 0.;

real totalflux = 0.;

int Nseg = 224; // or size of your arrays - 1

real[int] fluxsegmentice(Nseg); // store flux per segment

for (int i = 0; i < Nseg; i++) {

fluxsegmentice(i) = int1d(Thi, icefrontice(i))(-kice \* (dx(T)\*N.x + dy(T)\*N.y) \* (x + R)) \* 2 \* pi;

}

I think there is something wrong in doing it this way.

zhukunw · November 16, 2025, 8:44pm

Here’s how I defined the mesh by the way: // Define the borders and the meshes

border bottomsideice(t=0.,1.){ x=(1.-t)*xcoorsverticesice(0)+t*xcoorsverticesice(1); y=(1.-t)*ycoorsverticesice(0)+t*ycoorsverticesice(1);label=1;}

int[int] numelebottomsideice = [20];

border icefrontice(t=0.,1.;i){ x=(1.-t)*xcoorsicefrontverticesice(i)+t*xcoorsicefrontverticesice(i+1); y=(1.-t)*ycoorsicefrontverticesice(i)+t*ycoorsicefrontverticesice(i+1);label=2;}

border icesurfaceice(t=0.,1.){ x=(1.-t)*xcoorsverticesice(225)+t*xcoorsverticesice(226); y=(1.-t)*ycoorsverticesice(225)+t*ycoorsverticesice(226);label=3;}

int[int] numeleicesurfaceice = [10];

border leftsideice(t=0.,1.){ x=(1.-t)*xcoorsverticesice(226)+t*xcoorsverticesice(0); y=(1.-t)*ycoorsverticesice(226)+t*ycoorsverticesice(0);label=4;}

int[int] numeleleftsideice = [30];

mesh Thi=buildmesh(bottomsideice(numelebottomsideice) + icefrontice(numelementseachedgeice) + icesurfaceice(numeleicesurfaceice) + leftsideice(numeleleftsideice));

fb77 · November 19, 2025, 10:50am

Dear Zhukun Wang,
can you upload the code as a file, and the necessary auxiliary files, so that I can run your code? Otherwise I cannot understand well what you are doing…

zhukunw · November 30, 2025, 3:23am

Icefront_as_oneboundary.edp (4.6 KB)

These are files describing the coordinates but I am not able to upload. The ice has a shape shown in picture.

My question is essentially just on “border icefront(t=0.,1.;i){ x=(1.-t)*xcoorsicefrontvertices(i)+t*xcoorsicefrontvertices(i+1); y=(1.-t)*ycoorsicefrontvertices(i)+t*ycoorsicefrontvertices(i+1);label=2;}“ Here it is a short cut that I am able to combine all the small segments into one boundary. Now, I need to calculate the heat flux from each of the segments after solving the heat equation(I am interested in calculating local heat fluxes but I do not want go back to previous approach where I set each segment as a boundary and repeat hundreds of times), is there a way that I can do this without adding too much addition code? Thank you so much.

zhukunw · November 30, 2025, 3:31am

Icefront_as_oneboundary.edp (4.6 KB)

These are files describing the coordinates but I am not able to upload. The ice has a shape shown in picture.

My question is essentially just on “border icefront(t=0.,1.;i){ x=(1.-t)*xcoorsicefrontvertices(i)+t*xcoorsicefrontvertices(i+1); y=(1.-t)*ycoorsicefrontvertices(i)+t*ycoorsicefrontvertices(i+1);label=2;}“ Here it is a short cut that I am able to combine all the small segments into one boundary. Now, I need to calculate the heat flux from each of the segments after solving the heat equation(I am interested in calculating local heat fluxes but I do not want go back to previous approach where I set each segment as a boundary and repeat hundreds of times), is there a way that I can do this without adding too much addition code? Thank you so much. Please let me know if you need the original data so that I can find a way to send it to you.

fb77 · December 6, 2025, 1:49pm

You can identify the piece of boundary by a label set to “i” the multiborder index, as for example

int nbp=20;//number of points
real[int] bdrydatx(nbp),bdrydaty(nbp);//coordinates of points on the boundary
for (int j=0;j<nbp;j++){
 bdrydatx(j)=cos(j*2.*pi/nbp);//example of a circle
 bdrydaty(j)=sin(j*2.*pi/nbp);
}
int[int] discr(nbp);//number of split for each segment
discr=1;// 1 means no split
border a(t=0.,1.;i){ x=(1.-t)*bdrydatx(i)+t*bdrydatx((i+1)%nbp); y=(1.-t)*bdrydaty(i)+t*bdrydaty((i+1)%nbp);label=i;}
// label has been set to i the multiborder index
mesh Th=buildmesh(a(discr));
plot(Th,wait=1);

fespace Vh(Th,P1);
Vh u,v;

solve laplace(u,v)=int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
                -int2d(Th)(v)//right-hand side
                +int1d(Th)(1.e10*u*v)//Dirichlet BC u=1 by explicit penalty
                -int1d(Th)(1.e10*v)
                ;
plot(u);

real[int] flux(nbp);//flux on each segment
for (int j=0;j<nbp;j++){
 flux(j)=int1d(Th,j)(dx(u)*N.x+dy(u)*N.y);// in int1d(Th,j), j indicates the label of the piece of boundary
}
cout << "flux " << flux << endl;

Here I have set Dirichlet BD by explicit penalty, because I don’t want to list all the labels.

zhukunw · December 7, 2025, 9:16pm

Thanks a lot. I think this works if I design my entire mesh by only 1 boundary. Now, my newer cases are consists of several boundaries, let’s say in my example, I need boundaries with names: bottomside, rightside, watersurface, offsetside, and icefront. These five together formed the entire mesh. Now, since I have define icefront in a slightly different way: “border icefront(t=0.,1.;i){ x=(1.-t)*xcoorsicefrontvertices(i)+t*xcoorsicefrontvertices(i+1); y=(1.-t)*ycoorsicefrontvertices(i)+t*ycoorsicefrontvertices(i+1);label=5;}“ where this boundary should consists of ~200 pieces, how can I only calculate heat flux in this specific boundary? I did try “totalfluxleft = int1d(Th, icefront)(-k * (dx(T) * N.x + dy(T) * N.y) * (x)) * 2 * pi;“ but this give the integrated along all these ~200 pieces at the icefront boundary, but I need the flux piece by piece.

fb77 · December 9, 2025, 11:02am

You need to set a different label for each piece you want to consider alone for computing the flux.
Then you write int1d(Th,labelnumber)()

zhukunw · January 2, 2026, 8:31pm

Clear. A follow up question, I have noticed that there are something strange happening at the solution boundaries in the FreeFEM. For example, if I zoom into the solution picture attached, why do I see the rainbow part at the bottom right edge? If you look closer, there is a part of the boundary that is colored by multiple colors, I am wondering why this happens. Here’s my code that I set for solving it: //Solve the temperature distribution in ice THIS STAYS THE SAME AS NO PARAMETER TUNING IS INVOLVED

// Part 2 for ice

fespace Vhice(Thi, P1);

Vhice T, v; //uvel and wvel also need to be a P1 field

// Define the variables

real templeftice = 267; //266.55; // [K] Left boundary temperature all u change to T

real temprightice = 273.15;

real kice = 1.02; //0.134; 0.58 1.15

real rhoice = 917.0;

real cpice = 2100.0;

solve thermic(T, v)

= int2d(Thi)(

    + (x-2.2) \* (kice\*(

          dx(T) \* dx(v)

        + dy(T) \* dy(v)

        )

    ))

// Neumann Conditions

+ int1d(Thi, bottomsideice)(0\*v)

+ int1d(Thi, icesurfaceice)(0\*v)

// Dirchlet Conditions

//+ on(icefrontice, T=temprightice)

+ on(iceright1, T=temprightice) + on(iceright2, T=temprightice) + on(iceright3, T=temprightice) + on(iceright4, T=temprightice) + on(iceright5, T=temprightice) + on(iceright6, T=temprightice)

+ on(iceright7, T=temprightice) + on(iceright8, T=temprightice) + on(iceright9, T=temprightice) + on(iceright10, T=temprightice) + on(iceright11, T=temprightice) + on(iceright12, T=temprightice)

+ on(iceright13, T=temprightice) + on(iceright14, T=temprightice) + on(iceright15, T=temprightice) + on(iceright16, T=temprightice) + on(iceright17, T=temprightice) + on(iceright18, T=temprightice)

+ on(iceright19, T=temprightice) + on(iceright20, T=temprightice) + on(iceright21, T=temprightice) + on(iceright22, T=temprightice) + on(iceright23, T=temprightice) + on(iceright24, T=temprightice)

+ on(iceright25, T=temprightice) + on(iceright26, T=temprightice) + on(iceright27, T=temprightice) + on(iceright28, T=temprightice) + on(iceright29, T=temprightice) + on(iceright30, T=temprightice)

+ on(iceright31, T=temprightice) + on(iceright32, T=temprightice) + on(iceright33, T=temprightice) + on(iceright34, T=temprightice) + on(iceright35, T=temprightice) + on(iceright36, T=temprightice)

+ on(iceright37, T=temprightice) + on(iceright38, T=temprightice) + on(iceright39, T=temprightice) + on(iceright40, T=temprightice) + on(iceright41, T=temprightice) + on(iceright42, T=temprightice)

+ on(iceright43, T=temprightice) + on(iceright44, T=temprightice) + on(iceright45, T=temprightice) + on(iceright46, T=temprightice) + on(iceright47, T=temprightice) + on(iceright48, T=temprightice)

+ on(iceright49, T=temprightice) + on(iceright50, T=temprightice) + on(iceright51, T=temprightice) + on(iceright52, T=temprightice) + on(iceright53, T=temprightice) + on(iceright54, T=temprightice)

+ on(iceright55, T=temprightice) + on(iceright56, T=temprightice) + on(iceright57, T=temprightice) + on(iceright58, T=temprightice) + on(iceright59, T=temprightice) + on(iceright60, T=temprightice)

+ on(iceright61, T=temprightice) + on(iceright62, T=temprightice) + on(iceright63, T=temprightice) + on(iceright64, T=temprightice) + on(iceright65, T=temprightice) + on(iceright66, T=temprightice)

+ on(iceright67, T=temprightice) + on(iceright68, T=temprightice) + on(iceright69, T=temprightice) + on(iceright70, T=temprightice) + on(iceright71, T=temprightice) + on(iceright72, T=temprightice)

+ on(iceright73, T=temprightice) + on(iceright74, T=temprightice) + on(iceright75, T=temprightice) + on(iceright76, T=temprightice) + on(iceright77, T=temprightice) + on(iceright78, T=temprightice)

+ on(iceright79, T=temprightice) + on(iceright80, T=temprightice) + on(iceright81, T=temprightice) + on(iceright82, T=temprightice) + on(iceright83, T=temprightice) + on(iceright84, T=temprightice)

+ on(iceright85, T=temprightice) + on(iceright86, T=temprightice) + on(iceright87, T=temprightice) + on(iceright88, T=temprightice) + on(iceright89, T=temprightice) + on(iceright90, T=temprightice)

+ on(iceright91, T=temprightice) + on(iceright92, T=temprightice) + on(iceright93, T=temprightice) + on(iceright94, T=temprightice) + on(iceright95, T=temprightice) + on(iceright96, T=temprightice)

+ on(iceright97, T=temprightice) + on(iceright98, T=temprightice) + on(iceright99, T=temprightice) + on(iceright100, T=temprightice) + on(iceright101, T=temprightice) + on(iceright102, T=temprightice)

+ on(iceright103, T=temprightice) + on(iceright104, T=temprightice) + on(iceright105, T=temprightice) + on(iceright106, T=temprightice) + on(iceright107, T=temprightice) + on(iceright108, T=temprightice)

+ on(iceright109, T=temprightice) + on(iceright110, T=temprightice) + on(iceright111, T=temprightice) + on(iceright112, T=temprightice) + on(iceright113, T=temprightice) + on(iceright114, T=temprightice)

+ on(iceright115, T=temprightice) + on(iceright116, T=temprightice) + on(iceright117, T=temprightice) + on(iceright118, T=temprightice) + on(iceright119, T=temprightice) + on(iceright120, T=temprightice)

+ on(iceright121, T=temprightice) + on(iceright122, T=temprightice) + on(iceright123, T=temprightice) + on(iceright124, T=temprightice) + on(iceright125, T=temprightice) + on(iceright126, T=temprightice)

+ on(iceright127, T=temprightice) + on(iceright128, T=temprightice) + on(iceright129, T=temprightice) + on(iceright130, T=temprightice) + on(iceright131, T=temprightice) + on(iceright132, T=temprightice)

+ on(iceright133, T=temprightice) + on(iceright134, T=temprightice) + on(iceright135, T=temprightice) + on(iceright136, T=temprightice) + on(iceright137, T=temprightice) + on(iceright138, T=temprightice)

+ on(iceright139, T=temprightice) + on(iceright140, T=temprightice) + on(iceright141, T=temprightice) + on(iceright142, T=temprightice) + on(iceright143, T=temprightice) + on(iceright144, T=temprightice)

+ on(iceright145, T=temprightice) + on(iceright146, T=temprightice) + on(iceright147, T=temprightice) + on(iceright148, T=temprightice) + on(iceright149, T=temprightice) + on(iceright150, T=temprightice)

+ on(iceright151, T=temprightice) + on(iceright152, T=temprightice) + on(iceright153, T=temprightice) + on(iceright154, T=temprightice) + on(iceright155, T=temprightice) + on(iceright156, T=temprightice)

+ on(iceright157, T=temprightice) + on(iceright158, T=temprightice) + on(iceright159, T=temprightice) + on(iceright160, T=temprightice) + on(iceright161, T=temprightice) + on(iceright162, T=temprightice)

+ on(iceright163, T=temprightice) + on(iceright164, T=temprightice) + on(iceright165, T=temprightice) + on(iceright166, T=temprightice) + on(iceright167, T=temprightice) + on(iceright168, T=temprightice)

+ on(iceright169, T=temprightice) + on(iceright170, T=temprightice) + on(iceright171, T=temprightice) + on(iceright172, T=temprightice) + on(iceright173, T=temprightice) + on(iceright174, T=temprightice)

+ on(iceright175, T=temprightice) + on(iceright176, T=temprightice) + on(iceright177, T=temprightice) + on(iceright178, T=temprightice) + on(iceright179, T=temprightice) + on(iceright180, T=temprightice)

+ on(iceright181, T=temprightice) + on(iceright182, T=temprightice) + on(iceright183, T=temprightice) + on(iceright184, T=temprightice) + on(iceright185, T=temprightice) + on(iceright186, T=temprightice)

+ on(iceright187, T=temprightice) + on(iceright188, T=temprightice) + on(iceright189, T=temprightice) + on(iceright190, T=temprightice) + on(iceright191, T=temprightice) + on(iceright192, T=temprightice)

+ on(iceright193, T=temprightice) + on(iceright194, T=temprightice) + on(iceright195, T=temprightice) + on(iceright196, T=temprightice) + on(iceright197, T=temprightice) + on(iceright198, T=temprightice)

+ on(iceright199, T=temprightice) + on(iceright200, T=temprightice) + on(iceright201, T=temprightice) + on(iceright202, T=temprightice) + on(iceright203, T=temprightice) + on(iceright204, T=temprightice)

+ on(iceright205, T=temprightice) + on(iceright206, T=temprightice) + on(iceright207, T=temprightice) + on(iceright208, T=temprightice) + on(iceright209, T=temprightice) + on(iceright210, T=temprightice)

+ on(iceright211, T=temprightice) + on(iceright212, T=temprightice) + on(iceright213, T=temprightice) + on(iceright214, T=temprightice) + on(iceright215, T=temprightice) + on(iceright216, T=temprightice)

+ on(iceright217, T=temprightice) + on(iceright218, T=temprightice) + on(iceright219, T=temprightice) + on(iceright220, T=temprightice) + on(iceright221, T=temprightice) + on(iceright222, T=temprightice)

+ on(iceright223, T=temprightice) + on(iceright224, T=temprightice)



+ on(leftsideice, T=templeftice);

fb77 · January 3, 2026, 3:01pm

I can’t understand what is happening because I do not have a runnable code.
Tip: instead of writing
varf myT(T,v)=on(3,T=temp)+on(17,T=temp)+on(21,T=temp);
you can do

int[int] mylabs;
mylabs=[3,17,21];
varf myT(T,v)=on(mylabs,T=temp);

moreover you can write a loop in order to define mylabs.

zhukunw · January 7, 2026, 6:23pm

Thanks! Also one more question on the Neumann Boundary conditions in GENERAL.

So if I have a no flux boundary condition for example, if I am solving a heat equation (Steady State) and I want to set a boundary that is of no flux condition, should I do “// Neumann Conditions

+ int1d(Thi, copperbasetoicebase)(0\*v)“ or “int1d(Thi, copperbasetoicebase)(kice \* 0 \* v \* (x))“ or “int1d(Thi, copperbasetoicebase) (kice \* 0 \* v \* dT(x))“ or none of this?

Note that “Thi” is referring to the mesh, “copperbasetoicebase“ is referring to a vertical side that is of no flux boundary condition. The problem is: solve thermic(T, v)

= int2d(Thi)(

    + (x) \* (kice\*(

          dx(T) \* dx(v)

        + dy(T) \* dy(v)

        )

    )) + ALL CONDITIONS I HAVE

fb77 · January 7, 2026, 6:52pm

You propose several ways to write zero. They are all equivalent: the result vanishes, that’s the only important thing.
The most simple is not to write any boundary integral.
This gives the variational formulation for Neumann BC.

Topic		Replies	Views
Calculating the heat flux at two boundries General Discussion	1	47	April 17, 2025
Boundary condition on half edge General Discussion	17	897	September 16, 2020
Boundary condition General Discussion	19	934	April 25, 2024
Heat Flux Boundary Condition for 1D Line Element General Discussion	3	170	November 12, 2024
How to implement an operator at the boundary? General Discussion	6	357	January 14, 2024

Simplify boundary conditions

Related topics