fem_utils.cc

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
// ----------------------------------------------------------------------------

#include "SF_init.h"
#include "fem.h"

namespace opencarp {

void parse_comment_line(char* buff, const int buffsize, std::map<std::string,std::string> & metadata)
{
  // we skip comment lines with no key value pair
  if(! has_char(buff, buffsize, '=')) return;

  // there might be multiple comment characters and spaces preceding a key-value pair
  remove_preceding_char(buff, buffsize, COMMENT_CHAR);
  remove_preceding_char(buff, buffsize, ' ');

  const int wsize = 1024;
  char lword[wsize], cword[wsize], rword[wsize];

  int numread = sscanf(buff, "%s %s %s", lword, cword, rword);

  if(numread == 3 && cword[0] == '=') {
    // input was perfect
    std::string key = lword, value = rword;
    metadata[key] = value;
  }
  else if(numread == 1) {
    // we have key=value without spaces
    int ridx=0, widx=0;
    while(ridx < wsize && lword[ridx] != '=') cword[widx++] = lword[ridx++];
    cword[widx] = '\0';

    std::string key = cword;

    widx=0, ridx += 1;
    while(ridx < wsize && lword[ridx] != '\0') cword[widx++] = lword[ridx++];
    cword[widx] = '\0';

    std::string value = cword;

    metadata[key] = value;
  }
  else {
    fprintf(stderr, "%s warning: Malformed key-value pair:\n%s\n", __func__, buff);
  }
}

void read_metadata(const std::string filename, std::map<std::string,std::string> & metadata, MPI_Comm comm)
{
  // There are 3 ways to add metadata to a vtx file:
  // 1) the new way:
  //    # key1 = value1
  //    # key2 = value2
  //    # key3 = value3
  // 2) nvtx # mesh [intra|extra] nelem npts
  // 3) number of nodes
  //    [intra|extra]
  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);

  int err = 0;
  FILE_SPEC stream = NULL;

  if(rank == 0) {
    stream = f_open(filename.c_str(), "r");
    if(stream == NULL) err++;
  }

  // return on file open error
  if(get_global(err, MPI_SUM)) {
    log_msg(0,5,0, "%s error: Could not read vtx data. Aborting!", __func__);
    EXIT(1);
  }

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

  // parse type 1 metadata lines. for simplicity we parse on each rank.
  char* ptr = f_gets_par(readbuff, str_size, stream, comm);
  remove_preceding_char(readbuff, str_size, ' ');

  while(ptr != NULL && readbuff[0] == COMMENT_CHAR) {
    parse_comment_line(readbuff, str_size, metadata);
    ptr = f_gets_par(readbuff, str_size, stream, comm);
    remove_preceding_char(readbuff, str_size, ' ');
  }

  // parse type 2 metadata
  if(has_char(readbuff, str_size, '#')) {
    int ibuff, ne, nn;
    char iestr[128], nestr[128], nnstr[128];
    int nread = sscanf(readbuff, "%d # %*s %31s %d %d", &ibuff, iestr, &ne, &nn);

    if(nread == 4) {
      metadata["grid"]     = iestr;
      metadata["nelem"]    = std::to_string(ne);
      metadata["nnode"]    = std::to_string(nn);
    }
  }

  // parse type 3 metadata
  ptr = f_gets_par(readbuff, str_size, stream, comm);

  if(ptr) {
    remove_preceding_char(readbuff, str_size, ' ');
    char grid[128];
    sscanf(readbuff, "%s", grid);

    if(strcmp(grid, "intra") == 0) {
      metadata["grid"] = "intra";
    } else if(strcmp(grid, "extra") == 0) {
      metadata["grid"] = "extra";
    }
  }

  f_close(stream);
}

char* skip_comments(FILE_SPEC stream, char* readbuff, size_t buffsize, MPI_Comm comm)
{
  char* ptr = NULL;
  do {
    ptr = f_gets_par(readbuff, buffsize, stream, comm);
    remove_preceding_char(readbuff, buffsize, ' ');
  }
  while(ptr != NULL && readbuff[0] == COMMENT_CHAR);

  return ptr;
}

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

  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
    if(mesh.tag[eidx] == tag) {
      for(mesh_int_t i = mesh.dsp[eidx]; i<mesh.dsp[eidx + 1]; i++)
        idx_set.insert(mesh.con[i]);
    }
  }

  idx.assign(idx_set.begin(), idx_set.end());
}

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

  if(nodal) {
    idx.reserve(mesh.l_numpts);

    for(size_t i=0; i<mesh.l_numpts; i++) {
      p.x = mesh.xyz[i*3 + 0], p.y = mesh.xyz[i*3 + 1], p.z = mesh.xyz[i*3 + 2];
      if(point_in_shape(p, shape))
        idx.push_back(i);
    }
  }
  else {
    idx.reserve(mesh.l_numelem);

    for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
      p.x = p.y = p.z = 0.0;

      for(mesh_int_t i = mesh.dsp[eidx]; i<mesh.dsp[eidx + 1]; i++) {
        mesh_int_t c = mesh.con[i];
        p.x += mesh.xyz[c*3 + 0];
        p.y += mesh.xyz[c*3 + 1];
        p.z += mesh.xyz[c*3 + 2];
      }
      p /= double(mesh.dsp[eidx + 1] - mesh.dsp[eidx]);

      if(point_in_shape(p, shape))
        idx.push_back(eidx);
    }
  }
}

SF_real get_volume_from_nodes(sf_mat & mass, SF::vector<mesh_int_t> & local_idx)
{
  assert(mass.mesh_ptr() != NULL);

  sf_vec* weights; SF::init_vector(&weights, *mass.mesh_ptr(), 1, sf_vec::algebraic);
  sf_vec* vols; SF::init_vector(&vols, *mass.mesh_ptr(), 1, sf_vec::algebraic);

  const SF::vector<mesh_int_t> & petsc_nod = mass.mesh_ptr()->get_numbering(SF::NBR_PETSC);

  // generate global petsc indices
  SF::vector<SF_int> petsc_idx(local_idx.size());
  for(size_t i=0; i<local_idx.size(); i++) petsc_idx[i] = petsc_nod[local_idx[i]];

  weights->set(0);
  weights->set(petsc_idx, 1.0);

  mass.mult(*weights, *vols);
  SF_real V = vols->sum();

  return V;
}

void warn_when_passing_intra_vtx(const std::string filename)
{
  std::map<std::string,std::string> metadata;
  read_metadata(filename, metadata, PETSC_COMM_WORLD);

  if(metadata.count("grid")) {
    if(metadata["grid"].compare("intra") == 0) {
      log_msg(0,3,0, "Warning: openCARP requires input vtx indices to be of from the input mesh, i.e. of type \"extra\".");
    }
  } else {
    log_msg(0,3,0, "%s warning: Could not derive grid type info from file %s", __func__, filename.c_str());
  }
}

}  // namespace opencarp