From 99b46ad7ec92285527d833118250c33fe64be4ee Mon Sep 17 00:00:00 2001 From: "jakob.schratter" Date: Sat, 24 Jan 2026 14:26:22 +0100 Subject: [PATCH] added subdomain indexing for outer edges, for use in robin BC --- mgrid_2/geom.cpp | 28 +++++++------ mgrid_2/geom.h | 9 ++-- mgrid_2/getmatrix.cpp | 98 +++++++++++++++++++++++++++++++++++++++++-- mgrid_2/getmatrix.h | 14 +++++++ mgrid_2/main.cpp | 27 ++++++------ 5 files changed, 143 insertions(+), 33 deletions(-) diff --git a/mgrid_2/geom.cpp b/mgrid_2/geom.cpp index 6bc1482..318f9b3 100644 --- a/mgrid_2/geom.cpp +++ b/mgrid_2/geom.cpp @@ -761,24 +761,25 @@ const std::vector Mesh::BoundaryEdges() const // Only the outer edges for Robin BC -const std::vector Mesh::OuterEdges() const +void Mesh::InitializeOuterEdges() { - vector outerEdges; - - for (int k = 0; k < _ebedges.size()/2; ++k) + assert(EdgeSubdomains.size() == 2*_ebedges.size()); + for (size_t k = 0; k < _ebedges.size(); ++k) // loop over all boundary edges { - if (EdgeSubdomains[k] == 0 || EdgeSubdomains[k + 1] == 0) - outerEdges.push_back(_ebedges[k]); + if (EdgeSubdomains[2*k] == 0) + { + 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 << _ebedges; - // cout << endl; - // cout << outerEdges << endl; - - return outerEdges; + cout << "Outer boundary edges: " << OuterEdges.size() << endl; } @@ -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(); } diff --git a/mgrid_2/geom.h b/mgrid_2/geom.h index ac7787f..3dce492 100644 --- a/mgrid_2/geom.h +++ b/mgrid_2/geom.h @@ -326,10 +326,15 @@ public: private: void Write_ascii_paraview_2D(std::string const &fname, std::vector const &v) const; void Write_ascii_paraview_3D(std::string const &fname, std::vector const &v) const; + + void InitializeOuterEdges(); public: const std::vector BoundaryEdges() const; - const std::vector OuterEdges() const; + + std::vector OuterEdgesSubdomains; + std::vector OuterEdges; + /** @@ -565,8 +570,6 @@ public: */ [[nodiscard]] const std::vector ReadElementSubdomains(std::string const &dname) const; - [[nodiscard]] const std::vector ReadEdgeSubdomains(std::string const &filename) const; - /** * Calculates the largest inner angle in element @p idx. diff --git a/mgrid_2/getmatrix.cpp b/mgrid_2/getmatrix.cpp index 40fe2eb..bd00fb1 100644 --- a/mgrid_2/getmatrix.cpp +++ b/mgrid_2/getmatrix.cpp @@ -412,7 +412,7 @@ void FEM_Matrix::CalculateLaplaceMult(vector &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 &f) cout << "finished in " << duration << " sec. ########\n"; // ToDo: change to systemclock //Debug(); - - - return; } @@ -583,6 +580,99 @@ void FEM_Matrix::ApplyDirichletBC(std::vector const &u, std::vector const &u, 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) + { + cout << "Edge number " << RobinEdges[i] << ", subdomain: " << RobinEdgesSubdomains[i] << endl; + + } + + + + + + + + + // 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; +} + +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 &f) { diff --git a/mgrid_2/getmatrix.h b/mgrid_2/getmatrix.h index 7b39b08..6761ef3 100644 --- a/mgrid_2/getmatrix.h +++ b/mgrid_2/getmatrix.h @@ -394,6 +394,20 @@ class FEM_Matrix: public CRS_Matrix */ void ApplyDirichletBC(std::vector const &u, std::vector &f); + /** + * Applies Robin boundary conditions to stiffness matrix and to load vector @p f + * for multiple subdomains + * 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); + + double Heat_transfer_coefficient(const int subdomain); + /** * Extracts the diagonal elements of the sparse matrix. * diff --git a/mgrid_2/main.cpp b/mgrid_2/main.cpp index 79d2ca7..576b0a9 100644 --- a/mgrid_2/main.cpp +++ b/mgrid_2/main.cpp @@ -20,24 +20,24 @@ int main(int argc, char **argv ) #undef MG #ifndef MG // Jacobi iteration - int nrefine = 0; - if (argc > 1) nrefine = atoi(argv[1]); + //int nrefine = 0; + //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); - const Mesh &mesh = ggm.finest(); + //gMesh_Hierarchy ggm(mesh_c, nrefine); + //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 - 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) { // rhs + // 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 uv(SK.Nrows(), 0.0); // temperature - mesh.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.Init_Solution_mult(uv, 2, [](double x, double y) -> double { return 18; }); // air + mesh_c.Init_Solution_mult(uv, 0, [](double x, double y) -> double { return 18; }); // mug + mesh_c.Init_Solution_mult(uv, 1, [](double x, double y) -> double { return 80; }); // fluid + 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