From d7a12583f1915af26ba90ab2aa7a95b0a962ddfc Mon Sep 17 00:00:00 2001 From: Adeeb Arif Kor Date: Tue, 10 Sep 2024 15:35:36 +0800 Subject: [PATCH] Update fenicsx-pc codes for latest DOLFINx --- cpp/fenicsx-pc/common/Lossy.hpp | 102 +++++++++++------- .../examples/lossy_planewave2d_2/main.cpp | 22 +++- cpp/fenicsx/common/Lossy.hpp | 2 +- 3 files changed, 84 insertions(+), 42 deletions(-) diff --git a/cpp/fenicsx-pc/common/Lossy.hpp b/cpp/fenicsx-pc/common/Lossy.hpp index 0c9d7bd..e71375c 100644 --- a/cpp/fenicsx-pc/common/Lossy.hpp +++ b/cpp/fenicsx-pc/common/Lossy.hpp @@ -53,7 +53,8 @@ void axpy(la::Vector& r, T alpha, const la::Vector& x, const la::Vector template class Lossy2D { public: - Lossy2D(std::shared_ptr> Mesh, + Lossy2D(basix::FiniteElement element, + std::shared_ptr> Mesh, std::shared_ptr> FacetTags, std::shared_ptr> speedOfSound, std::shared_ptr> density, @@ -82,7 +83,7 @@ class Lossy2D { // Define function space V = std::make_shared>( - fem::create_functionspace(functionspace_form_forms_a, "u", mesh)); + fem::create_functionspace(mesh, element)); // Define field functions index_map = V->dofmap()->index_map; @@ -103,22 +104,35 @@ class Lossy2D { std::fill(u_.begin(), u_.end(), 1.0); // Compute exterior facets - std::map>>> fd; - auto facet_domains = fem::compute_integration_domains( - fem::IntegralType::exterior_facet, *V->mesh()->topology_mutable(), - ft->indices(), mesh->topology()->dim() - 1, ft->values()); - for (auto& facet : facet_domains) { - std::vector>> x; - x.emplace_back(facet.first, std::span(facet.second.data(), facet.second.size())); - fd.insert({fem::IntegralType::exterior_facet, std::move(x)}); - } + std::vector ft_unique(ft->values().size()); + std::copy(ft->values().begin(), ft->values().end(), ft_unique.begin()); + std::sort(ft_unique.begin(), ft_unique.end()); + auto it = std::unique(ft_unique.begin(), ft_unique.end()); + ft_unique.erase(it, ft_unique.end()); + + std::map>>> fd; + std::map>>> fd_view; + + std::vector facet_domains; + for (auto& tag : ft_unique) { + facet_domains = fem::compute_integration_domains( + fem::IntegralType::exterior_facet, *V->mesh()->topology_mutable(), + ft->find(tag), mesh->topology()->dim()-1); + fd[fem::IntegralType::exterior_facet].push_back( + {tag, facet_domains}); + } + + for (auto const& [key, val] : fd) { + for (auto const& [tag, vec] : val) { + fd_view[key].push_back({tag, std::span(vec.data(), vec.size())}); + } + } // Define LHS form - a = std::make_shared>(fem::create_form( + a = std::make_shared>(fem::create_form( *form_forms_a, {V}, {{"u", u}, {"c0", c0}, {"rho0", rho0}, {"delta0", delta0}}, - {}, fd)); + {}, fd_view, {})); m = std::make_shared>(index_map, bs); m_ = m->mutable_array(); @@ -127,15 +141,15 @@ class Lossy2D { m->scatter_rev(std::plus()); // Define RHS form - L = std::make_shared>( - fem::create_form(*form_forms_L, {V}, + L = std::make_shared>( + fem::create_form(*form_forms_L, {V}, {{"g", g}, {"dg", dg}, {"v_n", v_n}, {"c0", c0}, {"rho0", rho0}, {"delta0", delta0}}, - {}, fd)); + {}, fd_view, {}, {})); b = std::make_shared>(index_map, bs); b_ = b->mutable_array(); @@ -196,11 +210,13 @@ class Lossy2D { dwindow = 0.0; } + /* // Update boundary condition (homogenous domain) - // std::fill(g_.begin(), g_.end(), window * p0 * w0 / s0 * cos(w0 * t)); - // std::fill(dg_.begin(), dg_.end(), - // dwindow * p0 * w0 / s0 * cos(w0 * t) - // - window * p0 * w0 * w0 / s0 * sin(w0 * t)); + std::fill(g_.begin(), g_.end(), window * p0 * w0 / s0 * cos(w0 * t)); + std::fill(dg_.begin(), dg_.end(), + dwindow * p0 * w0 / s0 * cos(w0 * t) + - window * p0 * w0 * w0 / s0 * sin(w0 * t)); + */ // Update boundary condition (heterogenous domain) std::fill(g_.begin(), g_.end(), window * 2.0 * p0 * w0 / s0 * cos(w0 * t)); @@ -375,7 +391,8 @@ class Lossy2D { template class Lossy3D { public: - Lossy3D(std::shared_ptr> Mesh, + Lossy3D(basix::FiniteElement element, + std::shared_ptr> Mesh, std::shared_ptr> FacetTags, std::shared_ptr> speedOfSound, std::shared_ptr> density, @@ -405,7 +422,7 @@ class Lossy3D { // Define function space V = std::make_shared>( - fem::create_functionspace(functionspace_form_forms_a, "u", mesh)); + fem::create_functionspace(mesh, element)); // Define field functions index_map = V->dofmap()->index_map; @@ -426,21 +443,34 @@ class Lossy3D { std::fill(u_.begin(), u_.end(), 1.0); // Compute exterior facets - std::map>>> fd; - auto facet_domains = fem::compute_integration_domains( - fem::IntegralType::exterior_facet, *V->mesh()->topology_mutable(), - ft->indices(), mesh->topology()->dim() - 1, ft->values()); - for (auto& facet : facet_domains) { - std::vector>> x; - x.emplace_back(facet.first, std::span(facet.second.data(), facet.second.size())); - fd.insert({fem::IntegralType::exterior_facet, std::move(x)}); - } + std::vector ft_unique(ft->values().size()); + std::copy(ft->values().begin(), ft->values().end(), ft_unique.begin()); + std::sort(ft_unique.begin(), ft_unique.end()); + auto it = std::unique(ft_unique.begin(), ft_unique.end()); + ft_unique.erase(it, ft_unique.end()); + + std::map>>> fd; + std::map>>> fd_view; + + std::vector facet_domains; + for (auto& tag : ft_unique) { + facet_domains = fem::compute_integration_domains( + fem::IntegralType::exterior_facet, *V->mesh()->topology_mutable(), + ft->find(tag), mesh->topology()->dim()-1); + fd[fem::IntegralType::exterior_facet].push_back( + {tag, facet_domains}); + } + + for (auto const& [key, val] : fd) { + for (auto const& [tag, vec] : val) { + fd_view[key].push_back({tag, std::span(vec.data(), vec.size())}); + } + } // Define LHS form - a = std::make_shared>(fem::create_form( + a = std::make_shared>(fem::create_form( *form_forms_a, {V}, {{"u", u}, {"c0", c0}, {"rho0", rho0}, {"delta0", delta0}}, {}, - fd)); + fd_view, {})); m = std::make_shared>(index_map, bs); m_ = m->mutable_array(); @@ -452,7 +482,7 @@ class Lossy3D { L = std::make_shared>(fem::create_form( *form_forms_L, {V}, {{"g", g}, {"dg", dg}, {"v_n", v_n}, {"c0", c0}, {"rho0", rho0}, {"delta0", delta0}}, {}, - fd)); + fd_view, {}, {})); b = std::make_shared>(index_map, bs); b_ = b->mutable_array(); diff --git a/cpp/fenicsx-pc/examples/lossy_planewave2d_2/main.cpp b/cpp/fenicsx-pc/examples/lossy_planewave2d_2/main.cpp index 7b27786..8b155a0 100644 --- a/cpp/fenicsx-pc/examples/lossy_planewave2d_2/main.cpp +++ b/cpp/fenicsx-pc/examples/lossy_planewave2d_2/main.cpp @@ -51,10 +51,10 @@ int main(int argc, char* argv[]) { const int degreeOfBasis = 4; // Read mesh and mesh tags - auto element = fem::CoordinateElement(mesh::CellType::quadrilateral, 1); + auto coord_element = fem::CoordinateElement(mesh::CellType::quadrilateral, 1); io::XDMFFile fmesh(MPI_COMM_WORLD, "../mesh.xdmf", "r"); auto mesh = std::make_shared>( - fmesh.read_mesh(element, mesh::GhostMode::none, "planewave_2d_4")); + fmesh.read_mesh(coord_element, mesh::GhostMode::none, "planewave_2d_4")); mesh->topology()->create_connectivity(1, 2); auto mt_cell = std::make_shared>( fmesh.read_meshtags(*mesh, "planewave_2d_4_cells")); @@ -75,9 +75,21 @@ int main(int argc, char* argv[]) { MPI_Reduce(&meshSizeMinLocal, &meshSizeMinGlobal, 1, T_MPI, MPI_MIN, 0, MPI_COMM_WORLD); MPI_Bcast(&meshSizeMinGlobal, 1, T_MPI, 0, MPI_COMM_WORLD); + // Finite element + basix::FiniteElement element = basix::create_element( + basix::element::family::P, basix::cell::type::quadrilateral, degreeOfBasis, + basix::element::lagrange_variant::gll_warped, + basix::element::dpc_variant::unset, false + ); + // Define DG function space for the physical parameters of the domain + basix::FiniteElement element_DG = basix::create_element( + basix::element::family::P, basix::cell::type::quadrilateral, 0, + basix::element::lagrange_variant::gll_warped, + basix::element::dpc_variant::unset, true + ); auto V_DG = std::make_shared>( - fem::create_functionspace(functionspace_form_forms_a, "c0", mesh)); + fem::create_functionspace(mesh, element_DG)); auto c0 = std::make_shared>(V_DG); auto rho0 = std::make_shared>(V_DG); auto delta0 = std::make_shared>(V_DG); @@ -129,7 +141,7 @@ int main(int argc, char* argv[]) { } // Model - auto model = Lossy2D(mesh, mt_facet, c0, rho0, delta0, sourceFrequency, + auto model = Lossy2D(element, mesh, mt_facet, c0, rho0, delta0, sourceFrequency, sourceAmplitude, speedOfSound); // Solve @@ -140,7 +152,7 @@ int main(int argc, char* argv[]) { auto u_n = model.u_sol(); // Output to VTX - dolfinx::io::VTXWriter u_out(mesh->comm(), "output_final.bp", {u_n}, "BP4"); + dolfinx::io::VTXWriter u_out(mesh->comm(), "output_final.bp", {u_n}, "bp5"); u_out.write(0.0); } } \ No newline at end of file diff --git a/cpp/fenicsx/common/Lossy.hpp b/cpp/fenicsx/common/Lossy.hpp index 662353a..2d02a56 100644 --- a/cpp/fenicsx/common/Lossy.hpp +++ b/cpp/fenicsx/common/Lossy.hpp @@ -126,7 +126,7 @@ class LossySpectral { } // Define LHS form - a = std::make_shared>(fem::create_form( + a = std::make_shared>(fem::create_form( *form_forms_a, {V}, {{"u", u}, {"c0", c0}, {"rho0", rho0}, {"delta0", delta0}}, {}, fd_view, {}));