fem_utils.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: Apache-2.0
 
#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_region_tags(SF::vector<mesh_int_t> & idx, const sf_mesh & mesh, const hashmap::unordered_set<int> & tags);
 
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