Skip to content

Commit

Permalink
Update fenicsx-pc codes for latest DOLFINx
Browse files Browse the repository at this point in the history
  • Loading branch information
@adeebkor
adeebkor committed Sep 10, 2024
1 parent 41e8362 commit d7a1258
Show file tree
Hide file tree
Showing 3 changed files with 84 additions and 42 deletions.
102 changes: 66 additions & 36 deletions cpp/fenicsx-pc/common/Lossy.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,8 @@ void axpy(la::Vector<T>& r, T alpha, const la::Vector<T>& x, const la::Vector<T>
template <typename T, int P, int Q>
class Lossy2D {
public:
Lossy2D(std::shared_ptr<mesh::Mesh<T>> Mesh,
Lossy2D(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,
Expand Down Expand Up @@ -82,7 +83,7 @@ class Lossy2D {

// 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;
Expand All @@ -103,22 +104,35 @@ class Lossy2D {
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, T>(
*form_forms_a, {V},
{{"u", u}, {"c0", c0}, {"rho0", rho0}, {"delta0", delta0}},
{}, fd));
{}, fd_view, {}));

m = std::make_shared<la::Vector<T>>(index_map, bs);
m_ = m->mutable_array();
Expand All @@ -127,15 +141,15 @@ class Lossy2D {
m->scatter_rev(std::plus<T>());

// Define RHS form
L = std::make_shared<fem::Form<T, T>>(
fem::create_form<T, T>(*form_forms_L, {V},
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();

Expand Down Expand Up @@ -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));
Expand Down Expand Up @@ -375,7 +391,8 @@ class Lossy2D {
template <typename T, int P, int Q>
class Lossy3D {
public:
Lossy3D(std::shared_ptr<mesh::Mesh<T>> Mesh,
Lossy3D(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,
Expand Down Expand Up @@ -405,7 +422,7 @@ class Lossy3D {

// 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;
Expand All @@ -426,21 +443,34 @@ class Lossy3D {
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, T>(
*form_forms_a, {V}, {{"u", u}, {"c0", c0}, {"rho0", rho0}, {"delta0", delta0}}, {},
fd));
fd_view, {}));

m = std::make_shared<la::Vector<T>>(index_map, bs);
m_ = m->mutable_array();
Expand All @@ -452,7 +482,7 @@ class Lossy3D {
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();

Expand Down
22 changes: 17 additions & 5 deletions cpp/fenicsx-pc/examples/lossy_planewave2d_2/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,10 +51,10 @@ int main(int argc, char* argv[]) {
const int degreeOfBasis = 4;

// Read mesh and mesh tags
auto element = fem::CoordinateElement<T>(mesh::CellType::quadrilateral, 1);
auto coord_element = fem::CoordinateElement<T>(mesh::CellType::quadrilateral, 1);
io::XDMFFile fmesh(MPI_COMM_WORLD, "../mesh.xdmf", "r");
auto mesh = std::make_shared<mesh::Mesh<T>>(
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<mesh::MeshTags<std::int32_t>>(
fmesh.read_meshtags(*mesh, "planewave_2d_4_cells"));
Expand All @@ -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<T>(
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<T>(
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::FunctionSpace<T>>(
fem::create_functionspace(functionspace_form_forms_a, "c0", mesh));
fem::create_functionspace(mesh, element_DG));
auto c0 = std::make_shared<fem::Function<T>>(V_DG);
auto rho0 = std::make_shared<fem::Function<T>>(V_DG);
auto delta0 = std::make_shared<fem::Function<T>>(V_DG);
Expand Down Expand Up @@ -129,7 +141,7 @@ int main(int argc, char* argv[]) {
}

// Model
auto model = Lossy2D<T, 4, 5>(mesh, mt_facet, c0, rho0, delta0, sourceFrequency,
auto model = Lossy2D<T, 4, 5>(element, mesh, mt_facet, c0, rho0, delta0, sourceFrequency,
sourceAmplitude, speedOfSound);

// Solve
Expand All @@ -140,7 +152,7 @@ int main(int argc, char* argv[]) {
auto u_n = model.u_sol();

// Output to VTX
dolfinx::io::VTXWriter<T> u_out(mesh->comm(), "output_final.bp", {u_n}, "BP4");
dolfinx::io::VTXWriter<T> u_out(mesh->comm(), "output_final.bp", {u_n}, "bp5");
u_out.write(0.0);
}
}
2 changes: 1 addition & 1 deletion cpp/fenicsx/common/Lossy.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,7 +126,7 @@ class LossySpectral {
}

// 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, T>(
*form_forms_a, {V},
{{"u", u}, {"c0", c0}, {"rho0", rho0}, {"delta0", delta0}},
{}, fd_view, {}));
Expand Down

0 comments on commit d7a1258

Please sign in to comment.