added subdomain indexing for outer edges, for use in robin BC

mgrid_2/geom.cpp

@@ -761,24 +761,25 @@ const std::vector<int> Mesh::BoundaryEdges() const
 
 
 // Only the outer edges for Robin BC
-const std::vector<int> Mesh::OuterEdges() const
+void Mesh::InitializeOuterEdges()
 {
-    vector<int> outerEdges;
+    assert(EdgeSubdomains.size() == 2*_ebedges.size());
-
+    for (size_t k = 0; k < _ebedges.size(); ++k) // loop over all boundary edges
-    for (int k = 0; k < _ebedges.size()/2; ++k)
     {
-        if (EdgeSubdomains[k] == 0 || EdgeSubdomains[k + 1] == 0)
+        if (EdgeSubdomains[2*k] == 0)
-            outerEdges.push_back(_ebedges[k]);
+        {
+            OuterEdges.push_back(_ebedges[k]);
+            OuterEdgesSubdomains.push_back(EdgeSubdomains[2*k + 1]);
+        }
+        if (EdgeSubdomains[2*k + 1] == 0)
+        {
+            OuterEdges.push_back(_ebedges[k]);
+            OuterEdgesSubdomains.push_back(EdgeSubdomains[2*k]);
+        }
     }
 
     cout << "All boundary edges: " << _ebedges.size() << endl;
-    cout << "Outer boundary edges: " << outerEdges.size() << endl;
+    cout << "Outer boundary edges: " << OuterEdges.size() << endl;
-
-    // cout << _ebedges;
-    // cout << endl;
-    // cout << outerEdges << endl;
-
-    return outerEdges;
 }
 
 
@@ -1017,6 +1018,7 @@ Mesh::Mesh(std::string const &fname)
 Mesh::Mesh(std::string const &filename, std::string const &subdomain_filename) : Mesh(filename)
 {
     ElementSubdomains = ReadElementSubdomains(subdomain_filename);
+    InitializeOuterEdges();
 }
 
 

mgrid_2/geom.h

@@ -327,9 +327,14 @@ private:
     void Write_ascii_paraview_2D(std::string const &fname, std::vector<double> const &v) const;
     void Write_ascii_paraview_3D(std::string const &fname, std::vector<double> const &v) const;
 
+    void InitializeOuterEdges();
+    
 public:
     const std::vector<int> BoundaryEdges() const;
-    const std::vector<int> OuterEdges() const;
+
+    std::vector<int> OuterEdgesSubdomains;
+    std::vector<int> OuterEdges; 
+
 
 
      /**
@@ -565,8 +570,6 @@ public:
     */
     [[nodiscard]] const std::vector<int> ReadElementSubdomains(std::string const &dname) const;
 
-    [[nodiscard]] const std::vector<int> ReadEdgeSubdomains(std::string const &filename) const;
-
 
     /**
      * Calculates the largest inner angle in element @p idx.

mgrid_2/getmatrix.cpp

@@ -412,7 +412,7 @@ void FEM_Matrix::CalculateLaplaceMult(vector<double> &f)
         auto subdomain = sd_vec[i];
         double lambda = ThermalConductivity(subdomain);
 
-        cout << subdomain << endl;
+        //cout << subdomain << endl;
 
         CalcElemSpecific(ia.data() + 3 * i, xc.data(), lambda, ske);
         //AddElem(ia.data()+3 * i, ske, fe, _id.data(), _ik.data(), _sk.data(), f.data()); // GH: deprecated
@@ -423,9 +423,6 @@ void FEM_Matrix::CalculateLaplaceMult(vector<double> &f)
     cout << "finished in  " <<  duration  << " sec.    ########\n";  // ToDo: change to  systemclock
     //Debug();
 
-
-
-
     return;
 }
 
@@ -583,6 +580,99 @@ 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)
+{
+    auto const RobinEdges = _mesh.OuterEdges;
+    auto const RobinEdgesSubdomains = _mesh.OuterEdgesSubdomains;
+
+    for (int i = 0; i < RobinEdges.size(); ++i)
+    {
+        cout << "Edge number " << RobinEdges[i] << ", subdomain: " << RobinEdgesSubdomains[i] << endl;
+
+    }
+
+
+
+
+
+    
+
+
+    // 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;
+}
+
+double FEM_Matrix::Heat_transfer_coefficient(const int subdomain)
+{
+    int matlab_sd_index = subdomain - 1;
+    double alpha = 0.0;
+
+    switch (matlab_sd_index)
+    {
+        // outside
+        case 0:
+            alpha = 1.0;
+            break;
+
+        // ceramic
+        case 1:
+            alpha = 1.0;
+            break;
+
+        // water
+        case 2:
+            alpha = 1.0;
+            break;
+
+        // air
+        case 3:
+            alpha = 1.0;
+            break;
+
+        default:
+            alpha = 1.0;
+            break;
+    }
+
+    return alpha;
+}
+
+
+
 
 void FEM_Matrix::AddElem_3(int const ial[3], double const ske[3][3], double const fe[3], vector<double> &f)
 {

mgrid_2/getmatrix.h

@@ -394,6 +394,20 @@ class FEM_Matrix: public CRS_Matrix
         */
        void ApplyDirichletBC(std::vector<double> const &u, std::vector<double> &f);
 
