mechanics.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 "mechanics.h"
#include "mechanic_integrators.h"
#include "petsc_utils.h"
#include "constitutive_model.h"
#include "constitutive_model_library.h"
#include "SF_init.h"
 
namespace opencarp {
 
 
PetscErrorCode equilibrium_solver_internal_forces_petsc_helperfunction(SNES snes, Vec du, Vec resid, void * solver_pointer)
{
  Mechanics * mech = reinterpret_cast<Mechanics *>(user_globals::physics_reg[mech_phys]);
  equilibrium_solver * solver = reinterpret_cast<equilibrium_solver *>(solver_pointer);
  SF::petsc_vector<int, double> * position = reinterpret_cast<SF::petsc_vector<int, double> *>(solver->position);
  
  PetscErrorCode ierr;
  //change nodal coords based on du
  ierr = VecAssemblyBegin(du); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(du); CHKERRQ(ierr);
  ierr = VecAXPY(position->data, 1.0, du); CHKERRQ(ierr);
  ierr = VecAssemblyBegin(du); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(du); CHKERRQ(ierr);
  
  const SF::vector<mesh_int_t> & alg_nod = get_mesh(elasticity_msh).pl.algebraic_nodes();
  const SF::vector<mesh_int_t> & nbr_ref = get_mesh(elasticity_msh).get_numbering(SF::NBR_REF);
  const SF::vector<mesh_int_t> & nbr_petsc = get_mesh(elasticity_msh).get_numbering(SF::NBR_PETSC);
  
  solver->update_node_coords();
  solver->f_internal->set(0.0); // set f_internal to 0 before calculating new force
  solver->calc_internal_forces(mech->mtype, mech->logger);
  
  mech->surfmesh.xyz_Euler = get_mesh(elasticity_msh).xyz_Euler;
  
  solver->f_external->set(0.0);
  
  for (int nmbc_idx = 0; nmbc_idx < param_globals::num_mech_nbcs; nmbc_idx++) {
    solver->apply_external_forces(mech->mtype, mech->logger, mech->surfmesh, solver->nmbcs[nmbc_idx]);
  }
  SF::petsc_vector<int, double> * f_internal = reinterpret_cast<SF::petsc_vector<int, double> *>(solver->f_internal);
  SF::petsc_vector<int, double> * residuum = reinterpret_cast<SF::petsc_vector<int, double> *>(solver->residuum);
  SF::petsc_vector<int, double> * f_external = reinterpret_cast<SF::petsc_vector<int, double> *>(solver->f_external);
  residuum->set(0.0);
  residuum->add_scaled(*f_internal, 1.0);
  residuum->add_scaled(*f_external, -1.0);
  ierr = VecAssemblyBegin(resid); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(resid); CHKERRQ(ierr);
  ierr = VecCopy(residuum->data, resid); CHKERRQ(ierr);
  ierr = VecAssemblyBegin(resid); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(resid); CHKERRQ(ierr);
  
  // reset nodal coords to the ones before this iteration
  ierr = VecAssemblyBegin(du); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(du); CHKERRQ(ierr);
  ierr = VecAXPY(position->data, -1.0, du); CHKERRQ(ierr);
  ierr = VecAssemblyBegin(du); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(du); CHKERRQ(ierr);
  ierr = VecAssemblyBegin(resid); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(resid); CHKERRQ(ierr);
  solver->update_node_coords();
#ifdef DEBUG_INT
  PetscScalar resid_norm;
  VecNorm(resid, NORM_2, &resid_norm);
  printf("\n Residuum on rank %d: %e", PetscGlobalRank, resid_norm);
  ierr = VecAssemblyBegin(resid); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(resid); CHKERRQ(ierr);
#endif
  return ierr;
}
 
PetscErrorCode equilibrium_solver_jacobian_petsc_helperfunction(SNES snes, Vec du, Mat jacobian, Mat preconditioner_mat, void * solver_pointer)
{
  PetscErrorCode ierr;
  Mechanics * mech = reinterpret_cast<Mechanics *>(user_globals::physics_reg[mech_phys]);
  equilibrium_solver * solver = reinterpret_cast<equilibrium_solver *>(solver_pointer);
  SF::petsc_vector<int, double> * position = reinterpret_cast<SF::petsc_vector<int, double> *>(solver->position);
  ierr = VecAssemblyBegin(du); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(du); CHKERRQ(ierr);
  ierr = VecAXPY(position->data, 1.0, du); CHKERRQ(ierr);
  ierr = VecAssemblyBegin(du); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(du); CHKERRQ(ierr);
  solver->K_eq->zero();
  solver->update_node_coords();
  solver->rebuild_stiffness(mech->mtype, mech->logger);
  SF::petsc_matrix<mesh_int_t, mesh_real_t> * K_eq = reinterpret_cast<SF::petsc_matrix<mesh_int_t, mesh_real_t> *>(solver->K_eq);
  ierr = MatAssemblyBegin(preconditioner_mat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatAssemblyEnd(preconditioner_mat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatCopy(K_eq->data, preconditioner_mat, DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);
  ierr = MatAssemblyBegin(preconditioner_mat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatAssemblyEnd(preconditioner_mat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  
  // set position again to where it was before calculating the jacobian
  ierr = VecAssemblyBegin(du); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(du); CHKERRQ(ierr);
  ierr = VecAXPY(position->data, -1.0, du); CHKERRQ(ierr);
  ierr = VecAssemblyBegin(du); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(du); CHKERRQ(ierr);
  
  solver->update_node_coords();
  
  
  if(jacobian != preconditioner_mat)
  {
    ierr = MatAssemblyBegin(jacobian, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(jacobian, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  }
  return ierr;
}
 
void Mechanics::initialize()
{
  double t1, t2;
  get_time(t1);
    
  set_dir(OUTPUT);
    
  // open logger
  logger = f_open("mechanics.log", param_globals::experiment != 4 ? "w" : "r");
  
  setup_mappings();
    
  double global_time = user_globals::tm_manager->time;
  timer_idx = user_globals::tm_manager->add_eq_timer(global_time, param_globals::tend, 0,
                                                       param_globals::dt, 0, "elec::ref_dt", "TS");
    
  // generate the surface mesh
  sf_mesh mesh = get_mesh(elasticity_msh);
  mesh.xyz_Euler = mesh.xyz;
  compute_surface_mesh(mesh, SF::NBR_PETSC, surfmesh);
  // give surfmesh the same nodes, numberings etc
  surfmesh.g_numpts = mesh.g_numpts;
  surfmesh.l_numpts = mesh.l_numpts;
  surfmesh.xyz = mesh.xyz;
  surfmesh.xyz_Euler = mesh.xyz_Euler;
  surfmesh.nbr = mesh.nbr;
  surfmesh.pl = mesh.pl;
  // register the mechanics material model
  register_material_model();
  // set up the linear equation systems.
  log_msg(logger,0,0, "\n Setting up Mechanical Solvers \n");
  setup_solvers();
    
  // prepare the mechanics output. we skip it if we do post-processing
  if(param_globals::experiment != EXP_POSTPROCESS)
      setup_output();
  
  this->initialize_time += timing(t2, t1);
}
 
void Mechanics::setup_mappings()  //currently mostly copied from electrics setup_mapping
{
  bool mech_mesh_exits = mesh_is_registered(elasticity_msh);
  assert(mech_mesh_exits);
  const int dpn = 3;
  
  // It may be that another physic (e.g. electrics) has already computed the mappings,
  // thus we first test their existence
  if(get_scattering(elasticity_msh, ALG_TO_NODAL, dpn) == NULL) {
    log_msg(logger, 0, 0, "%s: Setting up mechanics algebraic-to-nodal scattering.", __func__);
    register_scattering(elasticity_msh, ALG_TO_NODAL, dpn);
  }
  if(get_permutation(elasticity_msh, PETSC_TO_CANONICAL, dpn) == NULL) {
    log_msg(logger, 0, 0, "%s: Setting up mechanics PETSc to canonical permutation.", __func__);
    register_permutation(elasticity_msh, PETSC_TO_CANONICAL, dpn);
  }
}
 
void Mechanics::compute_step()
{
  double t1, t2;
  get_time(t1);
  sf_mesh & mesh = get_mesh(elasticity_msh);
//  checkpointing(); // not yet implemented
 
  if(&eq_solver)
  {
    size_t nelem = mesh.g_numelem;
    surfmesh.xyz_Euler = mesh.xyz_Euler; //updates nodal coordinates of surface mesh
    eq_solver.rebuild_matrices(mtype, logger);
    
    const SF::vector<mesh_int_t>& alg_nod_surf  = surfmesh.pl.algebraic_nodes();
    const SF::vector<mesh_int_t>& alg_dsp_surf  = surfmesh.pl.algebraic_layout();
    const SF::vector<mesh_int_t>& ref_nbr_surf  = surfmesh.get_numbering(SF::NBR_REF);
    const SF::vector<mesh_int_t>& alg_nod_vol  = mesh.pl.algebraic_nodes();
    const SF::vector<mesh_int_t>& alg_dsp_vol  = mesh.pl.algebraic_layout();
    const SF::vector<mesh_int_t>& ref_nbr_vol  = mesh.get_numbering(SF::NBR_REF);
    for (int nmbc_idx = 0; nmbc_idx < param_globals::num_mech_nbcs; nmbc_idx++) {
      eq_solver.apply_external_forces(mtype, logger, surfmesh, eq_solver.nmbcs[nmbc_idx]);
    }
    eq_solver.solve();
    // update nodal coordinates with calculated du
    eq_solver.position->add_scaled(*eq_solver.du, 1.0);
    eq_solver.displacement->add_scaled(*eq_solver.du, 1.0);
    eq_solver.initial_guess->set(0.0);
    eq_solver.initial_guess->add_scaled(*eq_solver.du, 1.0);
    eq_solver.update_node_coords();
  }
  this->compute_time += timing(t2, t1);
}
 
void Mechanics::output_step()
{
    double t1, t2;
    get_time(t1);
    output_manager.write_data();
    this->output_time += timing(t2, t1);
}
 
void Mechanics::destroy()
{
    // TODO
    // close logger
    f_close(logger);
}
 
void Mechanics::setup_output()
{
  set_dir(OUTPUT);
  output_manager.register_output(eq_solver.displacement, elasticity_msh, 3, param_globals::ufile, "mm");
  output_manager.register_output(eq_solver.position, elasticity_msh, 3, param_globals::dynptfile, "mm");
  output_manager.register_output(eq_solver.f_internal, elasticity_msh, 3, param_globals::Ffile, "kN");
  output_manager.register_output(eq_solver.residuum, elasticity_msh, 3, param_globals::Residfile, "kN");
}
 
void Mechanics::dump_matrices()
{
    std::string bsname = param_globals::dump_basename;
    std::string fn;
    // TODO
    set_dir(OUTPUT);
}
 
 
double Mechanics::timer_val(const int timer_id)
{
    double val = 0.0;
    // TODO
    return val;
}
 
std::string Mechanics::timer_unit(const int timer_id)
{
    std::string s_unit;
    // TODO
    
    return s_unit;
}
 
void Mechanics::setup_solvers()
{
  set_dir(OUTPUT);
  eq_solver.init();
  eq_solver.rebuild_matrices(mtype, logger);
}
 
void Mechanics::checkpointing()
{
    // TODO
  const timer_manager & tm = *user_globals::tm_manager;
}
 
void Mechanics::register_material_model()
{
  mtype.regions.resize(param_globals::num_mech_mat_regions);
  RegionSpecs* reg = mtype.regions.data();
  for (std::size_t mreg = 0; mreg < param_globals::num_mech_mat_regions; mreg++) {
    mechMaterial* mMat = new mechMaterial();
    mMat->material_type = MechMat;
    mMat->set_constitutive_model(param_globals::mech_mat_region[mreg].constitutive_model);
    reg[mreg].material = mMat;
  }
}
 
void equilibrium_solver::init()
{
  FILE_SPEC logger = get_physics(mech_phys)->logger;
  log_msg(logger, 0, 0,"Equilibrium solver initialized.");
  
  SF::init_solver(&nonlin_solver);
  
  sf_mesh & mechanics_mesh = get_mesh(elasticity_msh);
  mechanics_mesh.xyz_Euler = mechanics_mesh.xyz;
  sf_vec::ltype alg_type = sf_vec::algebraic;
  sf_vec::ltype ele_type = sf_vec::elemwise;
  sf_vec::ltype nodal_type = sf_vec::nodal;
  sf_vec::ltype ip_type = sf_vec::ipwise;
  const int dpn = 3; //dof per node
  
  SF::init_vector(&position, mechanics_mesh, dpn, alg_type);
  SF::init_vector(&reference_pos, mechanics_mesh, dpn, alg_type);
  SF::init_vector(&displacement, mechanics_mesh, dpn, alg_type);
  SF::init_vector(&du, displacement);
  SF::init_vector(&initial_guess, du);
  SF::init_vector(&f_internal, mechanics_mesh, dpn, alg_type);
  SF::init_vector(&f_external, mechanics_mesh, dpn, alg_type);
  SF::init_vector(&residuum, mechanics_mesh, dpn, alg_type);
  SF::init_vector(&dirichlet_bcs, mechanics_mesh, dpn, alg_type);
  SF::init_vector(&dirichlet_bcs_1, mechanics_mesh, dpn, alg_type);
  
  SF::element_view<mesh_int_t, mesh_real_t> elem(mechanics_mesh, SF::NBR_ELEM_REF);
  
  int max_row_entries = max_nodal_edgecount(mechanics_mesh);
  
  SF::init_matrix(&K_eq);
  K_eq->init(mechanics_mesh, dpn, dpn, max_row_entries);
  
  displacement->set(0.0);
  
  const SF::vector<mesh_int_t>& alg_nod  = mechanics_mesh.pl.algebraic_nodes();
  const SF::vector<mesh_int_t>& alg_dsp  = mechanics_mesh.pl.algebraic_layout();
  const SF::vector<mesh_int_t> & nbr = mechanics_mesh.get_numbering(SF::NBR_REF);
  std::size_t loc_node;
  
  double * pos = position->ptr();
  double * ref_pos = reference_pos->ptr();
  for(std::size_t i = 0; i < alg_nod.size(); i++){
    loc_node = alg_nod[i];
    for(std::size_t k = 0; k < dpn; k++){
      pos[3 * i + k] = mechanics_mesh.xyz_Euler[loc_node * 3 + k];
      ref_pos[3 * i + k] = mechanics_mesh.xyz[loc_node * 3 + k];
    }
  }
  position->release_ptr(pos);
  reference_pos->release_ptr(ref_pos);
  
  std::size_t num_mnbcs = param_globals::num_mech_nbcs;
  std::size_t num_mdbcs = param_globals::num_mech_dbcs;
  
  
  for(std::size_t dbc_idx = 0; dbc_idx < num_mdbcs; dbc_idx++)
  {
    dmbcs[dbc_idx].definition = mech_bcs::Dirichlet;
    dmbcs[dbc_idx].setup(dbc_idx);
  }
  
  
  for(std::size_t nbc_idx = 0; nbc_idx < num_mnbcs; nbc_idx++)
  {
    nmbcs[nbc_idx].definition = mech_bcs::Neumann;
    nmbcs[nbc_idx].setup(nbc_idx);
  }
  
  dirichlet_bcs->set(1); //set all positions of dirichlet_boundary_conditions to 1 as default.
  dirichlet_bcs_1->set(0);
  
  double * dbc_ptr = dirichlet_bcs->ptr();
  double * dbc1_ptr = dirichlet_bcs_1->ptr();
 
  for(size_t dbc_idx = 0; dbc_idx < num_mdbcs; dbc_idx++)
  {
    switch (dmbcs[dbc_idx].dbc_type) {
      case mech_bcs::dbc_t::all:
        for (std::size_t lidx = 0; lidx < alg_nod.size(); lidx++) {
          std::size_t gidx = nbr[alg_nod[lidx]];
          bool found = FALSE;
          std::size_t vtx = 0;
          while (!found && (vtx < dmbcs[dbc_idx].vertices.size()))
          {
            found = (gidx == dmbcs[dbc_idx].vertices[vtx]);
            vtx++;
          }
          if (found) {
            for(size_t k = 0; k < dpn; k++)
            {
              dbc_ptr[3 * lidx + k] = 0;
              dbc1_ptr[3 * lidx + k] = 1;
            }
          }
        }
        break;
      case mech_bcs::dbc_t::x:
        for (std::size_t lidx = 0; lidx < alg_nod.size(); lidx++) {
          std::size_t gidx = nbr[alg_nod[lidx]];
          bool found = FALSE;
          std::size_t vtx = 0;
          while (!found && (vtx < dmbcs[dbc_idx].vertices.size()))
          {
            found = (gidx == dmbcs[dbc_idx].vertices[vtx]);
            vtx++;
          }
          if (found)
          {
            dbc_ptr[3 * lidx] = 0;
            dbc1_ptr[3 * lidx] = 1;
          }
        }
        break;
      case mech_bcs::dbc_t::y:
        for (std::size_t lidx = 0; lidx < alg_nod.size(); lidx++) {
          std::size_t gidx = nbr[alg_nod[lidx]];
          bool found = FALSE;
          std::size_t vtx = 0;
          while (!found && (vtx < dmbcs[dbc_idx].vertices.size()))
          {
            found = (gidx == dmbcs[dbc_idx].vertices[vtx]);
            vtx++;
          }
          if (found)
          {
            dbc_ptr[3 * lidx + 1] = 0;
            dbc1_ptr[3 * lidx + 1] = 1;
          }
        }
        break;
      case mech_bcs::dbc_t::z:
        for (std::size_t lidx = 0; lidx < alg_nod.size(); lidx++) {
          std::size_t gidx = nbr[alg_nod[lidx]];
          bool found = FALSE;
          std::size_t vtx = 0;
          while (!found && (vtx < dmbcs[dbc_idx].vertices.size()))
          {
            found = (gidx == dmbcs[dbc_idx].vertices[vtx]);
            vtx++;
          }
          if (found)
          {
            dbc_ptr[3 * lidx + 2] = 0;
            dbc1_ptr[3 * lidx + 2] = 1;
          }
        }
        break;
      case mech_bcs::dbc_t::xy:
        for (std::size_t lidx = 0; lidx < alg_nod.size(); lidx++) {
          std::size_t gidx = nbr[alg_nod[lidx]];
          bool found = FALSE;
          std::size_t vtx = 0;
          while (!found && (vtx < dmbcs[dbc_idx].vertices.size()))
          {
            found = (gidx == dmbcs[dbc_idx].vertices[vtx]);
            vtx++;
          }
          if (found)
          {
            dbc_ptr[3 * lidx] = 0;
            dbc_ptr[3 * lidx + 1] = 0;
            dbc1_ptr[3 * lidx] = 1;
            dbc1_ptr[3 * lidx + 1] = 1;
          }
        }
        break;
      case mech_bcs::dbc_t::xz:
        for (std::size_t lidx = 0; lidx < alg_nod.size(); lidx++) {
          std::size_t gidx = nbr[alg_nod[lidx]];
          bool found = FALSE;
          std::size_t vtx = 0;
          while (!found && (vtx < dmbcs[dbc_idx].vertices.size()))
          {
            found = (gidx == dmbcs[dbc_idx].vertices[vtx]);
            vtx++;
          }
          if (found)
          {
            dbc_ptr[3 * lidx] = 0;
            dbc_ptr[3 * lidx + 2] = 0;
            dbc1_ptr[3 * lidx] = 1;
            dbc1_ptr[3 * lidx + 2] = 1;
          }
        }
        break;
      case mech_bcs::dbc_t::yz:
        for (std::size_t lidx = 0; lidx < alg_nod.size(); lidx++) {
          std::size_t gidx = nbr[alg_nod[lidx]];
          bool found = FALSE;
          std::size_t vtx = 0;
          while (!found && (vtx < dmbcs[dbc_idx].vertices.size()))
          {
            found = (gidx == dmbcs[dbc_idx].vertices[vtx]);
            vtx++;
          }
          if (found)
          {
            dbc_ptr[3 * lidx + 1] = 0;
            dbc_ptr[3 * lidx + 2] = 0;
            dbc1_ptr[3 * lidx + 1] = 1;
            dbc1_ptr[3 * lidx + 2] = 1;
          }
        }
        break;
 
      default:
        log_msg(logger, 0, 0, "No direction type defined. Default option to fix displacement in all directions at given nodes");
        break;
    }
  }
  
  dirichlet_bcs->release_ptr(dbc_ptr);
  dirichlet_bcs_1->release_ptr(dbc1_ptr);
  
  setup_nonlinear_solver(logger);
}
 
void equilibrium_solver::rebuild_matrices(MaterialType mtype, FILE_SPEC logger)
{
  double t0, t1, dur;
  get_time(t0);
 
  rebuild_stiffness(mtype, logger);
}
 
void equilibrium_solver::rebuild_stiffness(MaterialType mtype, FILE_SPEC logger)
{
  double t0, t1, dur;
  int log_flag = param_globals::output_level > 1 ? ECHO : 0;
    
  sf_mesh & mesh     = get_mesh(elasticity_msh);
    
  get_time(t0);
    
  mech_stiffness_integrator m_stfn_integ(mtype);
  K_eq->zero();
    
  SF::assemble_matrix(*K_eq, mesh, m_stfn_integ);
  K_eq->finish_assembly();
    
  apply_bcs_jacobian();
  K_eq->finish_assembly();
    
  dur = timing(t1, t0);
}
 
void equilibrium_solver::calc_internal_forces(MaterialType mt, FILE_SPEC logger)
{
  double to, t1;
  mech_internal_forces_integrator int_force_int(mt);
    
  // get mesh reference
  sf_mesh & mesh = get_mesh(elasticity_msh);
  SF::assemble_vector(*f_internal, mesh, int_force_int);
  apply_bcs_forces();
}
 
void equilibrium_solver::apply_external_forces(MaterialType mt, FILE_SPEC logger, sf_mesh surfmesh, mech_bcs nmbcs)
{
  const SF::vector<mesh_int_t>& alg_nod  = surfmesh.pl.algebraic_nodes();
  const SF::vector<mesh_int_t>& alg_dsp  = surfmesh.pl.algebraic_layout();
  const SF::vector<mesh_int_t>& ref_nbr  = surfmesh.get_numbering(SF::NBR_REF);
 
  timer_manager & tm = *user_globals::tm_manager;
  if(nmbcs.nbc_type == mech_bcs::pressure)
  {
    log_msg(logger, 0, 0, "Applying external forces based on Neumann boundary conditions using a constant pressure for all elements");
    double time = tm.time + param_globals::dt; //Test
    double pressure = nmbcs.scaling[1] * time / param_globals::tend; //calculate pressure as function of the time
    f_external->set(0.0); // then set f_external to zero before each new calculation
    SF::element_view<mesh_int_t, mesh_real_t> view(surfmesh, SF::NBR_PETSC);
    mesh_real_t area;
    SF::elem_t type;
    bool bc0, bc1, bc2, bc3;
    SF::Point normal;
    SF::dmat<double> shape;
    SF::Point lpts[EP_MAX_LPOINTS];
    SF::Point ipts[EP_MAX_IPOINTS];
    double w[EP_MAX_IPOINTS];
    double F[9];
    double iF[9];
    double detF;
    SF::dmat<double> gshape;
    SF::vector<mesh_int_t> nbc_idx_test;
    double forces[SF_MAX_ELEM_NODES*3];
    double * f_ext = f_external->ptr();
    for(size_t eidx = 0; eidx < surfmesh.l_numelem; eidx++) //go through local elements
    {
      view.set_elem(eidx);
      type = view.type();
      shape.set_size(4, view.num_nodes());
      SF::Point fib = {1, 0, 0}; // just to satisfy get_transformed_pts
      int nint;
      int int_order;
      switch (type) {
        case SF::Tri:
          bc0 = std::find(std::begin(nmbcs.vertices), std::end(nmbcs.vertices), ref_nbr[view.node(0)]) != std::end(nmbcs.vertices);
          bc1 = std::find(std::begin(nmbcs.vertices), std::end(nmbcs.vertices), ref_nbr[view.node(1)]) != std::end(nmbcs.vertices);
          bc2 = std::find(std::begin(nmbcs.vertices), std::end(nmbcs.vertices), ref_nbr[view.node(2)]) != std::end(nmbcs.vertices);
          if(bc0 && bc1 && bc2) //check if all three nodes of the triangle surface element are part of the boundary condition vertices
          {
            int_order = 1; //hard-coded for linear elements
            view.integration_points(int_order, ipts, w, nint);
            SF::get_transformed_pts(view, lpts, fib);
            for(std::size_t iidx =0; iidx < nint; iidx++)
            {
              SF::reference_shape(type, ipts[iidx], shape);
//            Calculate area of surface element
              area = mag(cross(view.coord_upd(2) - view.coord_upd(0), view.coord_upd(1) - view.coord_upd(0)));
//            Calculate normal of surface element
              normal = normalize(cross(view.coord_upd(1) - view.coord_upd(0), view.coord_upd(2) - view.coord_upd(0)));
              for(std::size_t a = 0; a < view.num_nodes(); a++)
              {
                std::size_t walker = 0;
                bool alg_nod_dx_found = FALSE;
                while(walker < alg_nod.size() && !alg_nod_dx_found)
                {
                  alg_nod_dx_found = (alg_nod[walker] == view.node(a));
                  walker++;
                }
                double sa = shape[0][a];
                if(alg_nod_dx_found)
                {
                  f_ext[3 * (walker - 1)]     += sa * w[iidx] * area * pressure * normal.x;
                  f_ext[3 * (walker - 1) + 1] += sa * w[iidx] * area * pressure * normal.y;
                  f_ext[3 * (walker - 1) + 2] += sa * w[iidx] * area * pressure * normal.z;
                }
              }
            }
          }
          
          break;
        case SF::Quad:
          // TODO, current implementation only works for rectangular quadriliterals
          bc0 = std::find(std::begin(nmbcs.vertices), std::end(nmbcs.vertices), ref_nbr[view.node(0)]) != std::end(nmbcs.vertices);
          bc1 = std::find(std::begin(nmbcs.vertices), std::end(nmbcs.vertices), ref_nbr[view.node(1)]) != std::end(nmbcs.vertices);
          bc2 = std::find(std::begin(nmbcs.vertices), std::end(nmbcs.vertices), ref_nbr[view.node(2)]) != std::end(nmbcs.vertices);
          bc3 = std::find(std::begin(nmbcs.vertices), std::end(nmbcs.vertices), ref_nbr[view.node(3)]) != std::end(nmbcs.vertices);
          if(bc0 && bc1 && bc2 && bc3) //check if all four nodes of the rectangle surface element are part of the boundary condition vertices
          {
            int_order = 2; //hard-coded for linear elements
            view.integration_points(int_order, ipts, w, nint);
            SF::get_transformed_pts(view, lpts, fib);
            for(std::size_t iidx =0; iidx < nint; iidx++)
            {
              SF::reference_shape(type, ipts[iidx], shape);
//            Calculate area of surface element and scale with 0.25 for Gauss integration 
              area = 0.125 * (mag(cross(view.coord_upd(2) - view.coord_upd(0), view.coord_upd(1) - view.coord_upd(0))) + mag(cross(view.coord_upd(2) - view.coord_upd(0), view.coord_upd(3) - view.coord_upd(0))));
//            Calculate normal of surface element
              normal = normalize(cross(view.coord_upd(1) - view.coord_upd(0), view.coord_upd(2) - view.coord_upd(0)));
              for(std::size_t a = 0; a < view.num_nodes(); a++)
              {
                std::size_t walker = 0;
                bool alg_nod_dx_found = FALSE;
                while(walker < alg_nod.size() && !alg_nod_dx_found)
                {
                  alg_nod_dx_found = (alg_nod[walker] == view.node(a));
                  walker++;
                }
                double sa = shape[0][a];
                if(alg_nod_dx_found)
                {
                  f_ext[3 * (walker - 1)]     += sa * w[iidx] * area * pressure * normal.x;
                  f_ext[3 * (walker - 1) + 1] += sa * w[iidx] * area * pressure * normal.y;
                  f_ext[3 * (walker - 1) + 2] += sa * w[iidx] * area * pressure * normal.z;
                }
              }
            }
          }
          break;
        default:
          break;
      }
    }
    f_external->release_ptr(f_ext);
  }
  else
  {
    sf_vec * f_buff;
    f_buff->init(*f_external);
    f_buff->set(nmbcs.vertices, nmbcs.scaling);
  }
  apply_bcs_forces();
}
 
 
void equilibrium_solver::solve()
{
  double t0,t1;
  get_time(t0);
  du->set(0.0);
  (* nonlin_solver)(* du, * residuum);
}
 
void equilibrium_solver::setup_nonlinear_solver(FILE_SPEC logger)
{
  tol = param_globals::cg_tol_eq;
  max_it = param_globals::cg_maxit_eq;
  std::string solver_file;
  solver_file = param_globals::eq_options_file;
 
  equilibrium_solver* eq_solve = this;
    
  bool has_nullspace = true;
  nonlin_solver->setup_solver(*residuum, *K_eq, tol, max_it, param_globals::cg_norm_eq, "equilibrium PDE", has_nullspace, logger, solver_file.c_str(), "", equilibrium_solver_internal_forces_petsc_helperfunction, equilibrium_solver_jacobian_petsc_helperfunction, (void *) this);
}
 
void mech_bcs::setup(int idx)
{
  if(definition == Neumann)
  {
    const mech_nbcs & cur_mnbc = param_globals::mech_nbc[idx];
    nbc_type = (cur_mnbc.nbc_type == 0) ? mech_bcs::pressure : mech_bcs::force;
    direction.x = cur_mnbc.direction[0];
    direction.y = cur_mnbc.direction[1];
    direction.z = cur_mnbc.direction[2];
    // get a logger from the physics
    FILE_SPEC logger = get_physics(mech_phys)->logger;
    
    // these are the criteria for choosing how we extract the nodes for the Nbcs
    bool vertex_file_given = strlen(cur_mnbc.vtx_file) > 0;
    bool tag_index_given   = cur_mnbc.geomID > -1;
    // the mesh we need for computing the local vertex indices.
    const sf_mesh & mesh   = get_mesh(elasticity_msh);
    
    if(vertex_file_given && tag_index_given)
      log_msg(0,3,0, "%s warning: More than one region definitions set in Neumann boundary conditions %d", __func__, idx);
    
    if(vertex_file_given) {
      input_filename = cur_mnbc.vtx_file;
      
      log_msg(logger, 0, 0, "Neumann bcs region %d: Selecting vertices from file %s", idx, input_filename.c_str());
      
      set_dir(INPUT);
      
      // we read the indices. they are being localized w.r.t. the provided numbering. In our
      // case this is always the reference numbering.
      if(cur_mnbc.vtx_fcn)
      {
        if(nbc_type == 0)
        {
          log_msg(logger, 0, 0, "Trying to read a vertex file with scaled Neumann boundary conditions but specified pressure type! Aborting!");
          EXIT(1);
        }
      read_indices_with_data(vertices, scaling, input_filename, mesh, SF::NBR_REF, true, 1, PETSC_COMM_WORLD);
      }
      else
      {
        read_indices(vertices, input_filename, mesh, SF::NBR_REF, true, PETSC_COMM_WORLD);
        scaling.assign(vertices.size(), cur_mnbc.strength);
      }
      int gnum_idx = get_global(vertices.size(), MPI_SUM);
      if(gnum_idx == 0) {
        log_msg(0, 5, 0, "Neumann bcs region %d: Specified vertices are not in domain! Aborting!", idx);
        EXIT(1);
      }
    }
    else if(tag_index_given) {
      
      int tag = cur_mnbc.geomID;
      log_msg(logger, 0, 0, "Neumann bc region %d: Selecting vertices from tag %d", idx, tag);
      
      indices_from_region_tag(vertices, mesh, tag);
      // we restrict the indices to the algebraic subset
      SF::restrict_to_set(vertices, mesh.pl.algebraic_nodes());
    }
    else {
      log_msg(logger, 0, 0, "Neumann bc region %d: Selecting vertices from shape.", idx);
      
      geom_shape shape;
      shape.type   = geom_shape::shape_t(cur_mnbc.geom_type);
      shape.p0     = cur_mnbc.p0;
      shape.p1     = cur_mnbc.p1;
      shape.radius = cur_mnbc.radius;
      
      bool nodal = true;
      indices_from_geom_shape(vertices, mesh, shape, nodal);
      // we restrict the indices to the algebraic subset
      SF::restrict_to_set(vertices, mesh.pl.algebraic_nodes());
      
      int gsize = vertices.size();
      if(get_global(gsize, MPI_SUM) == 0) {
        log_msg(0,5,0, "error: Empty Neumann bc region [%d] def! Aborting!", idx);
        EXIT(1);
      }
      scaling.assign(gsize, cur_mnbc.strength); //read in the strength
    }
    if(cur_mnbc.dump_vtx_file) {
      SF::vector<mesh_int_t> glob_idx(vertices);
      mesh.pl.globalize(glob_idx);
 
      SF::vector<mesh_int_t> srt_idx;
      SF::sort_parallel(mesh.comm, glob_idx, srt_idx);
 
      size_t num_vtx = get_global(srt_idx.size(), MPI_SUM);
      int rank = get_rank();
 
      FILE_SPEC f = NULL;
      int err = 0;
 
      if(rank == 0) {
        char dump_name[1024];
        sprintf(dump_name, "NBCs_%d.vtx", idx);
 
        f = f_open(dump_name, "w");
        if(!f) err++;
        else {
          fprintf(f->fd, "%zd\nextra\n", num_vtx);
        }
      }
 
      if(!get_global(err, MPI_SUM)) {
        print_vector(mesh.comm, srt_idx, 1, f ? f->fd : NULL);
      } else {
        log_msg(0, 4, 0, "error: Neumann bc region [%d] cannot be dumped!");
      }
 
      // only root really does that
      f_close(f);
    }
  }
 
  if(definition == Dirichlet)
  {
    const mech_dbcs & cur_mdbc = param_globals::mech_dbc[idx];
    
    // assigning dbc_type
    switch (cur_mdbc.dbc_type) {
      case 0:
        dbc_type = all;
        break;
      case 1:
        dbc_type = x;
        break;
      case 2:
        dbc_type = y;
        break;
      case 3:
        dbc_type = z;
        break;
      case 4:
        dbc_type = xy;
        break;
      case 5:
        dbc_type = xz;
        break;
      case 6:
        dbc_type = yz;
        break;
      default:
        dbc_type = all;
        break;
    }
    
    // get a logger from the physics
    FILE_SPEC logger = get_physics(mech_phys)->logger;
    
    // these are the criteria for choosing how we extract the Dbc nodes
    bool vertex_file_given = strlen(cur_mdbc.vtx_file) > 0;
    bool tag_index_given   = cur_mdbc.geomID > -1;
    // the mesh we need for computing the local vertex indices.
    const sf_mesh & mesh   = get_mesh(elasticity_msh);
    
    if(vertex_file_given && tag_index_given)
      log_msg(0,3,0, "%s warning: More than one region definitions set in Dirichlet boundary conditions %d", __func__, idx);
    
    if(vertex_file_given) {
      input_filename = cur_mdbc.vtx_file;
      
      log_msg(logger, 0, 0, "Dirichlet bcs region %d: Selecting vertices from file %s", idx, input_filename.c_str());
      
      set_dir(INPUT);
      
      // we read the indices. they are being localized w.r.t. the provided numbering. In our
      // case this is always the reference numbering.
      if(cur_mdbc.vtx_fcn)
        read_indices_with_data(vertices, scaling, input_filename, mesh, SF::NBR_REF, true, 1, PETSC_COMM_WORLD);
      else
        read_indices(vertices, input_filename, mesh, SF::NBR_REF, true, PETSC_COMM_WORLD);
      
      int gnum_idx = get_global(vertices.size(), MPI_SUM);
      if(gnum_idx == 0) {
        log_msg(0, 5, 0, "Dirichlet bcs region %d: Specified vertices are not in domain! Aborting!", idx);
        EXIT(1);
      }
    }
    else if(tag_index_given) {
      
      int tag = cur_mdbc.geomID;
      log_msg(logger, 0, 0, "Dirichlet bc region %d: Selecting vertices from tag %d", idx, tag);
      
      indices_from_region_tag(vertices, mesh, tag);
      // we restrict the indices to the algebraic subset
      SF::restrict_to_set(vertices, mesh.pl.algebraic_nodes());
    }
    else {
      log_msg(logger, 0, 0, "Dirichlet bc region %d: Selecting vertices from shape.", idx);
      
      geom_shape shape;
      shape.type   = geom_shape::shape_t(cur_mdbc.geom_type);
      shape.p0     = cur_mdbc.p0;
      shape.p1     = cur_mdbc.p1;
      shape.radius = cur_mdbc.radius;
      
      bool nodal = true;
      indices_from_geom_shape(vertices, mesh, shape, nodal);
      // we restrict the indices to the algebraic subset
      SF::restrict_to_set(vertices, mesh.pl.algebraic_nodes());
      
      int gsize = vertices.size();
      if(get_global(gsize, MPI_SUM) == 0) {
        log_msg(0,5,0, "error: Empty Dirichlet bc region [%d] def! Aborting!", idx);
        EXIT(1);
      }
    }
    if(cur_mdbc.dump_vtx_file) {
      SF::vector<mesh_int_t> glob_idx(vertices);
      mesh.pl.globalize(glob_idx);
 
      SF::vector<mesh_int_t> srt_idx;
      SF::sort_parallel(mesh.comm, glob_idx, srt_idx);
 
      size_t num_vtx = get_global(srt_idx.size(), MPI_SUM);
      int rank = get_rank();
 
      FILE_SPEC f = NULL;
      int err = 0;
 
      if(rank == 0) {
        char dump_name[1024];
        sprintf(dump_name, "DBCs_%d.vtx", idx);
 
        f = f_open(dump_name, "w");
        if(!f) err++;
        else {
          fprintf(f->fd, "%zd\nextra\n", num_vtx);
        }
      }
 
      if(!get_global(err, MPI_SUM)) {
        print_vector(mesh.comm, srt_idx, 1, f ? f->fd : NULL);
      } else {
        log_msg(0, 4, 0, "error: Dirichlet bc region [%d] cannot be dumped!");
      }
 
      // only root really does that
      f_close(f);
      }
    }
  }
}  // namespace opencarp