convert mesh to meters, more coefficients

generate_mesh/coffee_cup.m

@@ -2,12 +2,12 @@
 clear all
 clc
 
-diam_bottom = 50.0;
+diam_bottom = 0.05;
-diam_top = 83.0;
+diam_top = 0.083;
-wall_thickness = 3.0;
+wall_thickness = 0.003;
-bottom_thickness = 6.0;
+bottom_thickness = 0.006;
-top_level = 105.0;
+top_level = 0.105;
-fluid_level = 66.0;
+fluid_level = 0.066;
 floor_level = 0.0;
 
 inner_diam = @(height) diam_bottom - 2*wall_thickness + (diam_top - diam_bottom)/(top_level - bottom_thickness)*(height - bottom_thickness);
@@ -34,7 +34,7 @@ g=[2 -diam_bottom/2 diam_bottom/2  floor_level floor_level  1  0;        % #vert
    ]';
 
 
-[p,e,t] = initmesh(g,'hmax',1.5); 
+[p,e,t] = initmesh(g,'hmax',0.0015); 
 pdemesh(p,e,t)
 
 

generate_mesh/visualize_results.m

@@ -15,6 +15,14 @@ fname = 'uv.txt';
 
 [xc,ia,v] = ascii_read_meshvector(fname);
 
-h = trisurf(ia, xc(:,1), xc(:,2), v);
+%h = trisurf(ia, xc(:,1), xc(:,2), v);
+
+figure;
+h = patch('Faces', ia, 'Vertices', xc, 'FaceVertexCData', v, ...
+      'FaceColor', 'interp', 'EdgeColor', 'k');
+axis equal tight; 
+colorbar;
+colormap(jet);
+title('Heat distribution');
 
 waitfor(h)                     % wait for closing the figure

mgrid_2/getmatrix.cpp

@@ -410,7 +410,7 @@ void FEM_Matrix::CalculateLaplaceMult(vector<double> &f)
     for (int i = 0; i < nelem; ++i) {
 
         auto subdomain = sd_vec[i];
-        double lambda = Thermal_coefficient(subdomain);
+        double lambda = ThermalConductivity(subdomain);
 
         cout << subdomain << endl;
 
@@ -429,30 +429,23 @@ void FEM_Matrix::CalculateLaplaceMult(vector<double> &f)
     return;
 }
 
-double FEM_Matrix::Thermal_coefficient(const int subdomain)
+double FEM_Matrix::ThermalConductivity(const int subdomain)
 {
-    int matlab_sd_index = subdomain - 1;
     double lambda = 0.0;
-
+    switch (subdomain)
-    switch (matlab_sd_index)
     {
-        // outside
+        // ceramic mug
         case 0:
-            lambda = 0.026;
-            break;
-
-        // ceramic
-        case 1:
             lambda = 3.0;           // anything from 1 to 4
             break;
 
         // water
-        case 2:
+        case 1:
             lambda = 0.6;
             break;
 
         // air
-        case 3:
+        case 2:
             lambda = 0.026;         // depends on temperature actually
             break;
 
@@ -464,6 +457,33 @@ double FEM_Matrix::Thermal_coefficient(const int subdomain)
     return lambda;
 }
 
+double FEM_Matrix::VolumetricHeatCapacity(const int subdomain)
+{
+    double lambda = 0.0;
+    switch (subdomain)
+    {
+        // ceramic mug
+        case 1:
+            lambda = 2.0 * 1e6;
+            break;
+
+        // water
+        case 2:
+            lambda = 4.184 * 1e6;       
+            break;
+
+        // air
+        case 3:
+            lambda = 1.2 * 1e3;
+            break;
+
+        default:
+            lambda = 1.0;
+            break;
+    }
+
+    return lambda;
+}
 
 void FEM_Matrix::CalculateLaplace(vector<double> &f)
 {

mgrid_2/getmatrix.h

@@ -350,8 +350,12 @@ class FEM_Matrix: public CRS_Matrix
         */
        void CalculateLaplaceMult(std::vector<double> &f);
 
-       double Thermal_coefficient(const int subdomain);
+       // Returns thermal conductivity coefficient of subdomain [ kg * m / (s^3 * K) ]
+       double ThermalConductivity(const int subdomain);
 
+       // Returns volumetric heap capacity (specific heat capacity * density) coefficient of subdomain
+      //  c * rho [ J / (m^3 * K) ]
+       double VolumetricHeatCapacity(const int subdomain);
 
         /**
         * Calculates the entries of f.e. stiffness matrix for the Laplace operator

mgrid_2/visualize_results.m

@@ -18,7 +18,7 @@ fname = 'uv.txt';
 %h = trisurf(ia, xc(:,1), xc(:,2), v);
 
 figure;
-patch('Faces', ia, 'Vertices', xc, 'FaceVertexCData', v, ...
+h = patch('Faces', ia, 'Vertices', xc, 'FaceVertexCData', v, ...
       'FaceColor', 'interp', 'EdgeColor', 'k');
 axis equal tight; 
 colorbar;