mechanics.h

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
// ----------------------------------------------------------------------------
 
#ifndef _MECHANICS_H
#define _MECHANICS_H
#define WITH_MECH
 
#include "physics_types.h"
#include "sim_utils.h"
#include "stimulate.h"
#include "sf_interface.h"
#include "timers.h"
#include "constitutive_model_library.h"
#include "constitutive_model.h"
 
 
namespace opencarp {
 
class mech_bcs
{
  public:
 
  enum def_t {Neumann, Dirichlet};
  enum nbc_t {pressure, force};
  enum dbc_t {all, x, y, z, xy, xz, yz};
  
  def_t definition; 
  nbc_t nbc_type; 
  dbc_t dbc_type; 
  SF::Point direction; 
  SF::vector<mesh_int_t> vertices; 
  SF::vector<SF_real>  scaling; 
  std::string  input_filename; 
 
  void setup(int idx);
};
 
class equilibrium_solver
{
  public:
  sf_vec* position; 
  sf_vec* reference_pos; 
  sf_vec* displacement; 
  sf_vec* du; 
  sf_vec* initial_guess; 
  sf_vec* f_internal; 
  sf_vec* f_external; 
  sf_vec* f_surf; 
  sf_vec* residuum; 
  SF::vector<mesh_int_t> dbc_idx; 
  SF::vector<mesh_int_t> nbc_idx; 
  sf_vec* dirichlet_bcs; 
  sf_vec* dirichlet_bcs_1; 
  sf_mat* K_eq        = nullptr; 
  
  mech_bcs dmbcs[10]; 
  mech_bcs nmbcs[10]; 
 
  sf_nl_sol* nonlin_solver = nullptr; 
    
  // solver config
  double tol = 1e-8; 
  int max_it = 300;  
  
  void init();
  void rebuild_matrices(MaterialType mtype, FILE_SPEC logger);
  void rebuild_stiffness(MaterialType mtype, FILE_SPEC logger);
  void calc_internal_forces(MaterialType mtype, FILE_SPEC logger);
  void calc_internal_forces_jacobian(MaterialType mtype, FILE_SPEC logger);
  void apply_external_forces(MaterialType mtype, FILE_SPEC logger, sf_mesh surfmesh, mech_bcs nmbcs);
  void solve();
  inline void apply_bcs_forces()
  {
    double * f_int = f_internal->ptr();
    double * f_ext = f_external->ptr();
    double * dbcs = dirichlet_bcs->ptr();
    for(int i = 0; i<dirichlet_bcs->lsize(); i++)
    {
      if(!dbcs[i])
      {
        f_int[i] = 0.0;
        f_ext[i] = 0.0;
      }
    }
    f_internal->release_ptr(f_int);
    f_external->release_ptr(f_ext);
    dirichlet_bcs->release_ptr(dbcs);
  }
  
  inline void apply_bcs_jacobian()
  {
    K_eq->mult_LR(*dirichlet_bcs, *dirichlet_bcs);
    K_eq->diag_add(*dirichlet_bcs_1);
  }
 
  inline void update_node_coords()
  {
    sf_mesh & mesh = get_mesh(elasticity_msh);
    const SF::vector<mesh_int_t> & nbr = mesh.get_numbering(SF::NBR_REF);
    const SF::vector<mesh_int_t> & alg_nod = mesh.pl.algebraic_nodes();
    double * pos = position->ptr();
    for(std::size_t i = 0; i < alg_nod.size(); i++)
    {
      for(std::size_t k = 0; k < position->dpn; k++)
      {
        mesh.xyz_Euler[3 * alg_nod[i] + k] = pos[3 * i + k];
      }
    }
    position->release_ptr(pos);
  }
    private:
    void setup_nonlinear_solver(FILE_SPEC logger);
};
 
class Mechanics : public Basic_physic
{
  public:
  MaterialType mtype;
  equilibrium_solver eq_solver;
  sf_mesh surfmesh;
    
  igb_output_manager output_manager;
    
  Mechanics()
  {
      name = "Mechanics";
  }
  void initialize();
 
  void destroy();
 
  // This funcs from the Basic_physic interface are currently empty
  void compute_step();
  void output_step();
 
 
  ~Mechanics();
 
  double timer_val(const int timer_id);
 
  std::string timer_unit(const int timer_id);
 
 
  private:
 
  void setup_solvers();
 
  void setup_mappings();
 
  void setup_output();
 
  void dump_matrices();
 
  void checkpointing();
    
  void register_material_model();
};
}   // namespace opencarp
 
#endif