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"

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 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("%d ", 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)
{
  for(int i=0; i<npts; i++)
    printf(" ( %g %g %g ) \n", pts[i].x, pts[i].y, pts[i].z);

  printf("\n");
}

//#define DEBUG_INT

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);

  // if we have a region scale
  if(material.el_scale.size()) {
    eVals[0] *= material.el_scale[eidx];
    eVals[1] *= material.el_scale[eidx];
    eVals[2] *= material.el_scale[eidx];
  }

  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);

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

  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