test_emi_dof_numbering.cc

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: Apache-2.0
 
#include "SF_base.h"
#include "basics.h"
 
#include <mpi.h>
 
#include <iostream>
#include <map>
#include <string>
#include <utility>
#include <vector>
 
namespace {
 
using mesh_t = SF::meshdata<mesh_int_t, mesh_real_t>;
using vertex_tag = std::pair<mesh_int_t, mesh_int_t>;
 
constexpr mesh_int_t MYOCYTE = 1;
constexpr mesh_int_t EXTRA = 2;
 
int failures = 0;
 
void check(bool ok, const std::string& what)
{
  if (!ok) {
    std::cerr << "FAILED: " << what << '\n';
    failures++;
  }
}
 
void build_mesh(int rank, int size, MPI_Comm comm, mesh_t& mesh)
{
  const std::vector<mesh_int_t> myocyte{0, 1, 2, 3}, extra{1, 2, 3, 5};
  std::vector<std::pair<mesh_int_t, const std::vector<mesh_int_t>*>> elems;
  if (rank == 0) elems.push_back({MYOCYTE, &myocyte});
  if (rank == size - 1) elems.push_back({EXTRA, &extra});
 
  mesh.comm = comm;
  SF::vector<mesh_int_t> cnt;
  for (const auto& [tag, nodes] : elems) {
    for (mesh_int_t n : *nodes) mesh.con.push_back(n);
    mesh.tag.push_back(tag);
    mesh.type.push_back(SF::Tetra);
    cnt.push_back(4);
  }
  mesh.l_numelem = elems.size();
  mesh.g_numelem = 2;
  SF::dsp_from_cnt(cnt, mesh.dsp);
  SF::vector<mesh_int_t>& eref = mesh.register_numbering(SF::NBR_ELEM_REF);
  eref.resize(mesh.l_numelem);
  for (size_t i = 0; i < eref.size(); i++) eref[i] = (rank == 0 ? 0 : 1) + i;
  mesh.localize(SF::NBR_REF);
  // as emi.cc before decoupling: g_numpts counts the used points
  mesh.generate_par_layout();
}
 
}  // namespace
 
int main(int argc, char** argv)
{
  MPI_Init(&argc, &argv);
  MPI_Comm comm = MPI_COMM_WORLD;
  int size, rank;
  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);
 
  mesh_t mesh;
  build_mesh(rank, size, comm, mesh);
  check(mesh.g_numpts == 5, "g_numpts counts the 5 used points");
 
  // 1 on the membrane, 0 inside a region
  hashmap::unordered_map<mesh_int_t, mesh_int_t> ptsdata{{0, 0}, {1, 1}, {2, 1}, {3, 1}, {5, 0}};
  hashmap::unordered_set<int> extra_tags{EXTRA};
  hashmap::unordered_map<vertex_tag, mesh_int_t> dofs;
  SF::compute_map_vertex_to_dof(mesh, SF::NBR_REF, ptsdata, extra_tags, dofs);
 
  // gather (vertex, tag, dof) on rank 0
  std::vector<mesh_int_t> local;
  for (const auto& [key, dof] : dofs) local.insert(local.end(), {key.first, key.second, dof});
  int n = local.size();
  std::vector<int> counts(size), displs(size);
  MPI_Gather(&n, 1, MPI_INT, counts.data(), 1, MPI_INT, 0, comm);
  for (int r = 1; r < size; r++) displs[r] = displs[r - 1] + counts[r - 1];
  std::vector<mesh_int_t> all(rank == 0 ? displs[size - 1] + counts[size - 1] : 0);
  MPI_Datatype type = sizeof(mesh_int_t) == 8 ? MPI_LONG_LONG : MPI_INT;
  MPI_Gatherv(local.data(), n, type, all.data(), counts.data(), displs.data(), type, 0, comm);
 
  if (rank == 0) {
    std::map<vertex_tag, mesh_int_t> dof_of;
    for (size_t i = 0; i < all.size(); i += 3) {
      const vertex_tag key{all[i], all[i + 1]};
      auto [it, inserted] = dof_of.insert({key, all[i + 2]});
      check(inserted || it->second == all[i + 2], "one DOF per (vertex, tag)");
    }
    for (mesh_int_t v : {1, 2, 3}) {
      check(dof_of.count({v, MYOCYTE}) && dof_of.count({v, EXTRA}),
            "membrane vertex " + std::to_string(v) + " has a myocyte and an extracellular DOF");
      check(dof_of[{v, MYOCYTE}] != dof_of[{v, EXTRA}],
            "membrane vertex " + std::to_string(v) + " is decoupled");
    }
    std::map<mesh_int_t, mesh_int_t> vertex_of;
    for (const auto& [key, dof] : dof_of) {
      auto [it, inserted] = vertex_of.insert({dof, key.first});
      check(inserted || it->second == key.first,
            "DOF " + std::to_string(dof) + " belongs to vertices " + std::to_string(it->second) +
                " and " + std::to_string(key.first));
    }
  }
 
  int total = 0;
  MPI_Allreduce(&failures, &total, 1, MPI_INT, MPI_SUM, comm);
  MPI_Finalize();
  return total ? 1 : 0;
}