added Robin BC, disabled cuthill_mckee_reordering for now

mgrid_2/geom.cpp

@@ -759,6 +759,11 @@ const std::vector<int> Mesh::BoundaryEdges() const
     return _ebedges;
 }
 
+const std::vector<int> Mesh::BoundaryEdgeNodes() const
+{
+    return _bedges;
+}
+
 
 // Only the outer edges for Robin BC
 void Mesh::InitializeOuterEdges()
@@ -996,12 +1001,16 @@ Mesh::Mesh(std::string const &fname)
     //Debug();
     
     DeriveEdgeFromVertexBased();        // Generate also the edge based information
+
+
+
+
     
     {
 // Cuthill-McKee reordering
 //     Increases mesh generation time by factor 5 -  but solver is faster.
-            auto const perm = cuthill_mckee_reordering(_edges);
+            // auto const perm = cuthill_mckee_reordering(_edges);
-            PermuteVertices_EdgeBased(perm);
+            // PermuteVertices_EdgeBased(perm);
     }
     DeriveVertexFromEdgeBased();    
     

mgrid_2/geom.h

@@ -331,6 +331,7 @@ private:
     
 public:
     const std::vector<int> BoundaryEdges() const;
+    const std::vector<int> BoundaryEdgeNodes() const;
 
     std::vector<int> OuterEdgesSubdomains;
     std::vector<int> OuterEdges; 

mgrid_2/getmatrix.cpp

@@ -616,57 +616,48 @@ void FEM_Matrix::ApplyDirichletBC(std::vector<double> const &u, std::vector<doub
 }
 
 
