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”);
