fem_utils.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------

#ifndef _FEM_UTILS_H
#define _FEM_UTILS_H

#include "sf_interface.h"

#include "petsc_utils.h" // TODO: for EXIT

namespace opencarp {


#define COMMENT_CHAR '#'

void parse_comment_line(char* buff, const int buffsize, std::map<std::string,std::string> & metadata);

void read_metadata(const std::string filename, std::map<std::string,std::string> & metadata, MPI_Comm comm);

char* skip_comments(FILE_SPEC stream, char* readbuff, size_t buffsize, MPI_Comm comm);

template<class T> inline
void read_indices_global(SF::vector<T> & idx,
                         const std::string filename,
                         MPI_Comm comm)
{
  int numread = 0, err = 0;

  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);

  FILE_SPEC stream = NULL;
  const size_t str_size = 5000;
  char readbuff[str_size];

  if(rank == 0)
    stream = f_open(filename.c_str(), "r");

  // skip comments
  char* ptr = skip_comments(stream, readbuff, str_size, comm);

  long int num_inp_nod;
  if(ptr != NULL) {
    numread = sscanf(readbuff, "%ld", &num_inp_nod);
  }

  if(numread == 0) {
    log_msg(0, 5, 0, "%s error: Could not determine vtx data size! Aborting!\n", __func__);
    EXIT(1);
  }

  SF::vector<int> nodbuff(num_inp_nod, -1);

  if(rank == 0) {
    long int num_read_loc = 0;

    ptr = fgets(readbuff, str_size, stream->fd);
    while(ptr != NULL) {
      int n = sscanf(readbuff, "%d", &nodbuff[num_read_loc]);

      // only increment local read counter if we indeed read an int
      if(n == 1) num_read_loc++;

      if(num_read_loc < num_inp_nod)
        ptr = fgets(readbuff, str_size, stream->fd);
      else
        ptr = NULL;
    }
  }

  f_close(stream);

  MPI_Bcast(nodbuff.data(), nodbuff.size(), MPI_INT, 0, comm);

  idx.reserve(nodbuff.size());
  for(int n : nodbuff)
    if(n > -1) idx.push_back(n);
}

template<class T> inline
void read_indices(SF::vector<T> & idx,
                  const std::string filename,
                  const hashmap::unordered_map<mesh_int_t, mesh_int_t> & dd_map,
                  MPI_Comm comm)
{
  int numread = 0, err = 0;

  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);

  // generate a hashed set of the nodes in the local domain
  hashmap::unordered_set<mesh_int_t> have_read;

  idx.resize(0);
  idx.reserve(dd_map.size());

  FILE_SPEC stream = NULL;

  if(rank == 0)
    stream = f_open(filename.c_str(), "r");

  const size_t str_size = 5000;
  char readbuff[str_size];

  // skip comments
  char* ptr = skip_comments(stream, readbuff, str_size, comm);

  long int num_inp_nod;
  if(ptr != NULL) {
    numread = sscanf(readbuff, "%ld", &num_inp_nod);
  }

  if(numread == 0) {
    log_msg(0, 5, 0, "%s error: Could not determine vtx data size! Aborting!\n", __func__);
    EXIT(1);
  }

  int nodbuff_size = 50000;
  if(nodbuff_size > num_inp_nod) nodbuff_size = num_inp_nod;

  SF::vector<int> nodbuff;
  numread = 0;

  while(numread < num_inp_nod) {
    nodbuff.assign(size_t(nodbuff_size), int(-1));

    if(rank == 0) {
      int num_read_loc = 0;

      ptr = fgets(readbuff, str_size, stream->fd);
      while(ptr != NULL) {
        int n = sscanf(readbuff, "%d", &nodbuff[num_read_loc]);
        // only increment local read counter if we indeed read an int
        if(n == 1) num_read_loc++;

        if(num_read_loc < nodbuff_size)
          ptr = fgets(readbuff, str_size, stream->fd);
        else
          ptr = NULL;
      }
    }

    MPI_Bcast(nodbuff.data(), nodbuff_size*sizeof(int), MPI_BYTE, 0, comm);

    for(int ridx = 0; ridx < nodbuff_size; ridx++) {
      int n = nodbuff[ridx];
      auto it = dd_map.find(n);

      if(n > -1 && it != dd_map.end() && have_read.count(n) == 0) {
        have_read.insert(n);
        mesh_int_t ln = it->second;
        idx.push_back(ln);
      }
      numread++;
    }
  }

  f_close(stream);
}

template<class T> inline
void read_indices(SF::vector<T> & idx,
                  const std::string filename,
                  const sf_mesh & mesh,
                  const SF::SF_nbr nbr,
                  const bool algebraic,
                  MPI_Comm comm)
{
  // generate a hashed set of the nodes in the local domain
  hashmap::unordered_map<mesh_int_t, mesh_int_t> dd_map;
  const SF::vector<mesh_int_t> & dd_nbr = mesh.get_numbering(nbr);

  if(algebraic) {
    const SF::vector<mesh_int_t> & alg_nod = mesh.pl.algebraic_nodes();
    for(size_t i = 0; i<alg_nod.size(); i++) {
      mesh_int_t an = alg_nod[i];
      dd_map[dd_nbr[an]] = an;
    }
  } else {
    for(size_t i = 0; i<dd_nbr.size(); i++)
      dd_map[dd_nbr[i]] = i;
  }

  read_indices(idx, filename, dd_map, comm);
}

