fixed robin BC

generate_mesh/coffee_cup.m

@@ -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',0.1); 
+[p,e,t] = initmesh(g,'hmax',0.0015); 
 pdemesh(p,e,t)
 
 

mgrid_2/geom.cpp

@@ -775,11 +775,15 @@ void Mesh::InitializeOuterEdges()
         {
             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)
         {            
             OuterEdges.push_back(_ebedges[k]);
             OuterEdgesSubdomains.push_back(EdgeSubdomains[2*k]);
+            OuterEdgesNodes.push_back(_bedges[2*k]);
+            OuterEdgesNodes.push_back(_bedges[2*k + 1]);
         }
     }
 
@@ -1014,12 +1018,12 @@ Mesh::Mesh(std::string const &fname)
     }
     DeriveVertexFromEdgeBased();    
     
-     //// GH: Check permuted numbering
+    //  //// GH: Check permuted numbering
-    //vector<int> perm(Nnodes());
+    // vector<int> perm(Nnodes());
-    //iota(rbegin(perm),rend(perm),0);
+    // iota(rbegin(perm),rend(perm),0);
-    //random_shuffle(begin(perm),end(perm));
+    // random_shuffle(begin(perm),end(perm));
-    //PermuteVertices(perm);   
+    // PermuteVertices(perm);   
-    //cout << " P E R M U T E D !" << endl;
+    // cout << " P E R M U T E D !" << endl;
 }
 
 

mgrid_2/geom.h

@@ -333,8 +333,9 @@ public:
     const std::vector<int> BoundaryEdges() const;
     const std::vector<int> BoundaryEdgeNodes() const;
 
-    std::vector<int> OuterEdgesSubdomains;
     std::vector<int> OuterEdges;
+    std::vector<int> OuterEdgesSubdomains;
+    std::vector<int> OuterEdgesNodes;
 
 
 

mgrid_2/getmatrix.cpp

@@ -598,11 +598,11 @@ void FEM_Matrix::ApplyDirichletBC(std::vector<double> const &u, std::vector<doub
         int const row = idx[i];
         for (int ij = _id[row]; ij < _id[row + 1]; ++ij) {
             int const col = _ik[ij];
-            if (col == row) {
+            if (col == row) { // then col = _ik[_id[row]]
                 _sk[ij] = 1.0;
                 f[row]  = u[row];
             }
-            else {
+            else {  // col != row, then 
                 int const id1 = fetch(col, row); // Find entry (col,row)
                 assert(id1 >= 0);
                 f[col] -= _sk[id1] * u[row];
@@ -620,43 +620,46 @@ void FEM_Matrix::ApplyRobinBC_mult(std::vector<double> &f, const double u_out)
 {
     auto const RobinEdges = _mesh.OuterEdges;
     auto const RobinEdgesSubdomains = _mesh.OuterEdgesSubdomains;
+    auto const RobinEdgeNodes = _mesh.OuterEdgesNodes;
+    assert(RobinEdgeNodes.size() == 2 * RobinEdges.size());
 
-    auto const BoundaryEdges = _mesh.BoundaryEdges();
+    vector<double> const Coordinates = _mesh.GetCoords();
-    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;
+        int const subdomain = RobinEdgesSubdomains[i];
-        double alpha = Heat_transfer_coefficient(RobinEdgesSubdomains[i]);
+        double const alpha = Heat_transfer_coefficient(subdomain);
 
+        //int const EdgeNumber = RobinEdges[i];
 
-        int const EdgeNode1 = BoundaryEdgeNodes[2*i];
+        int const EdgeNode1 = RobinEdgeNodes[2*i];
-        int const EdgeNode2 = BoundaryEdgeNodes[2*i + 1];
+        int const EdgeNode2 = RobinEdgeNodes[2*i + 1];
+        
+        //cout << "Edge number " << EdgeNumber << ", subdomain: " << subdomain << "        ";
+        //cout << "Node1: " << EdgeNode1 << "     Node2: " << EdgeNode2 << "            " << endl;
+
+        double x_1 = Coordinates[2*EdgeNode1];
+        double y_1 = Coordinates[2*EdgeNode1 + 1];
+        double x_2 = Coordinates[2*EdgeNode2];
+        double y_2 = Coordinates[2*EdgeNode2 + 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));
+        //cout << EdgeLength << endl;
 
+        int ii = fetch(EdgeNode1, EdgeNode1);
+        int jj = fetch(EdgeNode2, EdgeNode2);
+        int ij = fetch(EdgeNode1, EdgeNode2);
+        int ji = fetch(EdgeNode2, EdgeNode1);
 
-
+        cout << "ii: " << ii << ", jj: " << jj << ", ij: " << ij << ", ji: " << ji << endl;
-        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[EdgeNode1] += EdgeLength*alpha*u_out/2;
-        f[jj] += EdgeLength*alpha*u_out/2;
+        f[EdgeNode2] += EdgeLength*alpha*u_out/2;
 
     }
 

mgrid_2/main.cpp

@@ -45,14 +45,14 @@ int main(int argc, char **argv )
     // Initializing FEM matrix !pattern! (only zero entries now)
     FEM_Matrix SK(mesh_c);                                    // CRS matrix
     //SK.writeBinary("sparseMatrix.bin");
-    // SK.Debug();
+    //SK.Debug();
 
     vector<double> fv(SK.Nrows(), 0.0);                    
     SK.CalculateRHS(fv, [](double x, double y) {return 0;});  // rhs (f=0)
 
     SK.CalculateLaplace_mult(fv);                             // stiffness matrix (+K)
     SK.AddMass_mult(fv, 1.0/dt);                              // add mass matrix (+M/dt)
-    // SK.ApplyRobinBC_mult(fv, u_out);                           // apply Robin-bnd (+C = +F)
+    SK.ApplyRobinBC_mult(fv, u_out);                           // apply Robin-bnd (+C = +F)
                                                               // SK = M/dt + K + C = F
     
     // SK.Debug();