From: Jean-Paul Pelteret Date: Tue, 21 Mar 2017 10:21:08 +0000 (+0000) Subject: Linear Elastic Active Skeletal Muscle Model X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=353c02d31d76fdd76851a7313681d4d8bd1c42f4;p=code-gallery.git Linear Elastic Active Skeletal Muscle Model In this example we present an implementation of a linearised active skeletal muscle model with application to the simulation concentric contraction of the human biceps brachii. --- diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/CMakeLists.txt b/Linear_Elastic_Active_Skeletal_Muscle_Model/CMakeLists.txt new file mode 100644 index 0000000..452c826 --- /dev/null +++ b/Linear_Elastic_Active_Skeletal_Muscle_Model/CMakeLists.txt @@ -0,0 +1,39 @@ +## +# CMake script for the step-8 tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "Linear_active_muscle_model") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) + +FIND_PACKAGE(deal.II 8.5 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/Linear_active_muscle_model.cc b/Linear_Elastic_Active_Skeletal_Muscle_Model/Linear_active_muscle_model.cc new file mode 100644 index 0000000..c632083 --- /dev/null +++ b/Linear_Elastic_Active_Skeletal_Muscle_Model/Linear_active_muscle_model.cc @@ -0,0 +1,1972 @@ +/* --------------------------------------------------------------------- + * Copyright (C) 2017 by the deal.II authors and + * Jean-Paul Pelteret and Tim Hamann + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + */ + +/* + * Authors: Jean-Paul Pelteret, Tim Hamann, + * University of Erlangen-Nuremberg, 2017 + * + * The support of this work by the European Research Council (ERC) through + * the Advanced Grant 289049 MOCOPOLY is gratefully acknowledged by the + * first author. + */ + +// @sect3{Include files} + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace LMM +{ + using namespace dealii; + +// @sect3{Run-time parameters} +// +// There are several parameters that can be set in the code so we set up a +// ParameterHandler object to read in the choices at run-time. + namespace Parameters + { +// @sect4{Finite Element system} + +// Here we specify the polynomial order used to approximate the solution. +// The quadrature order should be adjusted accordingly. + struct FESystem + { + unsigned int poly_degree; + unsigned int quad_order; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void FESystem::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Finite element system"); + { + prm.declare_entry("Polynomial degree", "1", + Patterns::Integer(0), + "Displacement system polynomial order"); + + prm.declare_entry("Quadrature order", "2", + Patterns::Integer(0), + "Gauss quadrature order"); + } + prm.leave_subsection(); + } + + void FESystem::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Finite element system"); + { + poly_degree = prm.get_integer("Polynomial degree"); + quad_order = prm.get_integer("Quadrature order"); + } + prm.leave_subsection(); + } + +// @sect4{Problem} + +// Choose which problem is going to be solved + struct Problem + { + std::string problem; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void Problem::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Problem"); + { + prm.declare_entry("Problem", "IsotonicContraction", + Patterns::Selection("IsotonicContraction|BicepsBrachii"), + "The problem that is to be solved"); + } + prm.leave_subsection(); + } + + void Problem::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Problem"); + { + problem = prm.get("Problem"); + } + prm.leave_subsection(); + } + +// @sect4{IsotonicContractionGeometry} + +// Make adjustments to the geometry and discretisation of the +// isotonic contraction model from Martins2006. + + struct IsotonicContraction + { + const double half_length_x = 10e-3/2.0; + const double half_length_y = 10e-3/2.0; + const double half_length_z = 1e-3/2.0; + const types::boundary_id bid_CC_dirichlet_symm_X = 1; + const types::boundary_id bid_CC_dirichlet_symm_Z = 2; + const types::boundary_id bid_CC_neumann = 10; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void IsotonicContraction::declare_parameters(ParameterHandler &prm) + { + + } + + void IsotonicContraction::parse_parameters(ParameterHandler &prm) + { + + } + +// @sect4{BicepsBrachiiGeometry} + +// Make adjustments to the geometry and discretisation of the +// biceps model. + + struct BicepsBrachii + { + double axial_length; + double radius_insertion_origin; + double radius_midpoint; + double scale; + unsigned int elements_along_axis; + unsigned int n_refinements_radial; + bool include_gravity; + double axial_force; + + const types::boundary_id bid_BB_dirichlet_X = 1; + const types::boundary_id bid_BB_neumann = 10; + const types::boundary_id mid_BB_radial = 100; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void BicepsBrachii::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Biceps Brachii geometry"); + { + prm.declare_entry("Axial length", "250", + Patterns::Double(0), + "Axial length of the muscle"); + + prm.declare_entry("Radius insertion and origin", "5", + Patterns::Double(0), + "Insertion and origin radius"); + + prm.declare_entry("Radius midpoint", "7.5", + Patterns::Double(0), + "Radius at the midpoint of the muscle"); + + prm.declare_entry("Grid scale", "1e-3", + Patterns::Double(0.0), + "Global grid scaling factor"); + + prm.declare_entry("Elements along axis", "32", + Patterns::Integer(2), + "Number of elements along the muscle axis"); + + prm.declare_entry("Radial refinements", "4", + Patterns::Integer(0), + "Control the discretisation in the radial direction"); + + prm.declare_entry("Gravity", "false", + Patterns::Bool(), + "Include the effects of gravity (in the y-direction; " + " perpendicular to the muscle axis)"); + + prm.declare_entry("Axial force", "1", + Patterns::Double(), + "Applied distributed axial force (in Newtons)"); + } + prm.leave_subsection(); + } + + void BicepsBrachii::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Biceps Brachii geometry"); + { + axial_length = prm.get_double("Axial length"); + radius_insertion_origin = prm.get_double("Radius insertion and origin"); + radius_midpoint = prm.get_double("Radius midpoint"); + scale = prm.get_double("Grid scale"); + elements_along_axis = prm.get_integer("Elements along axis"); + n_refinements_radial = prm.get_integer("Radial refinements"); + include_gravity = prm.get_bool("Gravity"); + axial_force = prm.get_double("Axial force"); + } + prm.leave_subsection(); + + AssertThrow(radius_midpoint >= radius_insertion_origin, + ExcMessage("Unrealistic geometry")); + } + +// @sect4{Neurological signal} + + struct NeurologicalSignal + { + double neural_signal_start_time; + double neural_signal_end_time; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void NeurologicalSignal::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Neurological signal"); + { + prm.declare_entry("Start time", "1.0", + Patterns::Double(0), + "Time at which to start muscle activation"); + + prm.declare_entry("End time", "2.0", + Patterns::Double(0), + "Time at which to remove muscle activation signal"); + } + prm.leave_subsection(); + } + + void NeurologicalSignal::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Neurological signal"); + { + neural_signal_start_time = prm.get_double("Start time"); + neural_signal_end_time = prm.get_double("End time"); + } + prm.leave_subsection(); + + Assert(neural_signal_start_time < neural_signal_end_time, + ExcMessage("Invalid neural signal times.")); + } + +// @sect4{Time} + +// Set the timestep size $ \varDelta t $ and the simulation end-time. + struct Time + { + double delta_t; + double end_time; + double end_ramp_time; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void Time::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + prm.declare_entry("End time", "3", + Patterns::Double(0), + "End time"); + + prm.declare_entry("End ramp time", "1", + Patterns::Double(0), + "Force ramp end time"); + + prm.declare_entry("Time step size", "0.1", + Patterns::Double(0), + "Time step size"); + } + prm.leave_subsection(); + } + + void Time::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + end_time = prm.get_double("End time"); + end_ramp_time = prm.get_double("End ramp time"); + delta_t = prm.get_double("Time step size"); + } + prm.leave_subsection(); + } + +// @sect4{All parameters} + +// Finally we consolidate all of the above structures into a single container +// that holds all of our run-time selections. + struct AllParameters : public FESystem, + public Problem, + public IsotonicContraction, + public BicepsBrachii, + public NeurologicalSignal, + public Time + { + AllParameters(const std::string &input_file); + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + AllParameters::AllParameters(const std::string &input_file) + { + ParameterHandler prm; + declare_parameters(prm); + prm.parse_input(input_file); + parse_parameters(prm); + } + + void AllParameters::declare_parameters(ParameterHandler &prm) + { + FESystem::declare_parameters(prm); + Problem::declare_parameters(prm); + IsotonicContraction::declare_parameters(prm); + BicepsBrachii::declare_parameters(prm); + NeurologicalSignal::declare_parameters(prm); + Time::declare_parameters(prm); + } + + void AllParameters::parse_parameters(ParameterHandler &prm) + { + FESystem::parse_parameters(prm); + Problem::parse_parameters(prm); + IsotonicContraction::parse_parameters(prm); + BicepsBrachii::parse_parameters(prm); + NeurologicalSignal::parse_parameters(prm); + Time::parse_parameters(prm); + + // Override time setting for test defined + // in the literature + if (problem == "IsotonicContraction") + { + end_time = 3.0; + end_ramp_time = 1.0; + delta_t = 0.1; + + neural_signal_start_time = 1.0; + neural_signal_end_time = 2.0; + } + } + } + + // @sect3{Body force values} + + template + class BodyForce : public Function + { + public: + BodyForce (const double rho, + const Tensor<1,dim> direction); + virtual ~BodyForce () {} + + virtual void vector_value (const Point &p, + Vector &values) const; + + virtual void vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + + const double rho; + const double g; + const Tensor<1,dim> M; + }; + + + template + BodyForce::BodyForce (const double rho, + const Tensor<1,dim> direction) + : + Function (dim), + rho (rho), + g (9.81), + M (direction) + { + Assert(M.norm() == 1.0, ExcMessage("Direction vector is not a unit vector")); + } + + + template + inline + void BodyForce::vector_value (const Point &p, + Vector &values) const + { + Assert (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + Assert (dim >= 2, ExcNotImplemented()); + for (unsigned int d=0; d + void BodyForce::vector_value_list (const std::vector > &points, + std::vector > &value_list) const + { + Assert (value_list.size() == points.size(), + ExcDimensionMismatch (value_list.size(), points.size())); + + const unsigned int n_points = points.size(); + + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); + } + + template + class Traction : public Function + { + public: + Traction (const double force, + const double area); + virtual ~Traction () {} + + virtual void vector_value (const Point &p, + Vector &values) const; + + virtual void vector_value_list (const std::vector > &points, + std::vector > &value_list) const; + + const double t; + }; + + + template + Traction::Traction (const double force, + const double area) + : + Function (dim), + t (force/area) + {} + + + template + inline + void Traction::vector_value (const Point &p, + Vector &values) const + { + Assert (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + Assert (dim == 3, ExcNotImplemented()); + + // Assume uniform distributed load + values(0) = t; + values(1) = 0.0; + values(2) = 0.0; + } + + + template + void Traction::vector_value_list (const std::vector > &points, + std::vector > &value_list) const + { + Assert (value_list.size() == points.size(), + ExcDimensionMismatch (value_list.size(), points.size())); + + const unsigned int n_points = points.size(); + + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); + } + + // @sect3{Utility functions} + + template + inline + Tensor<2,dim> get_deformation_gradient (std::vector > &grad) + { + Assert (grad.size() == dim, ExcInternalError()); + + Tensor<2,dim> F (unit_symmetric_tensor()); + for (unsigned int i=0; i + inline + SymmetricTensor<2,dim> get_small_strain (std::vector > &grad) + { + Assert (grad.size() == dim, ExcInternalError()); + + SymmetricTensor<2,dim> strain; + for (unsigned int i=0; i + struct MuscleFibre + { + MuscleFibre (void) + : alpha (0.0), + alpha_t1 (0.0), + epsilon_f (0.0), + epsilon_c (0.0), + epsilon_c_t1 (0.0), + epsilon_c_dot (0.0) + { + + } + + MuscleFibre(const Tensor<1,dim> &direction) + : M (direction), + alpha (0.0), + alpha_t1 (0.0), + epsilon_f (0.0), + epsilon_c (0.0), + epsilon_c_t1 (0.0), + epsilon_c_dot (0.0) + { + Assert(M.norm() == 1.0, + ExcMessage("Fibre direction is not a unit vector")); + } + + void update_alpha (const double u, + const double dt); + + void update_state(const SymmetricTensor<2,dim> &strain_tensor, + const double dt); + + const Tensor<1,dim> &get_M () const + { + return M; + } + double get_m_p () const; + double get_m_s () const; + double get_beta (const double dt) const; + double get_gamma (const double dt) const; + + // Postprocessing + const double &get_alpha() const + { + return alpha; + } + const double &get_epsilon_f() const + { + return epsilon_f; + } + const double &get_epsilon_c() const + { + return epsilon_c; + } + const double &get_epsilon_c_dot() const + { + return epsilon_c_dot; + } + + private: + Tensor<1,dim> M; // Direction + + double alpha; // Activation level at current timestep + double alpha_t1; // Activation level at previous timestep + + double epsilon_f; // Fibre strain at current timestep + double epsilon_c; // Contractile strain at current timestep + double epsilon_c_t1; // Contractile strain at previous timestep + double epsilon_c_dot; // Contractile velocity at previous timestep + + double get_f_c_L () const; + double get_m_c_V () const; + double get_c_c_V () const; + }; + + template + void MuscleFibre::update_alpha (const double u, + const double dt) + { + static const double tau_r = 0.15; // s + static const double tau_f = 0.15; // s + static const double alpha_min = 0; + + if (u == 1.0) + alpha = (alpha_t1*tau_r*tau_f + dt*tau_f) / (tau_r*tau_f + dt*tau_f); + else if (u == 0) + alpha = (alpha_t1*tau_r*tau_f + dt*alpha_min*tau_r) / (tau_r*tau_f + dt*tau_r); + else + { + const double b = 1.0/tau_r - 1.0/tau_f; + const double c = 1.0/tau_f; + const double d = alpha_min/tau_f; + const double f1 = 1.0/tau_r - alpha_min/tau_f; + const double p = b*u + c; + const double q = f1*u + d; + + alpha = (q*dt + alpha_t1)/(1.0 + p*dt); + } + } + + + template + double MuscleFibre::get_m_p () const + { + static const double A = 8.568e-4*convert_gf_per_cm2_to_N_per_m2; + static const double a = 12.43; + if (epsilon_f >= 0.0) + { + // 100 times more compliant than Martins2006 + static const double m_p = 2.0*A*a/1e2; + return m_p; + } + else + return 0.0; + } + + template + double MuscleFibre::get_m_s (void) const + { + const double epsilon_s = epsilon_f - epsilon_c; // Small strain assumption + if (epsilon_s >= -1e-6) // Tolerant check + return 10.0; + else + return 0.0; + } + + template + double MuscleFibre::get_f_c_L (void) const + { + if (epsilon_c <= 0.5 && epsilon_c >= -0.5) + return 1.0; + else + return 0.0; + } + + template + double MuscleFibre::get_m_c_V (void) const + { + if (epsilon_c_dot < -5.0) + return 0.0; + else if (epsilon_c_dot <= 3.0) + return 1.0/5.0; + else + return 0.0; + } + + template + double MuscleFibre::get_c_c_V (void) const + { + if (epsilon_c_dot < -5.0) + return 0.0; + else if (epsilon_c_dot <= 3.0) + return 1.0; + else + return 1.6; + } + + template + double MuscleFibre::get_beta(const double dt) const + { + return get_f_c_L()*get_m_c_V()*alpha/dt + get_m_s(); + } + + template + double MuscleFibre::get_gamma(const double dt) const + { + return get_f_c_L()*alpha*(get_m_c_V()*epsilon_c_t1/dt - get_c_c_V()); + } + + template + void MuscleFibre::update_state(const SymmetricTensor<2,dim> &strain_tensor, + const double dt) + { + // Values from previous state + // These were the values that were used in the assembly, + // so we must use them in the update step to be consistant. + // Need to compute these before we overwrite epsilon_c_t1 + const double m_s = get_m_s(); + const double beta = get_beta(dt); + const double gamma = get_gamma(dt); + + // Update current state + alpha_t1 = alpha; + epsilon_f = M*static_cast< Tensor<2,dim> >(strain_tensor)*M; + epsilon_c_t1 = epsilon_c; + epsilon_c = (m_s*epsilon_f + gamma)/beta; + epsilon_c_dot = (epsilon_c - epsilon_c_t1)/dt; + } + + + // @sect3{The LinearMuscleModelProblem class template} + + template + class LinearMuscleModelProblem + { + public: + LinearMuscleModelProblem (const std::string &input_file); + ~LinearMuscleModelProblem (); + void run (); + + private: + void make_grid (); + void setup_muscle_fibres (); + double get_neural_signal (const double time); + void update_fibre_activation (const double time); + void update_fibre_state (); + void setup_system (); + void assemble_system (const double time); + void apply_boundary_conditions (); + void solve (); + void output_results (const unsigned int timestep, + const double time) const; + + Parameters::AllParameters parameters; + + Triangulation triangulation; + DoFHandler dof_handler; + + FESystem fe; + QGauss qf_cell; + QGauss qf_face; + + ConstraintMatrix hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + // Time + const double t_end; + const double dt; + const double t_ramp_end; // Force ramp end time + + // Loading + const BodyForce body_force; + const Traction traction; + + // Local data + std::vector< std::vector > > fibre_data; + + // Constitutive functions for assembly + SymmetricTensor<4,dim> get_stiffness_tensor (const unsigned int cell, + const unsigned int q_point_cell) const; + SymmetricTensor<2,dim> get_rhs_tensor (const unsigned int cell, + const unsigned int q_point_cell) const; + }; + + // @sect4{LinearMuscleModelProblem::LinearMuscleModelProblem} + + template + LinearMuscleModelProblem::LinearMuscleModelProblem (const std::string &input_file) + : + parameters(input_file), + dof_handler (triangulation), + fe (FE_Q(parameters.poly_degree), dim), + qf_cell (parameters.quad_order), + qf_face (parameters.quad_order), + t_end (parameters.end_time), + dt (parameters.delta_t), + t_ramp_end(parameters.end_ramp_time), + body_force ((parameters.problem == "BicepsBrachii" &¶meters.include_gravity == true) ? + BodyForce(0.375*1000.0, Tensor<1,dim>({0,-1,0})) : // (reduced) Density and direction + BodyForce(0.0, Tensor<1,dim>({0,0,1})) ), + traction (parameters.problem == "BicepsBrachii" ? + Traction(parameters.axial_force, // Force, area + M_PI*std::pow(parameters.radius_insertion_origin *parameters.scale,2.0) ) : + Traction(4.9*convert_gf_to_N, // Force; Conversion of gf to N, + (2.0*parameters.half_length_y)*(2.0*parameters.half_length_z)) ) // Area + { + Assert(dim==3, ExcNotImplemented()); + } + + + // @sect4{LinearMuscleModelProblem::~LinearMuscleModelProblem} + + template + LinearMuscleModelProblem::~LinearMuscleModelProblem () + { + dof_handler.clear (); + } + + + // @sect4{LinearMuscleModelProblem::make_grid} + + template + struct BicepsGeometry + { + BicepsGeometry(const double axial_length, + const double radius_ins_orig, + const double radius_mid) + : + ax_lgth (axial_length), + r_ins_orig (radius_ins_orig), + r_mid (radius_mid) + {} + + // The radial profile of the muscle + // This provides the new coordinates for points @p pt + // on a cylinder of radius r_ins_orig and length + // ax_lgth to be moved to in order to create the + // physiologically representative geometry of + // the muscle + Point profile (const Point &pt_0) const + { + Assert(pt[0] > -1e-6, + ExcMessage("All points must have x-coordinate > 0")); + + const double r_scale = get_radial_scaling_factor(pt_0[0]); + return pt_0 + Point(0.0, r_scale*pt_0[1], r_scale*pt_0[2]); + } + + Point operator() (const Point &pt) const + { + return profile(pt); + } + + // Provides the muscle direction at the point @p pt + // in the real geometry (one that has undergone the + // transformation given by the profile() function) + // and subequent grid rescaling. + // The directions are given by the gradient of the + // transformation function (i.e. the fibres are + // orientated by the curvature of the muscle). + // + // So, being lazy, we transform the current point back + // to the original point on the completely unscaled + // cylindrical grid. We then evaluate the transformation + // at two points (axially displaced) very close to the + // point of interest. The normalised vector joining the + // transformed counterparts of the perturbed points is + // the gradient of the transformation function and, + // thus, defines the fibre direction. + Tensor<1,dim> direction (const Point &pt_scaled, + const double &grid_scale) const + { + const Point pt = (1.0/grid_scale)*pt_scaled; + const Point pt_0 = inv_profile(pt); + + static const double eps = 1e-6; + const Point pt_0_eps_p = pt_0 + Point(+eps,0,0); + const Point pt_0_eps_m = pt_0 + Point(-eps,0,0); + const Point pt_eps_p = profile(pt_0_eps_p); + const Point pt_eps_m = profile(pt_0_eps_m); + + static const double tol = 1e-9; + Assert(profile(pt_0).distance(pt) < tol, ExcInternalError()); + Assert(inv_profile(pt_eps_p).distance(pt_0_eps_p) < tol, ExcInternalError()); + Assert(inv_profile(pt_eps_m).distance(pt_0_eps_m) < tol, ExcInternalError()); + + Tensor<1,dim> dir = pt_eps_p-pt_eps_m; + dir /= dir.norm(); + return dir; + } + + private: + const double ax_lgth; + const double r_ins_orig; + const double r_mid; + + double get_radial_scaling_factor (const double &x) const + { + // Expect all grid points with X>=0, but we provide a + // tolerant location for points "on" the Cartesian plane X=0 + const double lgth_frac = std::max(x/ax_lgth,0.0); + const double amplitude = 0.25*(r_mid - r_ins_orig); + const double phase_shift = M_PI; + const double y_shift = 1.0; + const double wave_func = y_shift + std::cos(phase_shift + 2.0*M_PI*lgth_frac); + Assert(wave_func >= 0.0, ExcInternalError()); + return std::sqrt(amplitude*wave_func); + } + + Point inv_profile (const Point &pt) const + { + Assert(pt[0] > -1e-6, + ExcMessage("All points must have x-coordinate > 0")); + + const double r_scale = get_radial_scaling_factor(pt[0]); + const double trans_inv_scale = 1.0/(1.0+r_scale); + return Point(pt[0], trans_inv_scale*pt[1], trans_inv_scale*pt[2]); + } + }; + + template + void LinearMuscleModelProblem::make_grid () + { + Assert (dim == 3, ExcNotImplemented()); + + if (parameters.problem == "IsotonicContraction") + { + const Point p1(-parameters.half_length_x, + -parameters.half_length_y, + -parameters.half_length_z); + const Point p2( parameters.half_length_x, + parameters.half_length_y, + parameters.half_length_z); + + GridGenerator::hyper_rectangle (triangulation, p1, p2); + + typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(), endc = triangulation.end(); + for (; cell != endc; ++cell) + { + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() == true) + { + if (cell->face(face)->center()[0] == -parameters.half_length_x) // -X oriented face + cell->face(face)->set_boundary_id(parameters.bid_CC_dirichlet_symm_X); // Dirichlet + else if (cell->face(face)->center()[0] == parameters.half_length_x) // +X oriented face + cell->face(face)->set_boundary_id(parameters.bid_CC_neumann); // Neumann + else if (std::abs(cell->face(face)->center()[2]) == parameters.half_length_z) // -Z/+Z oriented face + cell->face(face)->set_boundary_id(parameters.bid_CC_dirichlet_symm_Z); // Dirichlet + } + } + } + + triangulation.refine_global (1); + } + else if (parameters.problem == "BicepsBrachii") + { + SphericalManifold<2> manifold_cap; + Triangulation<2> tria_cap; + GridGenerator::hyper_ball(tria_cap, + Point<2>(), + parameters.radius_insertion_origin); + for (typename Triangulation<2>::active_cell_iterator + cell = tria_cap.begin_active(); + cell != tria_cap.end(); ++cell) + { + for (unsigned int face = 0; face < GeometryInfo<2>::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() == true) + cell->face(face)->set_all_manifold_ids(0); + } + } + tria_cap.set_manifold (0, manifold_cap); + tria_cap.refine_global(parameters.n_refinements_radial); + + Triangulation<2> tria_cap_flat; + GridGenerator::flatten_triangulation(tria_cap, tria_cap_flat); + + GridGenerator::extrude_triangulation(tria_cap_flat, + parameters.elements_along_axis, + parameters.axial_length, + triangulation); + + struct GridRotate + { + Point operator() (const Point &in) const + { + static const Tensor<2,dim> rot_mat = Physics::Transformations::Rotations::rotation_matrix_3d(Point(0,1,0), M_PI/2.0); + return Point(rot_mat*in); + } + }; + + // Rotate grid so that the length is axially + // coincident and aligned with the X-axis + GridTools::transform (GridRotate(), triangulation); + + // Deform the grid into something that vaguely + // resemble's a Biceps Brachii + GridTools::transform (BicepsGeometry(parameters.axial_length, + parameters.radius_insertion_origin, + parameters.radius_midpoint), triangulation); + + // Set boundary IDs + typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(), endc = triangulation.end(); + for (; cell != endc; ++cell) + { + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() == true) + { + static const double tol =1e-6; + if (std::abs(cell->face(face)->center()[0]) < tol) // -X oriented face + cell->face(face)->set_boundary_id(parameters.bid_BB_dirichlet_X); // Dirichlet + else if (std::abs(cell->face(face)->center()[0] - parameters.axial_length) < tol) // +X oriented face + cell->face(face)->set_boundary_id(parameters.bid_BB_neumann); // Neumann + } + } + } + + // Finally resize the grid + GridTools::scale (parameters.scale, triangulation); + } + else + AssertThrow(false, ExcNotImplemented()); + } + + // @sect4{LinearMuscleModelProblem::setup_muscle_fibres} + + template + void LinearMuscleModelProblem::setup_muscle_fibres () + { + fibre_data.clear(); + const unsigned int n_cells = triangulation.n_active_cells(); + fibre_data.resize(n_cells); + const unsigned int n_q_points_cell = qf_cell.size(); + + if (parameters.problem == "IsotonicContraction") + { + MuscleFibre fibre_template (Tensor<1,dim>({1,0,0})); + + for (unsigned int cell_no=0; cell_no fe_values (fe, qf_cell, update_quadrature_points); + BicepsGeometry bicep_geom (parameters.axial_length, + parameters.radius_insertion_origin, + parameters.radius_midpoint); + + unsigned int cell_no = 0; + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); + ++cell, ++cell_no) + { + Assert(cell_no pt = fe_values.get_quadrature_points()[q_point_cell]; + fibre_data[cell_no][q_point_cell] = MuscleFibre(bicep_geom.direction(pt,parameters.scale)); + } + } + } + else + AssertThrow(false, ExcNotImplemented()); + } + + // @sect4{LinearMuscleModelProblem::update_fibre_state} + + template + double LinearMuscleModelProblem::get_neural_signal (const double time) + { + // Note: 40 times less force generated than Martins2006 + // This is necessary due to the (compliant) linear tissue model + return (time > parameters.neural_signal_start_time && time < parameters.neural_signal_end_time ? + 1.0/40.0 : + 0.0); + } + + template + void LinearMuscleModelProblem::update_fibre_activation (const double time) + { + const double u = get_neural_signal(time); + + const unsigned int n_cells = triangulation.n_active_cells(); + const unsigned int n_q_points_cell = qf_cell.size(); + for (unsigned int cell=0; cell &fibre = fibre_data[cell][q_point_cell]; + fibre.update_alpha(u,dt); + } + } + } + + template + void LinearMuscleModelProblem::update_fibre_state () + { + const unsigned int n_cells = triangulation.n_active_cells(); + const unsigned int n_q_points_cell = qf_cell.size(); + + FEValues fe_values (fe, qf_cell, update_gradients); + + // Displacement gradient + std::vector< std::vector< Tensor<1,dim> > > u_grads (n_q_points_cell, + std::vector >(dim)); + + unsigned int cell_no = 0; + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell!=dof_handler.end(); ++cell, ++cell_no) + { + Assert(cell_no strain_tensor = get_small_strain (u_grads[q_point_cell]); + MuscleFibre &fibre = fibre_data[cell_no][q_point_cell]; + fibre.update_state(strain_tensor, dt); + } + } + } + + // @sect4{LinearMuscleModelProblem::setup_system} + + template + void LinearMuscleModelProblem::setup_system () + { + dof_handler.distribute_dofs (fe); + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + hanging_node_constraints.close (); + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + + hanging_node_constraints.condense (sparsity_pattern); + + sparsity_pattern.compress(); + + system_matrix.reinit (sparsity_pattern); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + } + + // @sect4{LinearMuscleModelProblem::assemble_system} + + template + SymmetricTensor<4,dim> + LinearMuscleModelProblem::get_stiffness_tensor (const unsigned int cell, + const unsigned int q_point_cell) const + { + static const SymmetricTensor<2,dim> I = unit_symmetric_tensor(); + + Assert(cell &fibre = fibre_data[cell][q_point_cell]; + + // Matrix + const double lambda = MuscleMatrix::lambda; + const double mu = MuscleMatrix::mu; + // Fibre + const double m_p = fibre.get_m_p(); + const double m_s = fibre.get_m_s(); + const double beta = fibre.get_beta(dt); + const double gamma = fibre.get_gamma(dt); + AssertThrow(beta != 0.0, ExcInternalError()); + const double Cf = T0*(m_p + m_s*(1.0 - m_s/beta)); + const Tensor<1,dim> &M = fibre.get_M(); + + SymmetricTensor<4,dim> C; + for (unsigned int i=0; i < dim; ++i) + for (unsigned int j=i; j < dim; ++j) + for (unsigned int k=0; k < dim; ++k) + for (unsigned int l=k; l < dim; ++l) + { + // Matrix contribution + C[i][j][k][l] = lambda * I[i][j]*I[k][l] + + mu * (I[i][k]*I[j][l] + I[i][l]*I[j][k]); + + // Fibre contribution (Passive + active branches) + C[i][j][k][l] += Cf * M[i]*M[j]*M[k]*M[l]; + } + + return C; + } + + template + SymmetricTensor<2,dim> + LinearMuscleModelProblem::get_rhs_tensor (const unsigned int cell, + const unsigned int q_point_cell) const + { + Assert(cell &fibre = fibre_data[cell][q_point_cell]; + + const double m_s = fibre.get_m_s(); + const double beta = fibre.get_beta(dt); + const double gamma = fibre.get_gamma(dt); + AssertThrow(beta != 0.0, ExcInternalError()); + const double Sf = T0*(m_s*gamma/beta); + const Tensor<1,dim> &M = fibre.get_M(); + + SymmetricTensor<2,dim> S; + for (unsigned int i=0; i < dim; ++i) + for (unsigned int j=i; j < dim; ++j) + { + // Fibre contribution (Active branch) + S[i][j] = Sf * M[i]*M[j]; + } + + return S; + } + + // @sect4{LinearMuscleModelProblem::assemble_system} + + template + void LinearMuscleModelProblem::assemble_system (const double time) + { + // Reset system + system_matrix = 0; + system_rhs = 0; + + FEValues fe_values (fe, qf_cell, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + FEFaceValues fe_face_values (fe, qf_face, + update_values | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points_cell = qf_cell.size(); + const unsigned int n_q_points_face = qf_face.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + // Loading + std::vector > body_force_values (n_q_points_cell, + Vector(dim)); + std::vector > traction_values (n_q_points_face, + Vector(dim)); + + unsigned int cell_no = 0; + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell!=dof_handler.end(); ++cell, ++cell_no) + { + cell_matrix = 0; + cell_rhs = 0; + + fe_values.reinit (cell); + body_force.vector_value_list (fe_values.get_quadrature_points(), + body_force_values); + + for (unsigned int q_point_cell=0; q_point_cell C = get_stiffness_tensor (cell_no, q_point_cell); + const SymmetricTensor<2,dim> R = get_rhs_tensor(cell_no, q_point_cell); + + for (unsigned int I=0; I::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() == true && + ((parameters.problem == "IsotonicContraction" && + cell->face(face)->boundary_id() == parameters.bid_CC_neumann) || + (parameters.problem == "BicepsBrachii" && + cell->face(face)->boundary_id() == parameters.bid_BB_neumann)) ) + { + fe_face_values.reinit(cell, face); + traction.vector_value_list (fe_face_values.get_quadrature_points(), + traction_values); + + // Scale applied traction according to time + const double ramp = (time <= t_ramp_end ? time/t_ramp_end : 1.0); + Assert(ramp >= 0.0 && ramp <= 1.0, ExcMessage("Invalid force ramp")); + for (unsigned int q_point_face = 0; q_point_face < n_q_points_face; ++q_point_face) + traction_values[q_point_face] *= ramp; + + for (unsigned int q_point_face = 0; q_point_face < n_q_points_face; ++q_point_face) + { + for (unsigned int I=0; Iget_dof_indices (local_dof_indices); + for (unsigned int i=0; i + void LinearMuscleModelProblem::apply_boundary_conditions () + { + std::map boundary_values; + + if (parameters.problem == "IsotonicContraction") + { + // Symmetry condition on -X faces + { + ComponentMask component_mask_x (dim, false); + component_mask_x.set(0, true); + VectorTools::interpolate_boundary_values (dof_handler, + parameters.bid_CC_dirichlet_symm_X, + ZeroFunction(dim), + boundary_values, + component_mask_x); + } + // Symmetry condition on -Z/+Z faces + { + ComponentMask component_mask_z (dim, false); + component_mask_z.set(2, true); + VectorTools::interpolate_boundary_values (dof_handler, + parameters.bid_CC_dirichlet_symm_Z, + ZeroFunction(dim), + boundary_values, + component_mask_z); + } + // Fixed point on -X face + { + const Point fixed_point (-parameters.half_length_x,0.0,0.0); + std::vector fixed_dof_indices; + bool found_point_of_interest = false; + types::global_dof_index dof_of_interest = numbers::invalid_dof_index; + + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); cell != endc; ++cell) + { + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + { + // We know that the fixed point is on the -X Dirichlet boundary + if (cell->face(face)->at_boundary() == true && + cell->face(face)->boundary_id() == parameters.bid_CC_dirichlet_symm_X) + { + for (unsigned int face_vertex_index = 0; face_vertex_index < GeometryInfo::vertices_per_face; ++face_vertex_index) + { + if (cell->face(face)->vertex(face_vertex_index).distance(fixed_point) < 1e-6) + { + found_point_of_interest = true; + for (unsigned int index_component = 0; index_component < dim; ++index_component) + fixed_dof_indices.push_back(cell->face(face)->vertex_dof_index(face_vertex_index, + index_component)); + } + + if (found_point_of_interest == true) break; + } + } + if (found_point_of_interest == true) break; + } + if (found_point_of_interest == true) break; + } + + Assert(found_point_of_interest == true, ExcMessage("Didn't find point of interest")); + AssertThrow(fixed_dof_indices.size() == dim, ExcMessage("Didn't find the correct number of DoFs to fix")); + + for (unsigned int i=0; i < fixed_dof_indices.size(); ++i) + boundary_values[fixed_dof_indices[i]] = 0.0; + } + } + else if (parameters.problem == "BicepsBrachii") + { + if (parameters.include_gravity == false) + { + // Symmetry condition on -X surface + { + ComponentMask component_mask_x (dim, false); + component_mask_x.set(0, true); + VectorTools::interpolate_boundary_values (dof_handler, + parameters.bid_BB_dirichlet_X, + ZeroFunction(dim), + boundary_values, + component_mask_x); + } + + // Fixed central point on -X surface + { + const Point fixed_point (0.0,0.0,0.0); + std::vector fixed_dof_indices; + bool found_point_of_interest = false; + types::global_dof_index dof_of_interest = numbers::invalid_dof_index; + + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); cell != endc; ++cell) + { + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + { + // We know that the fixed point is on the -X Dirichlet boundary + if (cell->face(face)->at_boundary() == true && + cell->face(face)->boundary_id() == parameters.bid_BB_dirichlet_X) + { + for (unsigned int face_vertex_index = 0; face_vertex_index < GeometryInfo::vertices_per_face; ++face_vertex_index) + { + if (cell->face(face)->vertex(face_vertex_index).distance(fixed_point) < 1e-6) + { + found_point_of_interest = true; + for (unsigned int index_component = 0; index_component < dim; ++index_component) + fixed_dof_indices.push_back(cell->face(face)->vertex_dof_index(face_vertex_index, + index_component)); + } + + if (found_point_of_interest == true) break; + } + } + if (found_point_of_interest == true) break; + } + if (found_point_of_interest == true) break; + } + + Assert(found_point_of_interest == true, ExcMessage("Didn't find point of interest")); + AssertThrow(fixed_dof_indices.size() == dim, ExcMessage("Didn't find the correct number of DoFs to fix")); + + for (unsigned int i=0; i < fixed_dof_indices.size(); ++i) + boundary_values[fixed_dof_indices[i]] = 0.0; + } + } + else + { + // When we apply gravity, some additional constraints + // are required to support the load of the muscle, as + // the material response is more compliant than would + // be the case in reality. + + // Symmetry condition on -X surface + { + ComponentMask component_mask_x (dim, true); + VectorTools::interpolate_boundary_values (dof_handler, + parameters.bid_BB_dirichlet_X, + ZeroFunction(dim), + boundary_values, + component_mask_x); + } + // Symmetry condition on -X surface + { + ComponentMask component_mask_x (dim, false); + component_mask_x.set(1, true); + component_mask_x.set(2, true); + VectorTools::interpolate_boundary_values (dof_handler, + parameters.bid_BB_neumann, + ZeroFunction(dim), + boundary_values, + component_mask_x); + } + } + + // Roller condition at central point on +X face + { + const Point roller_point (parameters.axial_length*parameters.scale,0.0,0.0); + std::vector fixed_dof_indices; + bool found_point_of_interest = false; + types::global_dof_index dof_of_interest = numbers::invalid_dof_index; + + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); cell != endc; ++cell) + { + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + { + // We know that the fixed point is on the +X Neumann boundary + if (cell->face(face)->at_boundary() == true && + cell->face(face)->boundary_id() == parameters.bid_BB_neumann) + { + for (unsigned int face_vertex_index = 0; face_vertex_index < GeometryInfo::vertices_per_face; ++face_vertex_index) + { + if (cell->face(face)->vertex(face_vertex_index).distance(roller_point) < 1e-6) + { + found_point_of_interest = true; + for (unsigned int index_component = 1; index_component < dim; ++index_component) + fixed_dof_indices.push_back(cell->face(face)->vertex_dof_index(face_vertex_index, + index_component)); + } + + if (found_point_of_interest == true) break; + } + } + if (found_point_of_interest == true) break; + } + if (found_point_of_interest == true) break; + } + + Assert(found_point_of_interest == true, ExcMessage("Didn't find point of interest")); + AssertThrow(fixed_dof_indices.size() == dim-1, ExcMessage("Didn't find the correct number of DoFs to fix")); + + for (unsigned int i=0; i < fixed_dof_indices.size(); ++i) + boundary_values[fixed_dof_indices[i]] = 0.0; + } + } + else + AssertThrow(false, ExcNotImplemented()); + + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); + } + + + // @sect4{LinearMuscleModelProblem::solve} + + template + void LinearMuscleModelProblem::solve () + { + SolverControl solver_control (system_matrix.m(), 1e-12); + SolverCG<> cg (solver_control); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + hanging_node_constraints.distribute (solution); + } + + + // @sect4{LinearMuscleModelProblem::output_results} + + + template + void LinearMuscleModelProblem::output_results (const unsigned int timestep, + const double time) const + { + // Visual output: FEM results + { + std::string filename = "solution-"; + filename += Utilities::int_to_string(timestep,4); + filename += ".vtk"; + std::ofstream output (filename.c_str()); + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + + std::vector + data_component_interpretation(dim, + DataComponentInterpretation::component_is_part_of_vector); + std::vector solution_name(dim, "displacement"); + + data_out.add_data_vector (solution, solution_name, + DataOut::type_dof_data, + data_component_interpretation); + data_out.build_patches (); + data_out.write_vtk (output); + } + + // Visual output: FEM data + { + std::string filename = "fibres-"; + filename += Utilities::int_to_string(timestep,4); + filename += ".vtk"; + std::ofstream output (filename.c_str()); + + output + << "# vtk DataFile Version 3.0" << std::endl + << "# " << std::endl + << "ASCII"<< std::endl + << "DATASET POLYDATA"<< std::endl << std::endl; + + // Extract fibre data from quadrature points + const unsigned int n_cells = triangulation.n_active_cells(); + const unsigned int n_q_points_cell = qf_cell.size(); + + // Data that we'll be outputting + std::vector results_fibre_names; + results_fibre_names.push_back("alpha"); + results_fibre_names.push_back("epsilon_f"); + results_fibre_names.push_back("epsilon_c"); + results_fibre_names.push_back("epsilon_c_dot"); + + const unsigned int n_results = results_fibre_names.size(); + const unsigned int n_data_points = n_cells*n_q_points_cell; + std::vector< Point > output_points(n_data_points); + std::vector< Tensor<1,dim> > output_displacements(n_data_points); + std::vector< Tensor<1,dim> > output_directions(n_data_points); + std::vector< std::vector > output_values(n_results, std::vector(n_data_points)); + + // Displacement + std::vector< Vector > u_values (n_q_points_cell, + Vector(dim)); + // Displacement gradient + std::vector< std::vector< Tensor<1,dim> > > u_grads (n_q_points_cell, + std::vector >(dim)); + + FEValues fe_values (fe, qf_cell, + update_values | update_gradients | update_quadrature_points); + unsigned int cell_no = 0; + unsigned int fibre_no = 0; + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); + ++cell, ++cell_no) + { + fe_values.reinit (cell); + fe_values.get_function_values (solution, u_values); + fe_values.get_function_gradients (solution, u_grads); + + for (unsigned int q_point_cell=0; q_point_cell &fibre = fibre_data[cell_no][q_point_cell]; + output_points[fibre_no] = fe_values.get_quadrature_points()[q_point_cell]; // Position + for (unsigned int d=0; d meas_pt (parameters.problem == "IsotonicContraction" ? + Point(parameters.half_length_x, 0.0, 0.0) : + Point(parameters.axial_length*parameters.scale, 0.0, 0.0) ); + + + const unsigned int index_of_interest = 0; + bool found_point_of_interest = false; + types::global_dof_index dof_of_interest = numbers::invalid_dof_index; + + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); cell != endc; ++cell) + { + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + { + // We know that the measurement point is on the Neumann boundary + if (cell->face(face)->at_boundary() == true && + ((parameters.problem == "IsotonicContraction" && + cell->face(face)->boundary_id() == parameters.bid_CC_neumann) || + (parameters.problem == "BicepsBrachii" && + cell->face(face)->boundary_id() == parameters.bid_BB_neumann)) ) + { + for (unsigned int face_vertex_index = 0; face_vertex_index < GeometryInfo::vertices_per_face; ++face_vertex_index) + { + if (cell->face(face)->vertex(face_vertex_index).distance(meas_pt) < 1e-6) + { + found_point_of_interest = true; + dof_of_interest = cell->face(face)->vertex_dof_index(face_vertex_index, + index_of_interest); + } + + if (found_point_of_interest == true) break; + } + } + if (found_point_of_interest == true) break; + } + if (found_point_of_interest == true) break; + } + + Assert(found_point_of_interest == true, ExcMessage("Didn't find point of interest")); + Assert(dof_of_interest != numbers::invalid_dof_index, ExcMessage("Didn't find DoF of interest")); + Assert(dof_of_interest < dof_handler.n_dofs(), ExcMessage("DoF index out of range")); + + const std::string filename = "displacement_POI.csv"; + std::ofstream output; + if (timestep == 0) + { + output.open(filename.c_str(), std::ofstream::out); + output + << "Time [s]" << "," << "X-displacement [mm]" << std::endl; + } + else + output.open(filename.c_str(), std::ios_base::app); + + output + << time + << "," + << solution[dof_of_interest]*1e3 + << std::endl; + } + } + + + + // @sect4{LinearMuscleModelProblem::run} + + template + void LinearMuscleModelProblem::run () + { + make_grid(); + setup_system (); + setup_muscle_fibres (); + +// const bool do_grid_refinement = false; + double time = 0.0; + for (unsigned int timestep=0; time<=t_end; ++timestep, time+=dt) + { + std::cout + << "Timestep " << timestep + << " @ time " << time + << std::endl; + + // First we update the fibre activation level + // based on the current time + update_fibre_activation(time); + + // Next we assemble the system and enforce boundary + // conditions. + // Here we assume that the system and fibres have + // a fixed state, and we will assemble based on how + // epsilon_c will update given the current state of + // the body. + assemble_system (time); + apply_boundary_conditions (); + + // Then we solve the linear system + solve (); + + // Now we update the fibre state based on the new + // displacement solution and the constitutive + // parameters assumed to govern the stiffness of + // the fibres at the previous state. i.e. We + // follow through with assumed update conditions + // used in the assembly phase. + update_fibre_state(); + + // Output some values to file + output_results (timestep, time); + } + } +} + +// @sect3{The main function} + +int main () +{ + try + { + dealii::deallog.depth_console (0); + const unsigned int dim = 3; + + LMM::LinearMuscleModelProblem lmm_problem ("parameters.prm"); + lmm_problem.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/README.md b/Linear_Elastic_Active_Skeletal_Muscle_Model/README.md new file mode 100644 index 0000000..476065b --- /dev/null +++ b/Linear_Elastic_Active_Skeletal_Muscle_Model/README.md @@ -0,0 +1,170 @@ +## Overview +The complex interaction of muscles with their surrounding anatomy and +environment plays a vital role in the many activities that are required for +animals to live and survive. +Skeletal muscle composes a large portion of that musculo-skeletal +system, and is controlled by the central nervous system in a conscious or +unconscious manner. +For humans in particular, the construction of- and control mechanisms behind +skeletal muscle allows us to accomplish complex tasks ranging from those +that are physically exerting to those that are delicate and require great +dexterity. + +As an introduction into the biomechanics of the human muscular system, we +combine a well known load-activation pattern taken from well established +literature on the topic (both in the fields of human physiology and the +computational simulation thereof) with an idealised model of a part of the +human anatomy that most can easily identify with, namely the biceps brachii. + +### An idealised model of the human biceps brachii +To tackle this problem, we do not deviate particularly far from the approach +that is comprehensively documented in `step-8`. +The primary differences between this code-gallery example and the tutorial +is the alteration of the geometry and boundary conditions that accompany it, +as well as the extension of the constitutive law to include the transversely +isotropic active muscle model. +We thus present both the theory and associated nomenclature (which is +mimicked in the code itself) in the +[accompanying document.](./doc/theory/theory-linear_elastic_active_muscle_model.pdf) +There you can observe the additional contributions to both the left- and +right-hand sides of the linear system due to the integration of the +rate-dependent, linearised muscle model. +Although such linear model models are valid only under conditions of small +deformations, in this scenario we will (with a healthy dose of skepticism) +make a very coarse-grained first approximation of the muscles behaviour in +the large deformation regime. + +The basic problem configuration, including a depiction of the underlying +muscle microstructural orientation, is (loosely) summarised in the following +image. +![Problem geometry](./doc/geometry/geometry.png) +Note that the driver for the deformation of the muscle tissue are the applied +traction alone when the muscle is in a passive state. +However, during active contraction, as governed by the prescribed input +neural signal, the muscle works against the applied traction. +This condition, where the traction applied to a muscle is constant during +periods of muscle activation, is known as isotonic contraction. +More specifically, since overall the muscle shortens during contraction we +are in fact modelling concentric contraction of the biceps. + +As for the specific geometry of the problem, we consider an idealised human +biceps with a length of `250mm`, insertion and origin diameter of `20mm` and +a diameter of `80mm` at its mid-point. +We assume that there exists a single muscle fibre family orientated axially. +The orientation of the underlying muscle fibres is, however, not parallel, +but rather follows the curvature of the macroscopic anatomy. +The longitudinal profile of the muscle is generated using a trignometric +function, as opposed to being extracted from medical images. +The benefit to doing so is that the geometry can be (parametrically) created +in `deal.II` itself and the associated microstructural orientation can be +directly linked to the user-defined geometry. + +## Requirements +* Version `8.5` or greater of `deal.II` + +There are no other requirements with regards to the third-party packages that +`deal.II` can link to. + + +## Compiling and running +Similar to the example programs, run +``` +cmake -DDEAL_II_DIR=/path/to/deal.II . +``` +in this directory to configure the problem. +You can switch between debug and release mode by calling either +``` +make debug +``` +or +``` +make release +``` +The problem may then be run with +``` +make run +``` + +Some simulation parameters may be changed by adjusting the `parameters.prm` +file. +Notably, its possible to switch between the model of the biceps and the +reduced geometry used to reproduce the linearised counterpart of the isotonic +contraction numerical experiments conducted by Martins. + + +## Recommended Literature +* Kajee, Y. and Pelteret, J-P. V. and Reddy, B. D. (2013), +The biomechanics of the human tongue. +International Journal for Numerical Methods in Biomedical Engineering +29 , 4, 492-514. +DOI: [10.1002/cnm.2531](http://doi.org/10.1002/cnm.2531); + +* J-P. V. Pelteret, A computational neuromuscular model of the human upper airway with application to the study of obstructive sleep apnoea. PhD Thesis, University of Cape Town, 2013. [http://hdl.handle.net/11427/9519](http://hdl.handle.net/11427/9519); + +* Martins, J. A. C. and Pires, E. B. and Salvado, R. and Dinis, P. B. (1998), +A numerical model of passive and active behaviour of skeletal muscles. +Computer Methods in Applied Mechanics and Engineering +151 , 419-433. +DOI: [10.1016/S0045-7825(97)00162-X](http://doi.org/10.1016/S0045-7825(97)00162-X); + +* Martins, J. A. C. and Pato, M. P. M. and Pires, E. B. (2006), +A finite element model of skeletal muscles. Virtual and Physical Prototyping +1 , 159-170. +DOI: [10.1080/17452750601040626](http://doi.org/10.1080/17452750601040626); + +* Pandy, M. G. and Zajac, F. E. and Sim, E. and Levine, W. S. (1990), +An optimal control model for maximum-height human jumping. +Journal of Biomechanics +23 , 1185-1198. +DOI: [10.1016/0021-9290(90)90376-E](http://doi.org/10.1016/0021-9290(90)90376-E); + +* T.J.R. Hughes (2000), +The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover. +ISBN: [978-0486411811](http://store.doverpublications.com/0486411818.html) + + +## Results +The displacement of the central point on the insertion surface (i.e. the +traction boundary) is plotted against the simulation time when the muscle +is made to undergo concentric contraction. +Within the first second, when the muscle is completely passive, the +displacement increases linearly due to the applied pressure that ramps to a +maximum over this duration. +This response is not entirely unsurprising for this geometrically symmetric, +linear elastic body. +When the muscle is activated, it shortens considerably until during the `1s` +for which the neural signal is applied. +The activation level increases exponentially until is saturates near the +`2s` mark. +At this point the neural signal is removed and the muscle starts to relax. +The contractile level decreases exponentially and the muscle is nearly +completely relaxed by the end of the simulation. +![Problem geometry](./doc/results/results.png) + +As a supplement to the above, the following animation shows the concentric +contraction (under the assumption that it experiences no additional +gravitational loading is present). +All of the highlights that are discussed above can be observed in the +gross displacement of the body, as well as the activation level that is +visualised through the depiction of the underlying microstructure directions. +This also shows how the muscle's cross-section influences the shortening +along the length of the muscle. +![Problem geometry](./doc/results/biceps-no_gravity-concentric_contraction.gif) + + +### Influence of gravity +Just for fun, we can repeat the above experiment with a fraction of the full +gravitational loading applied in the transverse direction. +We apply only a fraction of the full load because the muscle is not sufficiently +well constrained and does not see the support of its surrounding anatomy. +The loading condition is thus somewhat unphysical and, due to the lack of +constraint, the application of the full load results in excessive deformation. + +Here we see the fully passive muscle with partial gravitational loading and a +full traction load +![Displaced solution](./doc/results/biceps-gravity-traction-passive.png) +and its counterpart solution when in the active stage. +![Displaced solution](./doc/results/biceps-gravity-traction-active.png) +The asymmetry of the solution is clearly observed, although the length change +that it exhibits curing the concentric contraction cycle remains somewhat +similar to that demonstrated before. diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/author b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/author new file mode 100644 index 0000000..bbe66f6 --- /dev/null +++ b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/author @@ -0,0 +1,2 @@ +Jean-Paul Pelteret +Tim Hamann diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/builds-on b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/builds-on new file mode 100644 index 0000000..850b582 --- /dev/null +++ b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/builds-on @@ -0,0 +1 @@ +step-8 diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/dependencies b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/dependencies new file mode 100644 index 0000000..e69de29 diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/entry-name b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/entry-name new file mode 100644 index 0000000..28b1833 --- /dev/null +++ b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/entry-name @@ -0,0 +1 @@ +Linear Elastic Active Skeletal Muscle Model diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/geometry/biceps-geometry.png b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/geometry/biceps-geometry.png new file mode 100644 index 0000000..1be1a2f Binary files /dev/null and b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/geometry/biceps-geometry.png differ diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/geometry/geometry.png b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/geometry/geometry.png new file mode 100644 index 0000000..9014dad Binary files /dev/null and b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/geometry/geometry.png differ diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/geometry/geometry.tex b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/geometry/geometry.tex new file mode 100644 index 0000000..4bb670f --- /dev/null +++ b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/geometry/geometry.tex @@ -0,0 +1,24 @@ +\documentclass{standalone} +\usepackage{amsmath} +\usepackage{tikz} + +\begin{document} +\begin{tikzpicture} +\node (P) at (0,0) {\includegraphics[width=10cm]{./biceps-geometry.png}}; +\begin{scope}[scale=0.075] + \foreach \y in {-5,-2.5,...,5} + { + \draw[->, red, thick] (+67,\y) -- (+75,\y) {}; + \draw[-, black, very thick] (-67,\y) -- (-71,\y-1.5) {}; + } + \draw[-, black, very thick] (-67.25,-5) -- (-67.25,+5) {}; +\end{scope} + +\begin{scope}[xshift=-160] + \node[rotate=90] at (0,0) {$\overline{\boldsymbol{\varphi}} \left( t \right)$}; +\end{scope} +\begin{scope}[xshift=170] + \node[rotate=-90] at (0,0) {$\mathbf{\overline{t}\phantom{}^{\textnormal{mech}}} \left( t \right)$}; +\end{scope} +\end{tikzpicture} +\end{document} \ No newline at end of file diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/biceps-gravity-traction-active.png b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/biceps-gravity-traction-active.png new file mode 100644 index 0000000..57af5f9 Binary files /dev/null and b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/biceps-gravity-traction-active.png differ diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/biceps-gravity-traction-passive.png b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/biceps-gravity-traction-passive.png new file mode 100644 index 0000000..d587d42 Binary files /dev/null and b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/biceps-gravity-traction-passive.png differ diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/biceps-no_gravity-concentric_contraction.gif b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/biceps-no_gravity-concentric_contraction.gif new file mode 100644 index 0000000..f39a632 Binary files /dev/null and b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/biceps-no_gravity-concentric_contraction.gif differ diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/displacement_POI.csv b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/displacement_POI.csv new file mode 100644 index 0000000..2bd12eb --- /dev/null +++ b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/displacement_POI.csv @@ -0,0 +1,122 @@ +Time [s],X-displacement [mm] +0,0 +0.025,1.75858 +0.05,3.51717 +0.075,5.27575 +0.1,7.03433 +0.125,8.79292 +0.15,10.5515 +0.175,12.3101 +0.2,14.0687 +0.225,15.8273 +0.25,17.5858 +0.275,19.3444 +0.3,21.103 +0.325,22.8616 +0.35,24.6202 +0.375,26.3788 +0.4,28.1373 +0.425,29.8959 +0.45,31.6545 +0.475,33.4131 +0.5,35.1717 +0.525,36.9303 +0.55,38.6888 +0.575,40.4474 +0.6,42.206 +0.625,43.9646 +0.65,45.7232 +0.675,47.4818 +0.7,49.2403 +0.725,50.9989 +0.75,52.7575 +0.775,54.5161 +0.8,56.2747 +0.825,58.0333 +0.85,59.7918 +0.875,61.5504 +0.9,63.309 +0.925,65.0676 +0.95,66.8262 +0.975,68.5848 +1,60.6516 +1.025,56.4862 +1.05,49.2317 +1.075,44.2888 +1.1,37.6052 +1.125,32.6229 +1.15,26.5684 +1.175,21.7924 +1.2,16.3313 +1.225,11.9289 +1.25,7.00336 +1.275,2.98047 +1.3,-1.38489 +1.325,-5.03264 +1.35,-8.92158 +1.375,-12.1541 +1.4,-15.5626 +1.425,-18.4495 +1.45,-21.4398 +1.475,-23.9507 +1.5,-26.5293 +1.525,-28.7202 +1.55,-30.9316 +1.575,-32.8176 +1.6,-34.7054 +1.625,-36.3108 +1.65,-37.9368 +1.675,-39.3237 +1.7,-40.6995 +1.725,-41.8727 +1.75,-43.0522 +1.775,-44.0535 +1.8,-45.0427 +1.825,-45.8216 +1.85,-46.5598 +1.875,-47.1755 +1.9,-47.6597 +1.925,-48.1319 +1.95,-48.5589 +1.975,-48.6884 +2,-49.2806 +2.025,-44.8251 +2.05,-40.9043 +2.075,-32.6727 +2.1,-25.0115 +2.125,-14.3161 +2.15,-4.50006 +2.175,6.6845 +2.2,16.6699 +2.225,26.5275 +2.25,35.0387 +2.275,42.6003 +2.3,48.64 +2.325,53.4443 +2.35,56.7957 +2.375,59.344 +2.4,61.2859 +2.425,62.8183 +2.45,64.041 +2.475,65.0333 +2.5,65.8564 +2.525,66.5414 +2.55,67.1135 +2.575,67.5952 +2.6,68.0023 +2.625,68.3474 +2.65,68.6404 +2.675,68.8891 +2.7,69.1002 +2.725,69.2803 +2.75,69.434 +2.775,69.5652 +2.8,69.6776 +2.825,69.7734 +2.85,69.8554 +2.875,69.9254 +2.9,69.9854 +2.925,70.0368 +2.95,70.0807 +2.975,70.1184 +3,70.1506 diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/results.png b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/results.png new file mode 100644 index 0000000..e029b7b Binary files /dev/null and b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/results.png differ diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/results.tex b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/results.tex new file mode 100644 index 0000000..053aede --- /dev/null +++ b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/results/results.tex @@ -0,0 +1,32 @@ +\documentclass{standalone} +\usepackage{tikz} +\usepackage{pgfplots} + +\begin{document} +\begin{tikzpicture} + \pgfplotstableread[col sep=comma]{displacement_POI.csv}\loadedtable; + \begin{axis}[ + width = 10cm, + height = 6cm, + xlabel={Time [s]}, + ylabel={X-displacement [mm]}, +% yticklabel style={ +% /pgf/number format/fixed, +% /pgf/number format/precision=2 +% }, +% axis x line* = bottom, +% axis y line* = left, + xmin = 0, + xmax = 3, + ymin = -75, + ymax = +75, + ymajorgrids=true, + ytick distance=25, +% clip=false, + legend style={at={(0.84,+0.175)},anchor=north,font=\footnotesize} + ] + \addplot[red, very thick] table [x={Time [s]}, y={X-displacement [mm]}] {\loadedtable}; + \addlegendentry{$\mathbf{u}^{\textnormal{insertion}}_{x}$}; + \end{axis} +\end{tikzpicture} +\end{document} \ No newline at end of file diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/theory/muscle_fibre-hill_model.png b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/theory/muscle_fibre-hill_model.png new file mode 100644 index 0000000..316dbef Binary files /dev/null and b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/theory/muscle_fibre-hill_model.png differ diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/theory/theory-linear_elastic_active_muscle_model.pdf b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/theory/theory-linear_elastic_active_muscle_model.pdf new file mode 100644 index 0000000..f5adb01 Binary files /dev/null and b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/theory/theory-linear_elastic_active_muscle_model.pdf differ diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/theory/theory-linear_elastic_active_muscle_model.tex b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/theory/theory-linear_elastic_active_muscle_model.tex new file mode 100644 index 0000000..a0c0c76 --- /dev/null +++ b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/theory/theory-linear_elastic_active_muscle_model.tex @@ -0,0 +1,354 @@ +\documentclass[]{scrartcl} +\usepackage{amsmath} +\usepackage{amsfonts} +\usepackage{amssymb} +\usepackage{graphicx} +\usepackage{hyperref} +\usepackage{cleveref} +\usepackage{filecontents} +\usepackage[sort]{natbib} + +\title{Theory: Linear elastic active muscle model} +\author{Jean-Paul Pelteret} + +\begin{filecontents}{\jobname.bib} +@Article{Kajee2013a, + author = {Kajee, Y. and Pelteret, J-P. V. and Reddy, B. D.}, + title = {The biomechanics of the human tongue}, + journal = {International Journal for Numerical Methods in Biomedical Engineering}, + year = {2013}, + volume = {29}, + number = {4}, + pages = {492--514}, + month = {April}, + doi = {10.1002/cnm.2531}, + keywords = {obstructive sleep apnoea ; human tongue ; hill model ; FEM} +} +@Article{Martins1998a, + Title = {A numerical model of passive and active behaviour of skeletal muscles}, + Author = {Martins, J. A. C. and Pires, E. B. and Salvado, R. and Dinis, P. B.}, + Journal = {Computer Methods in Applied Mechanics and Engineering}, + Year = {1998}, + Pages = {419--433}, + Volume = {151}, + Doi = {10.1016/S0045-7825(97)00162-X} +} +@Article{Martins2006a, + Title = {A finite element model of skeletal muscles}, + Author = {Martins, J. A. C. and Pato, M. P. M. and Pires, E. B.}, + Journal = {Virtual and Physical Prototyping}, + Year = {2006}, + Pages = {159--170}, + Volume = {1}, + Doi = {10.1080/17452750601040626} +} +@Article{Pandy1990a, + Title = {An optimal control model for maximum-height human jumping}, + Author = {Pandy, M. G. and Zajac, F. E. and Sim, E. and Levine, W. S.}, + Journal = {Journal of Biomechanics}, + Year = {1990}, + Pages = {1185--1198}, + Volume = {23}, + Doi = {10.1016/0021-9290(90)90376-E} +} +\end{filecontents} + +\begin{document} + +\maketitle + +\begin{abstract} +An introduction to the theory applied to the linear elastic active muscle model of the biceps brachii. +\end{abstract} + +\section{Governing equations for quasi-static linear elasticity} +The strong statement of the balance of linear momentum reads +\begin{gather} +\nabla \cdot \boldsymbol{\sigma} + \mathbf{b} + = \mathbf{0} +\quad \text{on} \quad \Omega \quad , +\end{gather} +where $\nabla = \frac{\partial}{\partial x}$ is a differential operator, +$\boldsymbol{\sigma}$ is the Cauchy stress tensor and +$\mathbf{b} = \rho \mathbf{g}$ is the body force density vector. +This is expressed in index notation as +\begin{gather} +\frac{\partial \sigma_{ij}}{\partial x_{j}} + b_{i} + = 0 +\quad \text{on} \quad \Omega \quad . +\end{gather} +Pre-multiplying the above by test function $\delta \mathbf{v}$ and integrating over the domain $\Omega$ renders +\begin{gather} +- \int\limits_{\Omega} \delta v_{i} \, \frac{\partial \sigma_{ij}}{\partial x_{j}} \, dv + = \int\limits_{\Omega} \delta v_{i} \, b_{i} \, dv +\end{gather} +that, by using the product rule for derivatives (i.e. integration by parts), becomes +\begin{gather} +\int\limits_{\Omega} \frac{\partial \delta v_{i}}{\partial x_{j}} \, \sigma_{ij} \, dv +- \int\limits_{\Omega} \frac{\partial}{\partial x_{j}} \left[ \delta v_{i} \, \sigma_{ij} \right] \, dv + = \int\limits_{\Omega} \delta v_{i} \, b_{i} \, dv +\quad . +\end{gather} +Finally, by applying divergence theorem to the second term in the above, we attain the weak form of the balance of linear momentum +\begin{gather} +\int\limits_{\Omega} \frac{\partial \delta v_{i}}{\partial x_{j}} \, \sigma_{ij} \, dv + = \int\limits_{\Omega} \delta v_{i} \, b_{i} \, dv + + \int\limits_{\partial\Omega} \delta v_{i} \, \underbrace{\sigma_{ij} \, n_{j}}_{\bar{t}_{i}} \, da +\quad , +\label{equ: Weak form of the balance of linear momentum} +\end{gather} +wherein $\mathbf{n}$ represents the outward facing normal on $\partial\Omega$, the boundary of the domain, and $\bar{\mathbf{t}}$ the prescribed traction on the Neumann boundary. + +\section{Constitutive law: A linearised Hill three-element active muscle model with surrounding matrix \citep{Kajee2013a}} + +The linear constitutive law used to model active muscle tissue is derived by \citep{Kajee2013a} from the nonlinear model developed by \citep{Martins1998a,Martins2006a}. +In the representation given here, we deviate slightly from the notation given in \citep{Kajee2013a} to facilitate its implementation. + +\subsubsection*{Embedding of one-dimensional fibre model into three-dimensional space} + +We begin by defining the decomposition of the Cauchy stress tensor into a matrix and fibre contribution as +\begin{gather} +\boldsymbol{\sigma} + = \boldsymbol{\sigma}_{m} + \boldsymbol{\sigma}_{f} +\end{gather} +where $m,f$ respectively denote contributions from the surrounding matrix and muscle fibres. +The isotropic linear constitutive law for the matrix surrounding the muscle fibres is +\begin{gather} +\boldsymbol{\sigma}_{m} + = \boldsymbol{\mathbb{C}}_{m} : \boldsymbol{\varepsilon} +\end{gather} +where +$\boldsymbol{\mathbb{C}}_{m}$ is the stiffness tensor for the matrix, and +the small strain tensor +\begin{gather} +\boldsymbol{\varepsilon} + = \frac{1}{2} \left[ \nabla \mathbf{u} + \left[ \nabla \mathbf{u} \right]^{T} \right] +\quad . +\end{gather} +The fibre stress and strain are computed by +\begin{gather} +\boldsymbol{\sigma}_{f} + = T_{f} \mathbf{m} \otimes \mathbf{m} +\quad , \quad +\varepsilon_{f} + = \left[ \mathbf{m} \otimes \mathbf{m} \right] : \boldsymbol{\varepsilon} +\end{gather} + +\subsubsection*{Linearised version of Martin's one-dimensional muscle model} + +\Cref{fig: Hill three-element model} shows an analogue for the sarcomere, the smallest building-block of active muscle fibres. +\begin{figure}[!ht] +\centering +\includegraphics[height=0.23\textheight]{./muscle_fibre-hill_model.png} +\caption{Schematic of the Hill-type muscle fibre \citep{Kajee2013a}. +\label{fig: Hill three-element model} +} +\end{figure} +The distributions of strains and stresses within the various elements of the representative model is determined by their arrangement with respect to one another. +In the linearised (small-strain) version of the Hill three-element model, the decomposition of stress in the fibre as a whole and the one parallel branch are +\begin{gather} +T_{f} + = T_{p} + T_{s} +\quad\textnormal{and}\quad +T_{c} + = T_{s} +\end{gather} +where $T$ is a measure of nominal stress, +and the subscripts $f,p,s,c$ respectively denote the fibre (as a whole), and the parallel, series and contractile element in the Hill model. +Similarly, the decomposition of the (small) strains in the Hill model are +\begin{gather} +\varepsilon_{f} + = \varepsilon_{p} + \equiv \varepsilon_{s} + \varepsilon_{c} +\quad . +\end{gather} + +The constitutive laws governing the response of each element are as follows: +\begin{gather} +T_{p} + = T_{0} \, f_{p} +\quad , \quad +T_{s} + = T_{0} \, f_{s} +\quad \textnormal{and} \quad +T_{c} + = f_{c}^{l} \left( \varepsilon_{c} \right) \, f_{c}^{v} \left( \dot{\varepsilon}_{c} \right) \, \alpha\left( u \left( t \right) \right) +\quad . +\end{gather} +where $T_{0}$ is the nominal stress, a physiological constant which defines to the maximum force of contraction under isometric conditions. +Here the driver functions for the passive parallel and series elements are +\begin{gather} +f_{p} \left( \varepsilon_{f} \right) + = \begin{cases} + m_{p} \varepsilon_{f} \quad &\textnormal{if} \quad \varepsilon_{f} > 0 \\ + 0 \quad &\textnormal{otherwise} + \end{cases} +\\ +f_{s} \left( \varepsilon_{s} \right) + = \begin{cases} + m_{s} \varepsilon_{s} \equiv m_{s} \left[ \varepsilon_{f} - \varepsilon_{c} \right] \quad &\textnormal{if} \quad \varepsilon_{s} \equiv \varepsilon_{f} - \varepsilon_{c} > 0 \\ + 0 \quad &\textnormal{otherwise} + \end{cases} +\quad . +\end{gather} +Here the strain relationship between the elements is used to remove $\varepsilon_{s}$ as an unknown. +For the active contractile element, the force-length and force-velocity relationships are approximated as +\begin{gather} +f_{c}^{l} \left( \varepsilon_{c} \right) + = \begin{cases} + 1 \quad &\textnormal{if} \quad -0.5 \leq \varepsilon_{c} \leq 0.5 \\ + 0 \quad &\textnormal{otherwise} + \end{cases} +\\ +f_{c}^{v} \left( \dot{\varepsilon}_{c} \right) + = \begin{cases} + 0 &\textnormal{if} \quad \dot{\varepsilon}_{c} < -5 \\ + \frac{1}{5} \dot{\varepsilon}_{c} + 1 &\textnormal{if} \quad -5 \leq \dot{\varepsilon}_{c} < 3 \\ + 1.6 \quad &\textnormal{otherwise} + \end{cases} +\quad , +\end{gather} +the latter of which we can write in general as +\begin{gather} +f_{c}^{v} \left( \dot{\varepsilon}_{c} \right) + = m_{c}^{v} \dot{\varepsilon}_{c} + c_{c}^{v} +\quad . +\end{gather} +Note that alternative linearisations for these terms are possible, and that the rate-dependence of the contractile element makes this model ``visco-elastic''. +The differential equation that defines the muscle activation model \citep{Pandy1990a} is expressed a function of the neural signal $u \left( t \right)$ by +\begin{gather} +\dot{\alpha}\left( u \left( t \right) \right) + =\frac{1}{\tau_{r}} \left[ 1 - \alpha \right] u + \frac{1}{\tau_{f}} \left[ \alpha_{\min} - \alpha \right] \left[ 1- u \right] +\quad . +\end{gather} +The parameters $\tau_{r}$ and $\tau_{f}$ control the rise and fall of the activation function with respect to the history of the neural signal, and $\alpha_{\min}$ is the minimum activation level (real muscles are never completely inactive; they always retain some degree of tetanisation). + +\subsubsection*{Time differentiation} +For all time derivatives we employ a first-order backward Euler scheme. +Therefore the contractile strain rate and rate of change of muscle activation at timestep $n$ are approximated as +\begin{gather} +\dot{\varepsilon}_{c} + \approx \frac{\varepsilon_{c}^{n} - \varepsilon_{c}^{n-1}}{\Delta t} +\\ +\dot{\alpha} + \approx \frac{{\alpha}^{n} - \alpha^{n-1}}{\Delta t} +\quad . +\end{gather} +Consequently the expression for the force-velocity relationship and activation level can be explicitly stated in terms of the history variables $\varepsilon_{c}^{n-1},\alpha^{n-1}$ and the remaining unknowns $\varepsilon_{c}^{n},\alpha^{n}$. + +\subsubsection*{Substitution of fibre constitutive laws into one-dimensional stress relationship} +From the equivalence of $T_{c}$ and $T_{s}$, substituting in all of the salient previously derived expressions and considering $\alpha > 0$, we can extract the explicit expression for $\varepsilon_{c}$ in terms of $\varepsilon_{f}$ by the following steps: +\begin{gather*} +f_{c}^{l} \, f_{c}^{v} \, \alpha = f_{s} \\ +%f_{c}^{l} \, \left[ m_{c}^{v} \dot{\varepsilon}_{c} + c_{c}^{v} \right] \, \alpha = m_{s} \left[ \varepsilon_{f} - \varepsilon_{c} \right] \\ +\Rightarrow \quad f_{c}^{l} \, \left[ m_{c}^{v} \frac{\varepsilon_{c} - \varepsilon_{c}^{n-1}}{\Delta t} + c_{c}^{v} \right] \, \alpha = m_{s} \left[ \varepsilon_{f} - \varepsilon_{c} \right] +\end{gather*} +that, with some further rearrangement, becomes +\begin{align} +\varepsilon_{c} + &= \underbrace{\left[ f_{c}^{l} m_{c}^{v} \frac{1}{\Delta t} \alpha + m_{s} \right]}_{\beta}\phantom{}^{-1} \left[ m_{s} \varepsilon_{f} + \underbrace{f_{c}^{l} \alpha \left[ m_{c}^{v} \varepsilon_{c}^{n-1} \frac{1}{\Delta t} - c_{c}^{v} \right]}_{\gamma}\right] \notag\\ + &= \frac{m_{s}}{\beta} \varepsilon_{f} + \frac{\gamma}{\beta} +\end{align} +Note that $\beta > 0$ under all conditions as $m_{s} > 0$ during contraction. + +\subsubsection*{Substitution of constitutive laws into three-dimensional stress relationship} +For the most general case, we can decompose the total Cauchy stress as +\begin{align} +\boldsymbol{\sigma} + &= \boldsymbol{\mathbb{C}}_{m} : \boldsymbol{\varepsilon} + \boldsymbol{\sigma}_{f} \notag\\ + &= \boldsymbol{\mathbb{C}}_{m} : \boldsymbol{\varepsilon} + T_{f} \mathbf{m} \otimes \mathbf{m} \notag\\ + &= \boldsymbol{\mathbb{C}}_{m} : \boldsymbol{\varepsilon} + T_{0} \left[ m_{p} \varepsilon_{f} + m_{s} \left[ \varepsilon_{f} - \varepsilon_{c} \right] \right] \mathbf{m} \otimes \mathbf{m} \notag\\ + &= \boldsymbol{\mathbb{C}}_{m} : \boldsymbol{\varepsilon} + T_{0} \left[ m_{p} \varepsilon_{f} + m_{s} \left[ \varepsilon_{f} - \left[ \frac{m_{s}}{\beta} \varepsilon_{f} + \frac{\gamma}{\beta} \right] \right] \right] \mathbf{m} \otimes \mathbf{m} \notag\\ + &= \boldsymbol{\mathbb{C}}_{m} : \boldsymbol{\varepsilon} + T_{0} \left[ m_{p} + m_{s} - \frac{m_{s}^{2}}{\beta} \right] \varepsilon_{f} \mathbf{m} \otimes \mathbf{m} - \left[ T_{0} m_{s} \frac{\gamma}{\beta} \right] \mathbf{m} \otimes \mathbf{m} \notag\\ + &= \left[ \boldsymbol{\mathbb{C}}_{m} + \underbrace{T_{0} \left[ m_{p} + m_{s} - \frac{m_{s}^{2}}{\beta} \right] \mathbf{m} \otimes \mathbf{m} \otimes \mathbf{m} \otimes \mathbf{m} }_{\boldsymbol{\mathbb{C}}_{f}^{\ast}} \right] : \boldsymbol{\varepsilon} \left( \mathbf{u} \right) - \underbrace{\left[ T_{0} m_{s} \frac{\gamma}{\beta} \right] \mathbf{m} \otimes \mathbf{m}}_{\boldsymbol{\sigma}_{f}^{\ast}} +\quad . +\label{equ: Full expansion of constitutive laws} +\end{align} +Note here that the first term on the right hand side ($\left[ \boldsymbol{\mathbb{C}}_{m} + \boldsymbol{\mathbb{C}}_{f}^{\ast} \right] : \boldsymbol{\varepsilon} \left( \mathbf{u} \right)$) is dependent on the solution, and the second term ($\boldsymbol{\sigma}_{f}^{\ast}$) depends only on local history variables. + +\section{Finite element discretisation} +Combining \cref{equ: Weak form of the balance of linear momentum,equ: Full expansion of constitutive laws} renders the complete expression of the balance of linear momentum, with accommodation of the muscle fibre model, namely +\begin{gather} +\int\limits_{\Omega} \frac{\partial \delta v_{i}}{\partial x_{j}} \, \left[ \boldsymbol{\mathbb{C}}_{m} + \boldsymbol{\mathbb{C}}_{f}^{\ast} \right]_{ijkl} \varepsilon_{kl} \, dv + = \int\limits_{\Omega} \delta v_{i} \, b_{i} \, dv + + \int\limits_{\partial\Omega} \delta v_{i} \, \underbrace{\sigma_{ij} \, n_{j}}_{\bar{t}_{i}} \, da + - \int\limits_{\Omega} \frac{\partial \delta v_{i}}{\partial x_{j}} \, \left[ \boldsymbol{\sigma}_{f}^{\ast} \right]_{ij} \, dv +\quad , +\label{equ: Weak form of the balance of linear momentum: Muscle model} +\end{gather} + +We discretise the trial solution and test function using finite element shape functions (ansatz) +\begin{gather} +\mathbf{u} \left( \mathbf{x} \right) + \approx \sum\limits_{I} \boldsymbol{\varPhi}^{I} \left( \mathbf{x} \right) u^{I} +\quad , \quad +\mathbf{v} \left( \mathbf{x} \right) + \approx \sum\limits_{I} \boldsymbol{\varPhi}^{I} \left( \mathbf{x} \right) v^{I} +\end{gather} +where $\mathbf{N}^{I} \left( \mathbf{x} \right)$ is the (position-dependent) vector-valued finite element shape function corresponding to the $I^{\textnormal{th}}$ degree-of-freedom, and $u^{I}, v^{I}$ are coefficients of the solution and trial function. +In \texttt{deal.II} nomenclature, the shape function is computed from a scalar base shape function and some expansion into higher-dimensional space by +\begin{gather} +\boldsymbol{\varPhi}^{I} \left( \mathbf{x} \right) + = N^{I} \left( \mathbf{x} \right) \mathbf{e}_{\textnormal{comp}(I)} +\end{gather} +where $N^{I}$ is a scalar shape function and $\mathbf{e}_{\textnormal{comp}(I)}$ is the basis direction associated with the $I^{\textnormal{th}}$ degree-of-freedom. +Therefore, the $j^{\textnormal{th}}$ local component of shape function $\boldsymbol{\varPhi}^{I} \left( \mathbf{x} \right)$ is given by +\begin{gather} +\left[\boldsymbol{\varPhi}^{I} \left( \mathbf{x} \right)\right]_{j} + = N^{I} \left( \mathbf{x} \right) \left[ \mathbf{e}_{\textnormal{comp}(I)} \right]_{j} + = N^{I} \left( \mathbf{x} \right) \delta_{\textnormal{comp}(I) j} +\quad . +\end{gather} +where $\delta_{ij}$ is the Kronecker delta. +Note that in this instance we use the same ansatz for the test and trial spaces, and the $0 \leq \textnormal{comp}(I), j < \textnormal{spacedim}$. + +We now use these shape functions to discretise the weak expression for the balance of linear momentum. +Starting on the right-hand side of \cref{equ: Weak form of the balance of linear momentum: Muscle model}, the body force and traction contributions are computed by +\begin{align} +\int\limits_{\Omega} \delta v_{i} \, b_{i} \, dv + &= \int\limits_{\Omega} \left[ \sum\limits_{I} \boldsymbol{\varPhi}^{I} \left( \mathbf{x} \right) \delta v^{I} \right]_{i} \, b_{i} \, dv + = \sum\limits_{I} \delta v^{I} \int\limits_{\Omega} \left[ \boldsymbol{\varPhi}^{I} \left( \mathbf{x} \right) \right]_{i} \, b_{i} \, dv \notag\\ + &= \sum\limits_{I} \delta v^{I} \int\limits_{\Omega} N^{I} \left( \mathbf{x} \right) \delta_{\textnormal{comp}(I) i} \, b_{i} \, dv + = \sum\limits_{I} \delta v^{I} \int\limits_{\Omega} N^{I} \, b_{\textnormal{comp}(I)} \, dv +\label{equ: Discretisation: Body force} +\\ +\int\limits_{\Omega} \delta v_{i} \, t_{i} \, dv + &= \sum\limits_{I} \delta v^{I} \int\limits_{\Omega} N^{I} \, t_{\textnormal{comp}(I)} \, dv +\label{equ: Discretisation: Traction} +\quad . +\end{align} +while the contribution to the right-hand side that arise from the history variables is +\begin{align} +- \int\limits_{\Omega} \frac{\partial}{\partial x_{j}} \left[ \delta v_{i} \right] \, \left[ \boldsymbol{\sigma}_{f}^{\ast} \right]_{ij} \, dv + &= - \int\limits_{\Omega} \frac{\partial}{\partial x_{j}} \left[ \sum\limits_{I} \boldsymbol{\varPhi}^{I} \left( \mathbf{x} \right) \delta v^{I} \right]_{i} \, \left[ \boldsymbol{\sigma}_{f}^{\ast} \right]_{ij} \, dv \notag\\ + &= - \sum\limits_{I} \delta v^{I} \int\limits_{\Omega} \frac{\partial}{\partial x_{j}} \left[ \boldsymbol{\varPhi}^{I} \left( \mathbf{x} \right) \right]_{i} \, \left[ \boldsymbol{\sigma}_{f}^{\ast} \right]_{ij} \, dv \notag\\ + &= - \sum\limits_{I} \delta v^{I} \int\limits_{\Omega} \frac{\partial}{\partial x_{j}} \left[ N^{I} \left( \mathbf{x} \right) \delta_{\textnormal{comp}(I) i} \right] \, \left[ \boldsymbol{\sigma}_{f}^{\ast} \right]_{ij} \, dv \notag\\ + &= - \sum\limits_{I} \delta v^{I} \int\limits_{\Omega} \frac{\partial N^{I} \left( \mathbf{x} \right)}{\partial x_{j}} \, \left[ \boldsymbol{\sigma}_{f}^{\ast} \right]_{\textnormal{comp}(I)j} \, dv +\label{equ: Discretisation: Stress history} +\quad . +\end{align} +The last component of \cref{equ: Weak form of the balance of linear momentum: Muscle model} that we wish to express in discrete form is the left-hand side of the equation. +Before we do, we observe that using the minor symmetry of the material stiffness tensor we can re-express the contraction of it and the small strain tensor as +\begin{gather} +\boldsymbol{\mathbb{C}} : \boldsymbol{\varepsilon} + = \boldsymbol{\mathbb{C}} : \frac{1}{2} \left[ \nabla \mathbf{u} + \left[ \nabla \mathbf{u} \right]^{T} \right] + \equiv \boldsymbol{\mathbb{C}} : \nabla \mathbf{u} +\end{gather} +Therefore, this contribution written in discrete form is +\begin{align} +&\int\limits_{\Omega} \frac{\partial \delta v_{i}}{\partial x_{j}} \, \left[ \boldsymbol{\mathbb{C}}_{m} + \boldsymbol{\mathbb{C}}_{f}^{\ast} \right]_{ijkl} \varepsilon_{kl} \, dv + \equiv \int\limits_{\Omega} \frac{\partial \delta v_{i}}{\partial x_{j}} \, \mathbb{C}_{ijkl} \, \frac{\partial \delta u_{k}}{\partial x_{l}} \, dv \notag\\ + &\equiv \int\limits_{\Omega} \frac{\partial}{\partial x_{j}} \left[ \sum\limits_{I} \boldsymbol{\varPhi}^{I} \left( \mathbf{x} \right) \delta v^{I} \right]_{i} \, \mathbb{C}_{ijkl} \, \frac{\partial}{\partial x_{l}} \left[ \sum\limits_{J} \boldsymbol{\varPhi}^{J} \left( \mathbf{x} \right) \delta u^{J} \right]_{k} \, dv \notag\\ + &\equiv \sum\limits_{I,J} \delta v^{I} \left[ \int\limits_{\Omega} \frac{\partial}{\partial x_{j}} \left[ \boldsymbol{\varPhi}^{I} \left( \mathbf{x} \right) \right]_{i} \, \mathbb{C}_{ijkl} \, \frac{\partial}{\partial x_{l}} \left[ \boldsymbol{\varPhi}^{J} \left( \mathbf{x} \right) \right]_{k} \, dv \right] \delta u^{J} \notag\\ + &\equiv \sum\limits_{I,J} \delta v^{I} \left[ \int\limits_{\Omega} \frac{\partial N^{I} \left( \mathbf{x} \right)}{\partial x_{j}} \delta_{\textnormal{comp}(I) i} \, \mathbb{C}_{ijkl} \, \frac{\partial N^{J} \left( \mathbf{x} \right)}{\partial x_{l}} \delta_{\textnormal{comp}(J) k} \, dv \right] \delta u^{J} \notag\\ + &\equiv \sum\limits_{I,J} \delta v^{I} \left[ \int\limits_{\Omega} \frac{\partial N^{I} \left( \mathbf{x} \right)}{\partial x_{j}} \, \mathbb{C}_{\textnormal{comp}(I) \,j \, \textnormal{comp}(J) \, l} \, \frac{\partial N^{J} \left( \mathbf{x} \right)}{\partial x_{l}} \, dv \right] \delta u^{J} +\label{equ: Discretisation: Material tangent} +\quad . +\end{align} + +\Cref{equ: Discretisation: Body force,equ: Discretisation: Traction,equ: Discretisation: Stress history,equ: Discretisation: Material tangent} are collectively used to develop the system of linear equations that are solved at each time step. + +\bibliographystyle{plain} +\bibliography{\jobname} + +\end{document} diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/tooltip b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/tooltip new file mode 100644 index 0000000..a70a9e0 --- /dev/null +++ b/Linear_Elastic_Active_Skeletal_Muscle_Model/doc/tooltip @@ -0,0 +1 @@ +Linear elasticity including active skeletal muscle model solving the concentric contraction of an idealized biceps brachii. diff --git a/Linear_Elastic_Active_Skeletal_Muscle_Model/parameters.prm b/Linear_Elastic_Active_Skeletal_Muscle_Model/parameters.prm new file mode 100644 index 0000000..a938a1e --- /dev/null +++ b/Linear_Elastic_Active_Skeletal_Muscle_Model/parameters.prm @@ -0,0 +1,66 @@ +# Listing of Parameters +# --------------------- +subsection Finite element system + # Displacement system polynomial order + set Polynomial degree = 1 + + # Gauss quadrature order + set Quadrature order = 2 +end + + +subsection Problem + # Choose the problem to be solved + # Options: IsotonicContraction ; BicepsBrachii + set Problem = IsotonicContraction +end + + +subsection Biceps Brachii geometry + # Axial length of the muscle + set Axial length = 250 # (in millimetres) + + # Insertion and origin radius + set Radius insertion and origin = 10 # (in millimetres) + + # Radius at the midpoint of the muscle + set Radius midpoint = 40 # (in millimetres) + + # Global grid scaling factor (metres -> millimetres) + set Grid scale = 1e-3 + + # Number of elements along the muscle axis + set Elements along axis = 40 + + # Control the discretisation in the radial direction + set Radial refinements = 3 + + # Choose whether to include gravitational + # effects (in the y-direction; perpendicular + # to the muscle axis) + set Gravity = false + + # Applied distributed axial force (in Newtons) + set Axial force = 20 +end + + +subsection Neurological signal + # Time at which to start muscle activation + set Start time = 1.0 + + # Time at which to remove muscle activation signal + set End time = 2.0 +end + + +subsection Time + # End time + set End time = 3.0 + + # Force ramp end time + set End ramp time = 1.0 + + # Time step size + set Time step size = 0.1 +end