template<class T> inline
void read_indices(SF::vector<T> & idx,
                  const std::string filename,
                  const SF::vector<mesh_int_t> & dd_nbr,
                  MPI_Comm comm)
{
  // generate a hashed set of the nodes in the local domain
  hashmap::unordered_map<mesh_int_t, mesh_int_t> dd_map;
  for(size_t i = 0; i<dd_nbr.size(); i++)
    dd_map[dd_nbr[i]] = i;

  read_indices(idx, filename, dd_map, comm);
}

template<class T, class S> inline
void read_indices_with_data(SF::vector<T> & idx, SF::vector<S> & dat,
                            const std::string filename,
                            const hashmap::unordered_map<mesh_int_t, mesh_int_t> & dd_map,
                            const int dpn,
                            MPI_Comm comm)
{
  int numread = 0, err = 0;

  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);

  // Currenty we support dpn <= 3
  if(dpn > 3) {
    if(rank == 0) {
      fprintf(stderr, "%s error: dpn > 3 not supported!\n", __func__);
    }
    return;
  }

  // generate a hashed set of the nodes in the local domain
  hashmap::unordered_set<mesh_int_t> have_read;

  idx.resize(dd_map.size());
  dat.resize(dd_map.size() * dpn);
  FILE_SPEC stream = NULL;

  if(rank == 0)
    stream = f_open(filename.c_str(), "r");

  const size_t str_size = 5000;
  char readbuff[str_size];

  // skip comments
  char* ptr = skip_comments(stream, readbuff, str_size, comm);

  // parse size of vtx file. for simplicity we parse on each rank.
  long int num_inp_nod;
  if(ptr != NULL) {
    numread = sscanf(readbuff, "%ld", &num_inp_nod);
  }

  if(numread == 0) {
    log_msg(0, 5, 0, "%s error: Could not determine vtx data size! Aborting!\n", __func__);
    EXIT(1);
  }

  int nodbuff_size = 50000;
  if(nodbuff_size > num_inp_nod) nodbuff_size = num_inp_nod;

  SF::vector<int>    nodbuff;
  SF::vector<double> doublebuff;

  size_t widx = 0;
  numread = 0;

  while(numread < num_inp_nod) {
    nodbuff   .assign(size_t(nodbuff_size), int(-1));
    doublebuff.resize(size_t(nodbuff_size * dpn));

    if(rank == 0) {
      int num_read_loc = 0;

      ptr = fgets(readbuff, str_size, stream->fd);
      while(ptr != NULL) {
        int n = sscanf(readbuff, "%d %lf %lf %lf", &nodbuff[num_read_loc],
                       &doublebuff[num_read_loc*dpn+0], &doublebuff[num_read_loc*dpn+1], &doublebuff[num_read_loc*dpn+2]);

        // only increment local read counter if we indeed read an int
        if(n == dpn+1) num_read_loc++;

        if(num_read_loc < nodbuff_size)
          ptr = fgets(readbuff, str_size, stream->fd);
        else
          ptr = NULL;
      }
    }

    MPI_Bcast(nodbuff.data(),    nodbuff_size,     MPI_INT,    0, comm);
    MPI_Bcast(doublebuff.data(), nodbuff_size*dpn, MPI_DOUBLE, 0, comm);

    for(int ridx = 0; ridx < nodbuff_size; ridx++) {
      int n = nodbuff[ridx];

      auto it = dd_map.find(n);
      if(n > -1 && it != dd_map.end() && have_read.count(n) == 0) {
        have_read.insert(n);
        idx[widx] = it->second;

        for(int j=0; j<dpn; j++)
          dat[widx*dpn+j] = doublebuff[ridx*dpn+j];

        widx++;
      }

      numread++;
    }
  }

  idx.resize(widx);
  dat.resize(widx*dpn);

  f_close(stream);
}

template<class T, class S> inline
void read_indices_with_data(SF::vector<T> & idx, SF::vector<S> & dat,
                            const std::string filename,
                            const sf_mesh & mesh,
                            const SF::SF_nbr nbr,
                            const bool algebraic,
                            const int dpn,
                            MPI_Comm comm)
{
  const SF::vector<mesh_int_t> & dd_nbr = mesh.get_numbering(nbr);
  hashmap::unordered_map<mesh_int_t, mesh_int_t> dd_map;

  if(algebraic) {
    const SF::vector<mesh_int_t> & alg_nod = mesh.pl.algebraic_nodes();
    for(size_t i = 0; i<alg_nod.size(); i++) {
      mesh_int_t an = alg_nod[i];
      dd_map[dd_nbr[an]] = an;
    }

  } else {
    for(size_t i = 0; i<dd_nbr.size(); i++)
      dd_map[dd_nbr[i]] = i;
  }

  read_indices_with_data(idx, dat, filename, dd_map, dpn, comm);
}

template<class T, class S> inline
void read_indices_with_data(SF::vector<T> & idx, SF::vector<S> & dat,
                            const std::string filename,
                            const SF::vector<mesh_int_t> & dd_nbr,
                            const int dpn,
                            MPI_Comm comm)
{
  hashmap::unordered_map<mesh_int_t, mesh_int_t> dd_map;

  for(size_t i = 0; i<dd_nbr.size(); i++)
    dd_map[dd_nbr[i]] = i;

  read_indices_with_data(idx, dat, filename, dd_map, dpn, comm);
}

void warn_when_passing_intra_vtx(const std::string filename);

void indices_from_region_tag(SF::vector<mesh_int_t> & idx, const sf_mesh & mesh, const int tag);

void indices_from_geom_shape(SF::vector<mesh_int_t> & idx, const sf_mesh & mesh, const geom_shape shape, const bool nodal);

// compute volume of a node set using a mass matrix
SF_real get_volume_from_nodes(sf_mat & mass, SF::vector<mesh_int_t> & local_idx);


}  // namespace opencarp

#endif