fixed subdomain indices, CHANGED thermalCondicutivity of ceramic

mgrid_2/geom.cpp

@@ -799,14 +799,14 @@ void Mesh::InitializeOuterEdges()
     assert(EdgeSubdomains.size() == 2*_ebedges.size());
     for (size_t k = 0; k < _ebedges.size(); ++k) // loop over all boundary edges
     {
-        if (EdgeSubdomains[2*k] == 0)
+        if (EdgeSubdomains[2*k] == -1)
         {
             OuterEdges.push_back(_ebedges[k]);
             OuterEdgesSubdomains.push_back(EdgeSubdomains[2*k + 1]);
             OuterEdgesNodes.push_back(_bedges[2*k]);
             OuterEdgesNodes.push_back(_bedges[2*k + 1]);
         }
-        if (EdgeSubdomains[2*k + 1] == 0)
+        if (EdgeSubdomains[2*k + 1] == -1)
         {            
             OuterEdges.push_back(_ebedges[k]);
             OuterEdgesSubdomains.push_back(EdgeSubdomains[2*k]);
@@ -1155,6 +1155,7 @@ void Mesh::ReadVertexBasedMesh(std::string const &fname)
         for (int k = 0; k < nbedges * 2; ++k)
         {
             ifs >> EdgeSubdomains[k];
+            EdgeSubdomains[k] -= OFFSET;
         }
     }
     else

mgrid_2/getmatrix.cpp

@@ -428,12 +428,13 @@ void FEM_Matrix::CalculateLaplace_mult(vector<double> &f)
 
 double FEM_Matrix::ThermalConductivity(const int subdomain)
 {
+    // source: https://www.engineeringtoolbox.com/thermal-conductivity-d_429.html
     double lambda = 0.0;
     switch (subdomain)
     {
         // ceramic mug
         case 0:
-            lambda = 4.0;           // anything from 1 to 4
+            lambda = 6.0;           // anything from 1 to 4
             break;
 
         // water
@@ -447,7 +448,7 @@ double FEM_Matrix::ThermalConductivity(const int subdomain)
             break;
 
         default:
-            lambda = 1.0;
+            assert(false && "BAD subdomain");
             break;
     }
 
@@ -456,6 +457,7 @@ double FEM_Matrix::ThermalConductivity(const int subdomain)
 
 double FEM_Matrix::VolumetricHeatCapacity(const int subdomain)
 {
+    // source? https://en.wikipedia.org/wiki/Table_of_specific_heat_capacities
     double c = 0.0;
     switch (subdomain)
     {
@@ -465,8 +467,8 @@ double FEM_Matrix::VolumetricHeatCapacity(const int subdomain)
             break;
 
         // water
-        case 1:
+        case 1:    
-            c = 4.184 * 1e6;       
+            c = 4.184 * 1e6;
             break;
 
         // air
@@ -475,7 +477,7 @@ double FEM_Matrix::VolumetricHeatCapacity(const int subdomain)
             break;
 
         default:
-            c = 1.0;
+            assert(false && "BAD subdomain");
             break;
     }
 
@@ -630,13 +632,13 @@ void FEM_Matrix::ApplyRobinBC_mult(std::vector<double> &f, const double u_out)
         int const subdomain = RobinEdgesSubdomains[i];
         double const alpha = Heat_transfer_coefficient(subdomain);
 
-        //int const EdgeNumber = RobinEdges[i];
+        int const EdgeNumber = RobinEdges[i];
 
         int const EdgeNode1 = RobinEdgeNodes[2*i];
         int const EdgeNode2 = RobinEdgeNodes[2*i + 1];
         
-        //cout << "Edge number " << EdgeNumber << ", subdomain: " << subdomain << "        ";
+        cout << "Edge number " << EdgeNumber << ", subdomain: " << subdomain << "        ";
-        //cout << "Node1: " << EdgeNode1 << "     Node2: " << EdgeNode2 << "            " << endl;
+        cout << "Node1: " << EdgeNode1 << "     Node2: " << EdgeNode2 << "            " << endl;
 
         double x_1 = Coordinates[2*EdgeNode1];
         double y_1 = Coordinates[2*EdgeNode1 + 1];
@@ -666,35 +668,25 @@ void FEM_Matrix::ApplyRobinBC_mult(std::vector<double> &f, const double u_out)
     return;
 }
 
+// source: https://www.engineeringtoolbox.com/overall-heat-transfer-coefficient-d_434.html
 double FEM_Matrix::Heat_transfer_coefficient(const int subdomain)
 {
-    int matlab_sd_index = subdomain - 1;
     double alpha = 0.0;
 
-    switch (matlab_sd_index)
+    switch (subdomain)
     {
-        // outside
+        // ceramic to air
         case 0:
-            alpha = 10.0;
-            break;
-
-        // ceramic
-        case 1:
             alpha = 1.0;
             break;
 
-        // water
+        // air to air (Free Convection Gas - Free Convection Gas : U = 1 - 2 W/m2K  (typical window, room to outside air through glass))
         case 2:
-            alpha = 500.0;
+            alpha = 1.0;
-            break;
-
-        // air
-        case 3:
-            alpha = 10.0;
             break;
 
         default:
-            alpha = 1.0;
+            assert(false && "BAD subdomain");
             break;
     }
 

mgrid_2/main.cpp

@@ -31,7 +31,7 @@ int main(int argc, char **argv )
     // ##########################################
 
     double dt = 1.0;           // time step
-    //int steps = 200;            // number of time iterations
+    int steps = 100;            // number of time iterations
 
     double u0_mug = 18.0;
     double u0_fluid = 80.0;
@@ -80,7 +80,9 @@ int main(int argc, char **argv )
 
     auto t3 = system_clock::now(); // start timer
     double average_cup_temperature = u0_mug;
-    while (average_cup_temperature < 67)
+
+    double time_count = 0;
+    while (average_cup_temperature < 67.0)
     //for (int step = 0; step < steps; ++step)
     {
         vector<double> G(Mdt.Nrows(), 0.0);        
@@ -96,7 +98,8 @@ int main(int argc, char **argv )
                                          //        ----- SK -----               ------ H -------
 
         average_cup_temperature = mesh_c.AverageVectorFunction_perSubdomain(uv, 0);
-        cout << "Average cup temperature: " << average_cup_temperature << endl;
+        cout << "Average cup temperature: " << average_cup_temperature << " after " << time_count << " seconds. " << endl;
+        time_count += dt;
     }
     auto t4 = system_clock::now();  // stop timer