electric_integrators.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 "electric_integrators.h"
#include "sf_interface.h"
#include "SF_init.h"
 
namespace opencarp {
 
#define ELEC_STIFFNESS 1
#define ELEC_MASS      2
 
inline int get_preferred_int_order(SF::elem_t & etype, int mat_type)
{
  switch(mat_type) {
    default:
    case ELEC_MASS: return 2;
 
    case ELEC_STIFFNESS: {
      switch(etype) {
        default:
        case SF::Quad:
        case SF::Hexa:
        case SF::Prism:
        case SF::Pyramid:
          return 2;
 
        case SF::Line:
        case SF::Tetra:
        case SF::Tri:
          return 1;
      }
    }
  }
}
 
void read_el_scale_vec(const char* file, mesh_t mt, SF::vector<double>& el_scale, int& el_scale_dpn)
{
  sf_mesh& mesh = get_mesh(mt);
  int rank; MPI_Comm_rank(mesh.comm, &rank);
 
  int dpn = 1;
  FILE* fd = NULL;
 
  if (rank == 0) {
    fd = fopen(file, "r");
    assert(fd != NULL);
 
    char buf[1024];
    fgets(buf, sizeof(buf), fd);
 
    // A bare integer (no decimal, no second token) of value 1 or 3 is the header.
    // Anything else is the first data row; infer dpn from its field count.
    int dpn_candidate; char rest[4];
    if (sscanf(buf, "%d %3s", &dpn_candidate, rest) == 1 && (dpn_candidate == 1 || dpn_candidate == 3)) {
      dpn = dpn_candidate;
    } else {
      float v[3];
      dpn = (sscanf(buf, "%f %f %f", v, v+1, v+2) == 3) ? 3 : 1;
      fseek(fd, 0, SEEK_SET);
    }
  }
 
  MPI_Bcast(&dpn, 1, MPI_INT, 0, mesh.comm);
 
  sf_vec* escale;
  SF::init_vector(&escale, mesh, dpn, sf_vec::elemwise);
  size_t nrd = escale->read_ascii(fd);
  if (rank == 0) fclose(fd);
 
  if (nrd != (size_t)(mesh.g_numelem * dpn)) {
    log_msg(0, 4, 0, "%s error: %s scale file has wrong size (expected %zu elements * dpn=%d); skipping",
            __func__, get_mesh_type_name(mt), mesh.g_numelem, dpn);
    delete escale;
    return;
  }
 
  if (get_size() > 1) {
    if (get_permutation(mt, ELEM_PETSC_TO_CANONICAL, dpn) == NULL)
      register_permutation(mt, ELEM_PETSC_TO_CANONICAL, dpn);
    SF::scattering& sc = *get_permutation(mt, ELEM_PETSC_TO_CANONICAL, dpn);
    sc(*escale, false);
  }
 
  SF_real* p = escale->ptr();
  el_scale.assign(p, p + escale->lsize());
  escale->release_ptr(p);
  delete escale;
 
  el_scale_dpn = dpn;
}
 
void print_element_info(const SF::element_view<mesh_int_t,mesh_real_t> & elem,
                        const int* reg,
                        const double vol,
                        const double* eVals,
                        const char* msg_start)
{
  if(msg_start)
    printf("\n%s", msg_start);
 
  if(reg)
    printf("gregion[%d], ", *reg);
 
  printf("element %ld: ", (long int) elem.global_element_index(SF::NBR_ELEM_REF));
 
  for(int i=0; i<elem.num_nodes(); i++)
    printf("%lld ", static_cast<long long>(elem.global_node(i, SF::NBR_REF)));
 
  printf("\nVolume: %g", vol);
 
  if(eVals) {
    SF::Point fib = elem.fiber(), she = elem.sheet();
    printf(", fib: (%g, %g, %g), she: (%g, %g, %g), G: (%g, %g, %g)",
           fib.x, fib.y, fib.z, she.x, she.y, she.z, eVals[0]*1e3, eVals[1]*1e3, eVals[2]*1e3);
  }
 
  printf("\n\n");
}
 
void get_conductivity_evals(const elecMaterial & emat,
                            const bool is_bath,
                            double* evals)
{
  if(is_bath == false) {
    switch(emat.g) {
      case intra_cond:
        evals[0] = emat.InVal[0], evals[1] = emat.InVal[1], evals[2] = emat.InVal[2];
        break;
 
      case extra_cond:
        evals[0] = emat.ExVal[0], evals[1] = emat.ExVal[1], evals[2] = emat.ExVal[2];
        break;
 
      case sum_cond:
        evals[0] = emat.InVal[0] + emat.ExVal[0],
        evals[1] = emat.InVal[1] + emat.ExVal[1],
        evals[2] = emat.InVal[2] + emat.ExVal[2];
        break;
 
      case para_cond:
        evals[0] = (emat.InVal[0] * emat.ExVal[0]) / (emat.InVal[0] + emat.ExVal[0]),
        evals[1] = (emat.InVal[1] * emat.ExVal[1]) / (emat.InVal[1] + emat.ExVal[1]),
        evals[2] = (emat.InVal[2] * emat.ExVal[2]) / (emat.InVal[2] + emat.ExVal[2]);
        break;
    }
  }
  else
    evals[0] = emat.BathVal[0], evals[1] = emat.BathVal[1], evals[2] = emat.BathVal[2];
}
 
void get_conductivity(const double* evals,
                      const SF::Point *f,
                      const SF::Point *s,
                      SF::dmat<double> & cond)
{
  cond.set_size(3,3);
 
  // longitudonal, transversal, and normal direction eigenvalues
  double gl = evals[0], gt = evals[1], gn = evals[2];
 
  if(s == nullptr) {
    // we compute the non-orthotropic conductivity:
    // sigma = gl f f^T + gt (I - f f^T)
    //       = gt I + (gl - gt) f f^T
    SF::outer_prod(*f, *f, gl - gt, cond[0]);
    cond[0][0] += gt, cond[1][1] += gt, cond[2][2] += gt;
  }
  else {
    // we compute the orthotropic conductivity:
    // sigma = gl f f^T + gt s s^T + gn n n^T
    SF::Point n = cross(*f, *s);
 
    SF::outer_prod(*f, *f, gl, cond[0], false);
    SF::outer_prod(*s, *s, gt, cond[0], true);
    SF::outer_prod( n,  n, gn, cond[0], true);
  }
}
 
void print_points(SF::Point* pts, int npts, int eidx)
{
  for(int i=0; i<npts; i++){
    if(pts[i].x!=pts[i].x or pts[i].y!=pts[i].y or pts[i].z!=pts[i].z){
      printf(" eidx = %d print_points ( %g %g %g ) \n", eidx, pts[i].x, pts[i].y, pts[i].z);
      printf("\n");
    }
  }
}
 
void elec_stiffness_integrator::operator() (const SF::element_view<mesh_int_t,mesh_real_t> & elem, SF::dmat<double> & buff)
{
  const size_t eidx      = elem.element_index();
  SF::elem_t   type      = elem.type();
  mesh_int_t   nnodes    = elem.num_nodes();
  const int    int_order = get_preferred_int_order(type, ELEC_STIFFNESS);
 
  const int    reg  = material.regionIDs.size() ? material.regionIDs[eidx] : 0;
  const RegionSpecs & reg_spec = material.regions[reg];
 
  assert(reg_spec.material->material_type == ElecMat);
  const elecMaterial & mat = *static_cast<elecMaterial*>(reg_spec.material);
  SF::Point fib = elem.fiber(), she = elem.sheet();
 
  bool is_bath   = SF::inner_prod(fib, fib) < 0.01; // bath elements have no fiber direction
  bool has_sheet = elem.has_sheet();            // do we have a sheet direciton
 
  double eVals[3];
  get_conductivity_evals(mat, is_bath, eVals);
 
  assert(material.el_scale_dpn == 0 || material.el_scale.size() > 0);
 
  if (material.el_scale_dpn == 1) {
    double s = material.el_scale[eidx];
    eVals[0] *= s; eVals[1] *= s; eVals[2] *= s;
  } else if (material.el_scale_dpn == 3) {
    if (is_bath) {
      double s = material.el_scale[3*eidx];
      eVals[0] *= s; eVals[1] *= s; eVals[2] *= s;
    } else {
      eVals[0] *= material.el_scale[3*eidx];
      eVals[1] *= material.el_scale[3*eidx + 1];
      eVals[2] *= material.el_scale[3*eidx + 2];
    }
  }
 
  if(is_bath) {
    fib = {1, 0, 0};
    she = {0, 1, 0};
  }
 
  SF::get_transformed_pts(elem, lpts, fib);
 
#ifdef DEBUG_INT
  print_points(lpts, nnodes);
#endif
 
  get_conductivity(eVals, &fib, has_sheet ? &she : nullptr, cond);
 
  int ndof = nnodes;
  shape.set_size(4, ndof), gshape.set_size(4, ndof);
  buff.assign(ndof, ndof, 0.0);
 
  int nint = 0;
  elem.integration_points(int_order, ipts, w, nint);
 
  double vol = 0.0;
  bool have_nan = false;
 
  // Integration loop
  for(int iidx=0; iidx < nint; iidx++)
  {
    double detj;
    // get shape function and its spatial derivatives
    SF::reference_shape(type, ipts[iidx], shape);
    SF::jacobian_matrix(shape, nnodes, lpts, J);
    SF::invert_jacobian_matrix(type, J, detj);
 
#ifdef DEBUG_INT
    SF::dmat<double> Jbuff(3,3);
    memcpy(Jbuff.data(), J, 9*sizeof(double));
    std::string j_name = "J_" + std::to_string(iidx);
    Jbuff.disp(j_name.c_str());
#endif
 
    SF::shape_deriv(J, shape, ndof, gshape);
 
    double dx = w[iidx] * detj;
    vol += dx;
 
    for (int k = 0; k < ndof; k++) {
      // derivatives of shape function k
      double skx = gshape[1][k], sky = gshape[2][k], skz = gshape[3][k];
      double cx  = cond[0][0] * skx + cond[0][1] * sky + cond[0][2] * skz;
      double cy  = cond[1][0] * skx + cond[1][1] * sky + cond[1][2] * skz;
      double cz  = cond[2][0] * skx + cond[2][1] * sky + cond[2][2] * skz;
 
      for (int l = 0; l < ndof; l++) {
        // derivatives of shape function l
        double slx = gshape[1][l], sly = gshape[2][l], slz = gshape[3][l];
        double val = - (cx*slx + cy * sly + cz * slz) * dx;
        buff[k][l] += val;
 
        if(val != val) have_nan = true;
      }
    }
  }
 
 
#ifdef DEBUG_INT
  printf("\n vol: %g \n\n", vol);
  buff.disp("Element Matrix:");
#endif
 
  if(have_nan) {
    print_element_info(elem, &reg, vol, eVals, "NaN detected: ");
    exit(EXIT_FAILURE);
  }
}
 
void elec_stiffness_integrator::dpn(mesh_int_t & row_dpn, mesh_int_t & col_dpn)
{
  row_dpn = 1;
  col_dpn = 1;
}
 
 
void mass_integrator::operator() (const SF::element_view<mesh_int_t,mesh_real_t> & elem, SF::dmat<double> & buff)
{
  SF::elem_t type      = elem.type();
  mesh_int_t nnodes    = elem.num_nodes();
  const int  int_order = get_preferred_int_order(type, ELEC_MASS);
 
  int ndof = nnodes;
  shape.set_size(4, ndof);
  buff.assign(ndof, ndof, 0.0);
 
  SF::Point fib = {1, 0, 0};  // this is a dummy fibre to satisfy get_transformed_pts()
 
  double vol = 0.0;
  bool have_nan = false;
  int nint;
 
  elem.integration_points(int_order, ipts, w, nint);
  SF::get_transformed_pts(elem, lpts, fib, phys_defined(PHYSREG_INTRA_ELEC));
 
  // Integration loop
  for(int iidx=0; iidx < nint; iidx++)
  {
    double detj;
    // get shape function and its spatial derivatives
    SF::reference_shape(type, ipts[iidx], shape);
    SF::jacobian_matrix(shape, nnodes, lpts, J);
    SF::invert_jacobian_matrix(type, J, detj);
 
    double dx = w[iidx] * detj;
    vol += dx;
 
    for (int k = 0; k < ndof; k++) {
      double sk = shape[0][k];
      for (int l = 0; l < ndof; l++) {
        double sl = shape[0][l];
        double val = sk * sl * dx;
        buff[k][l] += val;
 
        if(val != val) have_nan = true;
      }
    }
  }
 
#ifdef DEBUG_INT
  printf("\n vol: %g \n\n", vol);
  buff.disp("Mass Matrix:");
#endif
 
  if(have_nan) {
    print_element_info(elem, NULL, vol, NULL, "NaN detected: ");
    exit(EXIT_FAILURE);
  }
}
 
void mass_integrator::dpn(mesh_int_t & row_dpn, mesh_int_t & col_dpn)
{
  row_dpn = 1;
  col_dpn = 1;
}
 
}  // namespace opencarp