-void FEM_Matrix::ApplyRobinBC_mult(std::vector<double> const &u, std::vector<double> &f, const double u_out)
+void FEM_Matrix::ApplyRobinBC_mult(std::vector<double> &f, const double u_out)
 {
     auto const RobinEdges = _mesh.OuterEdges;
     auto const RobinEdgesSubdomains = _mesh.OuterEdgesSubdomains;
 
-    for (int i = 0; i < RobinEdges.size(); ++i)
+    auto const BoundaryEdges = _mesh.BoundaryEdges();
+    auto const BoundaryEdgeNodes = _mesh.BoundaryEdgeNodes();
+
+    assert (BoundaryEdgeNodes.size() == 2* BoundaryEdges.size());
+
+    vector<double> Coordinates = _mesh.GetCoords();
+
+    for (size_t i = 0; i < RobinEdges.size(); ++i)
     {
-        cout << "Edge number " << RobinEdges[i] << ", subdomain: " << RobinEdgesSubdomains[i] << endl;
+        //cout << "Edge number " << RobinEdges[i] << ", subdomain: " << RobinEdgesSubdomains[i] << "        " << endl;
+        double alpha = Heat_transfer_coefficient(RobinEdgesSubdomains[i]);
 
-    }
+
+        int const EdgeNode1 = BoundaryEdgeNodes[2*i];
+        int const EdgeNode2 = BoundaryEdgeNodes[2*i + 1];
+
+        double x_1 = Coordinates[EdgeNode1];
+        double y_1 = Coordinates[EdgeNode1 + 1];
+        double x_2 = Coordinates[EdgeNode2];
+        double y_2 = Coordinates[EdgeNode2 + 1];
+        double EdgeLength = sqrt((x_2 - x_1)*(x_2 - x_1) +  (y_2 - y_1)*(y_2 - y_1));
 
 
 
+        int ii = _id[EdgeNode1];
+        int jj = _id[EdgeNode2];
+        int ij = fetch(_ik[_id[EdgeNode1]], EdgeNode1);
+        int ji = fetch(_ik[_id[EdgeNode2]], EdgeNode2);
 
+        _sk[ii] += EdgeLength*alpha/3;
+        _sk[jj] += EdgeLength*alpha/3;
+        _sk[ij] += EdgeLength*alpha/6;
+        _sk[ji] += EdgeLength*alpha/6;
 
-    
+        f[ii] += EdgeLength*alpha*u_out/2;
+        f[jj] += EdgeLength*alpha*u_out/2;
 
-
-    // Jakob Todo
-    auto const idx = _mesh.Index_DirichletNodes(); // ALL boundary nodes
-
-    const vector<int> sd_vec = _mesh.ElementSubdomains;
-
-
-    for (int i = 0; i < idx.size(); ++i) {
-        int const row = idx[i];
-        int subdomain = sd_vec[row];
-        // cout << "row: " << row;
-        // cout << ", subdomain: " << subdomain << endl;
-        double alpha = Heat_transfer_coefficient(subdomain);
-
-
-
-        f[row] += u_out*alpha/2;
-
-
-        for (int ij = _id[row]; ij < _id[row + 1]; ++ij)
-        {
-            int const col = _ik[ij];
-            if(col == row)
-            {
-                _sk[ij] += alpha/3;
-            }
-            else
-            {
-                int const id1 = fetch(col, row); // Find entry (col,row)
-                assert(id1 >= 0);
-                _sk[id1] += alpha/6;
-            }
-
-        }
     }
 
     return;
@@ -681,7 +672,7 @@ double FEM_Matrix::Heat_transfer_coefficient(const int subdomain)
     {
         // outside
         case 0:
-            alpha = 1.0;
+            alpha = 10.0;
             break;
 
         // ceramic
@@ -691,12 +682,12 @@ double FEM_Matrix::Heat_transfer_coefficient(const int subdomain)
 
         // water
         case 2:
-            alpha = 1.0;
+            alpha = 500.0;
             break;
 
         // air
         case 3:
-            alpha = 1.0;
+            alpha = 10.0;
             break;
 
         default:

mgrid_2/getmatrix.h

@@ -410,11 +410,10 @@ class FEM_Matrix: public CRS_Matrix
         * The <a href="https://www.jstor.org/stable/2005611?seq=1#metadata_info_tab_contents">penalty method</a>
         * is used for incorporating the given values @p u.
         *
-        * @param[in]     u (global) vector with Dirichlet data
         * @param[in]     u_out outside temperature
         * @param[in,out] f load vector
         */
-       void ApplyRobinBC_mult(std::vector<double> const &u, std::vector<double> &f, const double u_out);
+       void ApplyRobinBC_mult(std::vector<double> &f, const double u_out);
 
        double Heat_transfer_coefficient(const int subdomain);
 

mgrid_2/main.cpp

@@ -33,7 +33,7 @@ int main(int argc, char **argv )
     //gMesh_Hierarchy ggm(mesh_c, nrefine);
     //const Mesh &mesh = ggm.finest();
     //mesh.Debug();
-    //mesh.DebugEdgeBased();
+    //mesh_c.DebugEdgeBased();
 
 
     // Initializing FEM matrix !pattern! (only zero entries now)
@@ -55,7 +55,7 @@ int main(int argc, char **argv )
     // Calculate RHS
     SK.CalculateRHS(fv, [](double x, double y) {return 0;});
     //SK.CheckRowSum();
-    SK.CheckMatrix();
+    //SK.CheckMatrix();
 
 
     // Initialize temperature
@@ -65,7 +65,7 @@ int main(int argc, char **argv )
     mesh_c.Init_Solution_mult(uv, 2, [](double x, double y) -> double { return 18; });    // air
 
     // Apply BC
-    SK.ApplyRobinBC_mult(uv, fv, 18.0);
+    SK.ApplyRobinBC_mult(fv, 18.0);
 
 
     // Solve
@@ -74,7 +74,7 @@ int main(int argc, char **argv )
 
     auto t3 = system_clock::now(); // start timer
 
-    // JacobiSolve(SK, fv, uv );          // solve the system of equations
+    //JacobiSolve(SK, fv, uv );          // solve the system of equations
 
     auto t4 = system_clock::now();  // stop timer
     auto duration = duration_cast<microseconds>(t4 - t3);        // duration in microseconds