+        /**
+        * Applies Robin boundary conditions to stiffness matrix and to load vector @p f
+        * for multiple subdomains
+        * 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);
+
+       double Heat_transfer_coefficient(const int subdomain);
+
        /**
         * Extracts the diagonal elements of the sparse matrix.
         *

mgrid_2/main.cpp

@@ -20,24 +20,24 @@ int main(int argc, char **argv )
 #undef MG
 #ifndef MG
     // Jacobi iteration
-    int nrefine = 0;
+    //int nrefine = 0;
-    if (argc > 1)  nrefine = atoi(argv[1]);
+    //if (argc > 1)  nrefine = atoi(argv[1]);
 
 
     // generating the mesh
     Mesh const mesh_c("../generate_mesh/coffee_cup.txt", "../generate_mesh/coffee_cup_sd.txt");
-    //Mesh const mesh_c("square_tiny.txt");
     bool ba = mesh_c.checkObtuseAngles();
     if (ba) cout << "mesh corrected" << endl;
+    //mesh_c.DebugEdgeBased();
 
-    gMesh_Hierarchy ggm(mesh_c, nrefine);
+    //gMesh_Hierarchy ggm(mesh_c, nrefine);
-    const Mesh &mesh = ggm.finest();
+    //const Mesh &mesh = ggm.finest();
     //mesh.Debug();
     //mesh.DebugEdgeBased();
 
 
     // Initializing FEM matrix !pattern! (only zero entries now)
-    FEM_Matrix SK(mesh);                   // CRS matrix
+    FEM_Matrix SK(mesh_c);                   // CRS matrix
     //SK.writeBinary("sparseMatrix.bin");
     //SK.Debug();
 
@@ -50,19 +50,20 @@ int main(int argc, char **argv )
     //SK.Debug();
 
     // Calculate RHS
-     SK.CalculateRHS(fv, [](double x, double y) {   // rhs
+    //  SK.CalculateRHS(fv, [](double x, double y) {   // rhs
-         return std::sin(M_PI * 2.5 * y) * (M_PI * M_PI * 2.5 * 2.5 * x * x - 2); });
+    //      return std::sin(M_PI * 2.5 * y) * (M_PI * M_PI * 2.5 * 2.5 * x * x - 2); });
+    SK.CalculateRHS(fv, [](double x, double y) {return 0;});
     //SK.CheckRowSum();
     SK.CheckMatrix();
 
     // Initialize temperature
     vector<double> uv(SK.Nrows(), 0.0);    // temperature
-    mesh.Init_Solution_mult(uv, 0, [](double x, double y) -> double { return 18; });    // mug
+    mesh_c.Init_Solution_mult(uv, 0, [](double x, double y) -> double { return 18; });    // mug
-    mesh.Init_Solution_mult(uv, 1, [](double x, double y) -> double { return 80; });    // fluid
+    mesh_c.Init_Solution_mult(uv, 1, [](double x, double y) -> double { return 80; });    // fluid
-    mesh.Init_Solution_mult(uv, 2, [](double x, double y) -> double { return 18; });    // air
+    mesh_c.Init_Solution_mult(uv, 2, [](double x, double y) -> double { return 18; });    // air
 
     // Apply BC
-    SK.ApplyDirichletBC(uv, fv);
+    SK.ApplyRobinBC_mult(uv, fv, 18.0);
 
 
     // Solve
@@ -83,7 +84,7 @@ int main(int argc, char **argv )
     auto [val, idx] = findLargestAbsError(exact_sol, uv, 1e+6, 100);
 
     //mesh.Visualize(getAbsError(exact_sol, uv));
-    mesh.Visualize(uv);
+    mesh_c.Visualize(uv);
 
 #else
     // multigrid iteration