From 7f442457ca7d27c568248df7a2c729282dd8507d Mon Sep 17 00:00:00 2001 From: Adeeb Arif Kor Date: Tue, 10 Sep 2024 17:10:10 +0800 Subject: [PATCH] Update fenicsx-sf codes for latest DOLFINx --- cpp/fenicsx-sf/common/Linear.hpp | 44 ++++++++++++++------- cpp/fenicsx-sf/common/Lossy.hpp | 58 ++++++++++++++++++---------- cpp/fenicsx-sf/common/Westervelt.hpp | 56 +++++++++++++++++---------- 3 files changed, 102 insertions(+), 56 deletions(-) diff --git a/cpp/fenicsx-sf/common/Linear.hpp b/cpp/fenicsx-sf/common/Linear.hpp index fdd7f06..f8cef6a 100644 --- a/cpp/fenicsx-sf/common/Linear.hpp +++ b/cpp/fenicsx-sf/common/Linear.hpp @@ -52,7 +52,8 @@ void axpy(la::Vector& r, T alpha, const la::Vector& x, const la::Vector template class LinearSpectral3D { public: - LinearSpectral3D(std::shared_ptr> Mesh, + LinearSpectral3D(basix::FiniteElement element, + std::shared_ptr> Mesh, std::shared_ptr> FacetTags, std::shared_ptr> speedOfSound, std::shared_ptr> density, @@ -79,7 +80,7 @@ class LinearSpectral3D { // 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; @@ -98,19 +99,32 @@ class LinearSpectral3D { 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}}, {}, {})); m = std::make_shared>(index_map, bs); @@ -120,10 +134,10 @@ class LinearSpectral3D { m->scatter_rev(std::plus()); // Define RHS form - L = std::make_shared>(fem::create_form( + L = std::make_shared>(fem::create_form( *form_forms_L, {V}, {{"g", g}, {"v_n", v_n}, {"c0", c0}, {"rho0", rho0}}, {}, - fd)); + fd_view, {}, {})); b = std::make_shared>(index_map, bs); b_ = b->mutable_array(); diff --git a/cpp/fenicsx-sf/common/Lossy.hpp b/cpp/fenicsx-sf/common/Lossy.hpp index 3cdd79c..cf13882 100644 --- a/cpp/fenicsx-sf/common/Lossy.hpp +++ b/cpp/fenicsx-sf/common/Lossy.hpp @@ -39,7 +39,6 @@ void axpy(la::Vector& r, T alpha, const la::Vector& x, const la::Vector } // namespace kernels - /// Solver for the 3D second order linear wave equation with attenuation. /// This solver uses GLL lattice and GLL quadrature such that it produces /// a diagonal mass matrix. @@ -54,7 +53,8 @@ void axpy(la::Vector& r, T alpha, const la::Vector& x, const la::Vector template class LossySpectral3D { public: - LossySpectral3D(std::shared_ptr> Mesh, + LossySpectral3D(basix::FiniteElement element, + std::shared_ptr> Mesh, std::shared_ptr> FacetTags, std::shared_ptr> speedOfSound, std::shared_ptr> density, @@ -62,6 +62,7 @@ class LossySpectral3D { const T& sourceFrequency, const T& sourceAmplitude, const T& sourceSpeed) { + // MPI MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); @@ -83,7 +84,7 @@ class LossySpectral3D { // 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; @@ -104,21 +105,34 @@ class LossySpectral3D { 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,10 +141,10 @@ class LossySpectral3D { m->scatter_rev(std::plus()); // Define RHS form - L = std::make_shared>(fem::create_form( + 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(); @@ -191,11 +205,13 @@ class LossySpectral3D { 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)); diff --git a/cpp/fenicsx-sf/common/Westervelt.hpp b/cpp/fenicsx-sf/common/Westervelt.hpp index 946f7fb..ee317ba 100644 --- a/cpp/fenicsx-sf/common/Westervelt.hpp +++ b/cpp/fenicsx-sf/common/Westervelt.hpp @@ -55,7 +55,8 @@ void axpy(la::Vector& r, T alpha, const la::Vector& x, const la::Vector template class WesterveltSpectral3D { public: - WesterveltSpectral3D(std::shared_ptr> Mesh, + WesterveltSpectral3D(basix::FiniteElement element, + std::shared_ptr> Mesh, std::shared_ptr> FacetTags, std::shared_ptr> speedOfSound, std::shared_ptr> density, @@ -86,7 +87,7 @@ class WesterveltSpectral3D { // 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; @@ -108,21 +109,34 @@ class WesterveltSpectral3D { 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, {})); m0 = std::make_shared>(index_map, bs); m0_ = m0->mutable_array(); @@ -133,10 +147,10 @@ class WesterveltSpectral3D { m_ = m->mutable_array(); // Define RHS form - L = std::make_shared>(fem::create_form( + 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(); @@ -211,11 +225,13 @@ class WesterveltSpectral3D { 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));