-- FreeFem++ v4.6 (Thu Apr 2 15:47:38 CEST 2020 - git v4.6) Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : // thermal_design.edp 2 : // 3 : // Discussion: 4 : // 5 : // A thermal design problem. The wall is a composite material 6 : // made of one cheap and one expensive material. 7 : // 8 : // Location: 9 : // 10 : // http://people.sc.fsu.edu/~jburkardt/freefem_src/thermal_design/thermal_design.edp 11 : // 12 : // Modified: 13 : // 14 : // 17 May 2015 15 : // 16 : // Author: 17 : // 18 : // Florian De Vuyst 19 : // 20 : // Reference: 21 : // 22 : // Florian De Vuyst, 23 : // Numerical modeling of transport problems using freefem++ software - 24 : // with examples in biology, CFD, traffic flow and energy transfer, 25 : // HAL id: cel-00842234 26 : // https://cel.archives-ouvertes.fr/cel-00842234 27 : // 28 : cout << "\n"; 29 : cout << "thermal_design\n"; 30 : cout << " FreeFem++ version\n"; 31 : 32 : real Lx = 4; 33 : real Ly = 5; 34 : real lr = 1; 35 : real lx = 0.3; 36 : real Tr = 20; 37 : real Text = -5.0; 38 : real hy = 0.1; 39 : real hx = 0.1; 40 : real mu = 0.2; 41 : real kappaa = 10.0; 42 : real kappac = 1.0; 43 : real kappae = 0.008; 44 : 45 : real V = lx * Ly; 46 : real Ve = mu * V; 47 : real xem = Lx + hx; 48 : real deltaxe = Ve / ( Ly - 2.0 * hy ); 49 : // 50 : // Define boundaries. 51 : // 52 : border c1 ( t = 0, Lx ) { x = t; y = 0.0; } 53 : border c1w ( t = Lx, Lx + lx ) { x = t; y = 0.0; } 54 : border cext ( t = 0, Ly ) { x = Lx + lx; y = t; } 55 : border c3w ( t = Lx+lx, Lx ) { x = t; y = Ly; } 56 : border c3 ( t = Lx, 0.0 ) { x = t; y = Ly; } 57 : border c4 ( t = Ly, Ly/2+lr/2 ) { x = 0; y = t; } 58 : border cr ( t = Ly/2+lr/2,Ly/2-lr/2 ) { x = 0; y = t; } 59 : border c5 ( t = Ly/2-lr/2,0.0 ) { x = 0; y = t; } 60 : border c0 ( t = 0, Ly ) { x = Lx; y = t; } 61 : 62 : border ce1 ( t = xem, xem + deltaxe ) { x = t; y = hy; } 63 : border ce2 ( t = hy, Ly - hy ) { x = xem+deltaxe; y = t; } 64 : border ce3 ( t = xem+deltaxe, xem ) { x = t; y = Ly-hy; } 65 : border ce4 ( t = Ly-hy, hy ) { x = xem; y = t; } 66 : // 67 : // Generate the mesh. 68 : // 69 : mesh Th = buildmesh ( 70 : c1(30) + c1w(10) + cext(100) 71 : + c3w(10) + c3(30) + c4(15) + cr(15) + c5(15) + c0(100) 72 : + ce1(15) + ce2(500) + ce3(15) + ce4(500) ); 73 : 74 : plot ( Th, wait = 1, ps = "thermal_design_mesh.ps" ); 75 : // 76 : // Generate the second mesh. 77 : // 78 : mesh Th2 = buildmesh ( 79 : c1(30) + c0(100) + c3(30) + c4(15) + cr(15) + c5(15) ); 80 : 81 : plot ( Th, wait = 1, ps = "thermal_design_mesh2.ps" ); 82 : // 83 : // Pick points that can identify the two regions. 84 : // 85 : real region1 = Th ( Lx/2, Ly/2).region; 86 : real region2 = Th ( xem+deltaxe/2,Ly/2).region; 87 : // 88 : // Define the finite element space. 89 : // 90 : fespace Vh ( Th, P1 ); 91 : Vh kappa; 92 : Vh th; 93 : Vh vh; 94 : // 95 : // Define a piecewise linear function KAPPA which has the values 96 : // 97 : // KAPPAA in region1, 98 : // KAPPAE in region2, 99 : // KAPPAC elsewhere. 100 : // 101 : kappa = kappac + ( kappaa - kappac ) * ( region == region1 ) 102 : + ( kappae - kappac ) * ( region == region2 ); 103 : 104 : plot ( kappa, nbiso = 60, fill = 1, value = 1, ps = "thermal_design_kappa.ps" ); 105 : // 106 : // Define the weak form of the thermal design problem. 107 : // 108 : problem thermaldesign ( th, vh ) = 109 : int2d ( Th ) ( kappa * dx ( th ) * dx ( vh ) 110 : + kappa * dy ( th ) * dy ( vh ) ) 111 : + on ( cext, th = Text ) 112 : + on ( cr, th = Tr ); 113 : // 114 : // Solve the problem. 115 : // 116 : thermaldesign; 117 : // 118 : // Define colors to be associated with contour plot levels. 119 : // 120 : // Since the name "colorhsv" occurs nowhere else, is this a 121 : // recognized keyword with a default value? 122 : // 123 : real[int] colorhsv = [ 124 : 4.0 / 6.0, 1.0, 0.5, 125 : 4.0 / 6.0, 1.0, 1.0, 126 : 5.0 / 6.0, 1.0, 1.0, 127 : 1.0, 1.0, 1.0, 128 : 1.0, 0.5, 1.0 ]; 129 : // 130 : // Define the levels of interest for a contour plot. 131 : // 132 : real[int] viso(26); 133 : for ( int i = 0; i < viso.n; i++ ) 134 : { 135 : viso[i] = i - 5; 136 : } 137 : // 138 : // Plot T. 139 : // 140 : plot ( th, viso = viso(0:viso.n-1), value = 1, fill = 1, 141 : ps = "thermal_design_tfield.ps" ); 142 : // 143 : // Compute performance indicators defined on the second mesh: 144 : // J1 = mean value of T; 145 : // J2 = minimum value of T; 146 : // J3 = standard deviation of T. 147 : // 148 : fespace Vh2 ( Th2, P1 ); 149 : Vh2 th2 = th; 150 : real J1; 151 : real J2; 152 : real J3; 153 : 154 : J1 = int2d ( Th2 ) ( th2 ) / Th2.area; 155 : J2 = th2[].min; 156 : J3 = sqrt ( int2d ( Th2 ) ( ( th2 - J1 )^2 ) / Th2.area ); 157 : 158 : cout << " Mean temperature = " << J1 << "\n"; 159 : cout << " Minimum temperature = " << J2 << "\n"; 160 : cout << " Standard deviation << " << J3 << "\n"; 161 : // 162 : // Terminate. 163 : // 164 : cout << "\n"; 165 : cout << "thermal_design:\n"; 166 : cout << " Normal end of execution.\n"; 167 : 168 : sizestack + 1024 =3272 ( 2248 ) thermal_design FreeFem++ version -- mesh: Nb of Triangles = 19207, Nb of Vertices 9717 -- mesh: Nb of Triangles = 6705, Nb of Vertices 3456 -- Solve : min -5 max 20 Mean temperature = 17.339 Minimum temperature = 14.1899 Standard deviation << 0.855644 thermal_design: Normal end of execution. times: compile 0.005435s, execution 1.50797s, mpirank:0 CodeAlloc : nb ptr 4086, size :495128 mpirank: 0 Ok: Normal End