diff --git a/mgrid_2/geom.cpp b/mgrid_2/geom.cpp index 318f9b3..eb2513f 100644 --- a/mgrid_2/geom.cpp +++ b/mgrid_2/geom.cpp @@ -759,6 +759,11 @@ const std::vector Mesh::BoundaryEdges() const return _ebedges; } +const std::vector 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); - PermuteVertices_EdgeBased(perm); + // auto const perm = cuthill_mckee_reordering(_edges); + // PermuteVertices_EdgeBased(perm); } DeriveVertexFromEdgeBased(); diff --git a/mgrid_2/geom.h b/mgrid_2/geom.h index 3dce492..aebbe84 100644 --- a/mgrid_2/geom.h +++ b/mgrid_2/geom.h @@ -331,6 +331,7 @@ private: public: const std::vector BoundaryEdges() const; + const std::vector BoundaryEdgeNodes() const; std::vector OuterEdgesSubdomains; std::vector OuterEdges; diff --git a/mgrid_2/getmatrix.cpp b/mgrid_2/getmatrix.cpp index 17deea4..d20e630 100644 --- a/mgrid_2/getmatrix.cpp +++ b/mgrid_2/getmatrix.cpp @@ -616,57 +616,48 @@ void FEM_Matrix::ApplyDirichletBC(std::vector const &u, std::vector const &u, std::vector &f, const double u_out) +void FEM_Matrix::ApplyRobinBC_mult(std::vector &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 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 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: diff --git a/mgrid_2/getmatrix.h b/mgrid_2/getmatrix.h index eb28851..39cef6f 100644 --- a/mgrid_2/getmatrix.h +++ b/mgrid_2/getmatrix.h @@ -410,11 +410,10 @@ class FEM_Matrix: public CRS_Matrix * The penalty method * 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 const &u, std::vector &f, const double u_out); + void ApplyRobinBC_mult(std::vector &f, const double u_out); double Heat_transfer_coefficient(const int subdomain); diff --git a/mgrid_2/main.cpp b/mgrid_2/main.cpp index 973245a..5414edc 100644 --- a/mgrid_2/main.cpp +++ b/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(t4 - t3); // duration in microseconds