From 0a3e5b6ff54ac187bf5fdf6b97d4201a21a22e0a Mon Sep 17 00:00:00 2001
From: Adeeb Arif Kor <adeebkor@um.edu.my>
Date: Tue, 10 Sep 2024 17:12:40 +0800
Subject: [PATCH] Update fenicsx-sf-fastor codes for latest DOLFINx

---
 cpp/fenicsx-sf-fastor/common/Westervelt.hpp | 56 +++++++++++++--------
 1 file changed, 36 insertions(+), 20 deletions(-)

diff --git a/cpp/fenicsx-sf-fastor/common/Westervelt.hpp b/cpp/fenicsx-sf-fastor/common/Westervelt.hpp
index 946f7fb..ee317ba 100644
--- a/cpp/fenicsx-sf-fastor/common/Westervelt.hpp
+++ b/cpp/fenicsx-sf-fastor/common/Westervelt.hpp
@@ -55,7 +55,8 @@ void axpy(la::Vector<T>& r, T alpha, const la::Vector<T>& x, const la::Vector<T>
 template <typename T, int P>
 class WesterveltSpectral3D {
 public:
-  WesterveltSpectral3D(std::shared_ptr<mesh::Mesh<T>> Mesh,
+  WesterveltSpectral3D(basix::FiniteElement<T> element,
+                       std::shared_ptr<mesh::Mesh<T>> Mesh,
                        std::shared_ptr<mesh::MeshTags<std::int32_t>> FacetTags,
                        std::shared_ptr<fem::Function<T>> speedOfSound,
                        std::shared_ptr<fem::Function<T>> density,
@@ -86,7 +87,7 @@ class WesterveltSpectral3D {
 
     // Define function space
     V = std::make_shared<fem::FunctionSpace<T>>(
-        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<fem::IntegralType,
-        std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>> 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<std::pair<std::int32_t, std::span<const std::int32_t>>> 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<std::int32_t> 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<fem::IntegralType, std::vector<std::pair<std::int32_t, std::vector<std::int32_t>>>> fd;
+    std::map<fem::IntegralType, std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>> fd_view;
+
+    std::vector<std::int32_t> 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::Form<T, T>>(fem::create_form<T, T>(
+    a = std::make_shared<fem::Form<T>>(fem::create_form<T>(
         *form_forms_a, {V}, {{"u", u}, {"c0", c0}, {"rho0", rho0}, {"delta0", delta0}}, {},
-        fd));
+        fd_view, {}));
 
     m0 = std::make_shared<la::Vector<T>>(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::Form<T, T>>(fem::create_form<T, T>(
+    L = std::make_shared<fem::Form<T>>(fem::create_form<T>(
         *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<la::Vector<T>>(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));