274 * <a name=
"Linear_active_muscle_model.cc-Includefiles"></a>
322 * <a name=
"Linear_active_muscle_model.cc-Runtimeparameters"></a>
336 * <a name=
"Linear_active_muscle_model.cc-FiniteElementsystem"></a>
359 *
prm.enter_subsection(
"Finite element system");
361 *
prm.declare_entry(
"Polynomial degree",
"1",
363 *
"Displacement system polynomial order");
365 *
prm.declare_entry(
"Quadrature order",
"2",
367 *
"Gauss quadrature order");
369 *
prm.leave_subsection();
374 *
prm.enter_subsection(
"Finite element system");
376 *
poly_degree = prm.get_integer(
"Polynomial degree");
377 *
quad_order = prm.get_integer(
"Quadrature order");
379 *
prm.leave_subsection();
385 * <a name=
"Linear_active_muscle_model.cc-Problem"></a>
406 *
prm.enter_subsection(
"Problem");
408 *
prm.declare_entry(
"Problem",
"IsotonicContraction",
410 *
"The problem that is to be solved");
412 *
prm.leave_subsection();
417 *
prm.enter_subsection(
"Problem");
419 *
problem = prm.get(
"Problem");
421 *
prm.leave_subsection();
427 * <a name=
"Linear_active_muscle_model.cc-IsotonicContractionGeometry"></a>
468 * <a name=
"Linear_active_muscle_model.cc-BicepsBrachiiGeometry"></a>
504 *
prm.enter_subsection(
"Biceps Brachii geometry");
506 *
prm.declare_entry(
"Axial length",
"250",
508 *
"Axial length of the muscle");
510 *
prm.declare_entry(
"Radius insertion and origin",
"5",
512 *
"Insertion and origin radius");
514 *
prm.declare_entry(
"Radius midpoint",
"7.5",
516 *
"Radius at the midpoint of the muscle");
518 *
prm.declare_entry(
"Grid scale",
"1e-3",
520 *
"Global grid scaling factor");
522 *
prm.declare_entry(
"Elements along axis",
"32",
524 *
"Number of elements along the muscle axis");
526 *
prm.declare_entry(
"Radial refinements",
"4",
528 *
"Control the discretisation in the radial direction");
530 *
prm.declare_entry(
"Gravity",
"false",
532 *
"Include the effects of gravity (in the y-direction; "
533 *
" perpendicular to the muscle axis)");
535 *
prm.declare_entry(
"Axial force",
"1",
537 *
"Applied distributed axial force (in Newtons)");
539 *
prm.leave_subsection();
544 *
prm.enter_subsection(
"Biceps Brachii geometry");
549 *
scale = prm.get_double(
"Grid scale");
555 *
prm.leave_subsection();
564 * <a name=
"Linear_active_muscle_model.cc-Neurologicalsignal"></a>
585 *
prm.enter_subsection(
"Neurological signal");
587 *
prm.declare_entry(
"Start time",
"1.0",
589 *
"Time at which to start muscle activation");
591 *
prm.declare_entry(
"End time",
"2.0",
593 *
"Time at which to remove muscle activation signal");
595 *
prm.leave_subsection();
600 *
prm.enter_subsection(
"Neurological signal");
605 *
prm.leave_subsection();
608 *
ExcMessage(
"Invalid neural signal times."));
614 * <a name=
"Linear_active_muscle_model.cc-Time"></a>
637 *
prm.enter_subsection(
"Time");
639 *
prm.declare_entry(
"End time",
"3",
643 *
prm.declare_entry(
"End ramp time",
"1",
645 *
"Force ramp end time");
647 *
prm.declare_entry(
"Time step size",
"0.1",
651 *
prm.leave_subsection();
656 *
prm.enter_subsection(
"Time");
658 *
end_time = prm.get_double(
"End time");
660 *
delta_t = prm.get_double(
"Time step size");
662 *
prm.leave_subsection();
668 * <a name=
"Linear_active_muscle_model.cc-Allparameters"></a>
696 *
declare_parameters(prm);
698 *
parse_parameters(prm);
703 *
FESystem::declare_parameters(prm);
704 *
Problem::declare_parameters(prm);
705 *
IsotonicContraction::declare_parameters(prm);
706 *
BicepsBrachii::declare_parameters(prm);
707 *
NeurologicalSignal::declare_parameters(prm);
708 *
Time::declare_parameters(prm);
713 *
FESystem::parse_parameters(prm);
714 *
Problem::parse_parameters(prm);
715 *
IsotonicContraction::parse_parameters(prm);
716 *
BicepsBrachii::parse_parameters(prm);
717 *
NeurologicalSignal::parse_parameters(prm);
718 *
Time::parse_parameters(prm);
726 *
if (
problem ==
"IsotonicContraction")
741 * <a name=
"Linear_active_muscle_model.cc-Bodyforcevalues"></a>
756 *
virtual void vector_value (
const Point<dim> &p,
759 *
virtual void vector_value_list (
const std::vector<
Point<dim> > &points,
777 *
Assert(
M.
norm() == 1.0, ExcMessage(
"Direction vector is not a unit vector"));
787 *
ExcDimensionMismatch (
values.size(), dim));
788 *
Assert (dim >= 2, ExcNotImplemented());
789 *
for (
unsigned int d=0;
d<dim; ++
d)
800 *
Assert (value_list.size() == points.size(),
801 *
ExcDimensionMismatch (value_list.size(), points.size()));
803 *
const unsigned int n_points = points.size();
805 *
for (
unsigned int p=0; p<n_points; ++p)
815 *
const double area);
818 *
virtual void vector_value (
const Point<dim> &p,
821 *
virtual void vector_value_list (
const std::vector<
Point<dim> > &points,
843 *
ExcDimensionMismatch (
values.size(), dim));
844 *
Assert (dim == 3, ExcNotImplemented());
861 *
Assert (value_list.size() == points.size(),
862 *
ExcDimensionMismatch (value_list.size(), points.size()));
864 *
const unsigned int n_points = points.size();
866 *
for (
unsigned int p=0; p<n_points; ++p)
874 * <a name=
"Linear_active_muscle_model.cc-Utilityfunctions"></a>
885 *
Assert (
grad.size() == dim, ExcInternalError());
888 *
for (
unsigned int i=0; i<dim; ++i)
889 *
for (
unsigned int j=0;
j<dim; ++
j)
898 *
Assert (
grad.size() == dim, ExcInternalError());
901 *
for (
unsigned int i=0; i<dim; ++i)
904 *
for (
unsigned int i=0; i<dim; ++i)
905 *
for (
unsigned int j=i+1;
j<dim; ++
j)
913 * <a name=
"Linear_active_muscle_model.cc-Propertiesformusclematrix"></a>
922 *
static const double E;
923 *
static const double nu;
925 *
static const double mu;
926 *
static const double lambda;
929 *
const double MuscleMatrix::E = 26
e3;
937 * <a name=
"Linear_active_muscle_model.cc-Localdataformusclefibres"></a>
977 *
ExcMessage(
"Fibre direction is not a unit vector"));
1033 *
template <
int dim>
1037 *
static const double tau_r = 0.15;
1038 *
static const double tau_f = 0.15;
1048 *
const double c = 1.0/
tau_f;
1051 *
const double p =
b*
u + c;
1052 *
const double q =
f1*
u +
d;
1059 *
template <
int dim>
1063 *
static const double a = 12.43;
1071 *
static const double m_p = 2.0*
A*a/1
e2;
1078 *
template <
int dim>
1088 *
template <
int dim>
1097 *
template <
int dim>
1108 *
template <
int dim>
1119 *
template <
int dim>
1125 *
template <
int dim>
1131 *
template <
int dim>
1163 * <a name=
"Linear_active_muscle_model.cc-ThecodeLinearMuscleModelProblemcodeclasstemplate"></a>
1170 *
template <
int dim>
1189 *
const double time)
const;
1213 *
const double t_end;
1230 *
std::vector< std::vector<MuscleFibre<dim> > >
fibre_data;
1246 * <a name=
"Linear_active_muscle_model.cc-LinearMuscleModelProblemLinearMuscleModelProblem"></a>
1253 *
template <
int dim>
1258 *
fe (
FE_Q<dim>(parameters.poly_degree), dim),
1259 *
qf_cell (parameters.quad_order),
1260 *
qf_face (parameters.quad_order),
1261 *
t_end (parameters.end_time),
1262 *
dt (parameters.delta_t),
1264 *
body_force ((parameters.problem ==
"BicepsBrachii" &¶meters.include_gravity ==
true) ?
1267 *
traction (parameters.problem ==
"BicepsBrachii" ?
1269 *
M_PI*
std::pow(parameters.radius_insertion_origin *parameters.scale,2.0) ) :
1280 * <a name=
"Linear_active_muscle_model.cc-LinearMuscleModelProblemLinearMuscleModelProblem"></a>
1287 *
template <
int dim>
1290 *
dof_handler.
clear ();
1297 * <a name=
"Linear_active_muscle_model.cc-LinearMuscleModelProblemmake_grid"></a>
1329 *
ExcMessage(
"All points must have x-coordinate > 0"));
1368 *
static const double eps = 1
e-6;
1374 *
static const double tol = 1
e-9;
1381 *
dir /= dir.
norm();
1388 *
const double r_mid;
1410 *
ExcMessage(
"All points must have x-coordinate > 0"));
1418 *
template <
int dim>
1423 *
if (parameters.problem ==
"IsotonicContraction")
1426 *
-parameters.half_length_y,
1427 *
-parameters.half_length_z);
1429 *
parameters.half_length_y,
1430 *
parameters.half_length_z);
1436 *
for (; cell !=
endc; ++cell)
1440 *
if (cell->face(face)->at_boundary() ==
true)
1442 *
if (cell->face(face)->center()[0] == -parameters.half_length_x)
1443 *
cell->face(face)->set_boundary_id(parameters.bid_CC_dirichlet_symm_X);
1444 *
else if (cell->face(face)->center()[0] == parameters.half_length_x)
1445 *
cell->face(face)->set_boundary_id(parameters.bid_CC_neumann);
1446 *
else if (
std::abs(cell->face(face)->center()[2]) == parameters.half_length_z)
1447 *
cell->face(face)->set_boundary_id(parameters.bid_CC_dirichlet_symm_Z);
1454 *
else if (parameters.problem ==
"BicepsBrachii")
1460 *
parameters.radius_insertion_origin);
1467 *
if (cell->face(face)->at_boundary() ==
true)
1468 *
cell->face(face)->set_all_manifold_ids(0);
1472 *
tria_cap.refine_global(parameters.n_refinements_radial);
1478 *
parameters.elements_along_axis,
1479 *
parameters.axial_length,
1505 * GridTools::transform (BicepsGeometry<dim>(parameters.axial_length,
1506 * parameters.radius_insertion_origin,
1507 * parameters.radius_midpoint), triangulation);
1514 * typename Triangulation<dim>::active_cell_iterator cell =
1515 * triangulation.begin_active(), endc = triangulation.end();
1516 * for (; cell != endc; ++cell)
1518 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
1520 * if (cell->face(face)->at_boundary() == true)
1522 * static const double tol =1e-6;
1523 * if (std::abs(cell->face(face)->center()[0]) < tol) // -X oriented face
1524 * cell->face(face)->set_boundary_id(parameters.bid_BB_dirichlet_X); // Dirichlet
1525 * else if (std::abs(cell->face(face)->center()[0] - parameters.axial_length) < tol) // +X oriented face
1526 * cell->face(face)->set_boundary_id(parameters.bid_BB_neumann); // Neumann
1533 * Finally resize the grid
1536 * GridTools::scale (parameters.scale, triangulation);
1539 * AssertThrow(false, ExcNotImplemented());
1545 * <a name="Linear_active_muscle_model.cc-LinearMuscleModelProblemsetup_muscle_fibres"></a>
1546 * <h4>LinearMuscleModelProblem::setup_muscle_fibres</h4>
1552 * template <int dim>
1553 * void LinearMuscleModelProblem<dim>::setup_muscle_fibres ()
1555 * fibre_data.clear();
1556 * const unsigned int n_cells = triangulation.n_active_cells();
1557 * fibre_data.resize(n_cells);
1558 * const unsigned int n_q_points_cell = qf_cell.size();
1560 * if (parameters.problem == "IsotonicContraction")
1562 * MuscleFibre<dim> fibre_template (Tensor<1,dim>({1,0,0}));
1564 * for (unsigned int cell_no=0; cell_no<triangulation.n_active_cells(); ++cell_no)
1566 * fibre_data[cell_no].resize(n_q_points_cell);
1567 * for (unsigned int q_point_cell=0; q_point_cell<n_q_points_cell; ++q_point_cell)
1569 * fibre_data[cell_no][q_point_cell] = fibre_template;
1573 * else if (parameters.problem == "BicepsBrachii")
1575 * FEValues<dim> fe_values (fe, qf_cell, update_quadrature_points);
1576 * BicepsGeometry<dim> bicep_geom (parameters.axial_length,
1577 * parameters.radius_insertion_origin,
1578 * parameters.radius_midpoint);
1580 * unsigned int cell_no = 0;
1581 * for (typename Triangulation<dim>::active_cell_iterator
1582 * cell = triangulation.begin_active();
1583 * cell != triangulation.end();
1584 * ++cell, ++cell_no)
1586 * Assert(cell_no<fibre_data.size(), ExcMessage("Trying to access fibre data not stored for this cell index"));
1587 * fe_values.reinit(cell);
1589 * fibre_data[cell_no].resize(n_q_points_cell);
1590 * for (unsigned int q_point_cell=0; q_point_cell<n_q_points_cell; ++q_point_cell)
1592 * const Point<dim> pt = fe_values.get_quadrature_points()[q_point_cell];
1593 * fibre_data[cell_no][q_point_cell] = MuscleFibre<dim>(bicep_geom.direction(pt,parameters.scale));
1598 * AssertThrow(false, ExcNotImplemented());
1604 * <a name="Linear_active_muscle_model.cc-LinearMuscleModelProblemupdate_fibre_state"></a>
1605 * <h4>LinearMuscleModelProblem::update_fibre_state</h4>
1611 * template <int dim>
1612 * double LinearMuscleModelProblem<dim>::get_neural_signal (const double time)
1616 * Note: 40 times less force generated than Martins2006
1617 * This is necessary due to the (compliant) linear tissue model
1620 * return (time > parameters.neural_signal_start_time && time < parameters.neural_signal_end_time ?
1625 * template <int dim>
1626 * void LinearMuscleModelProblem<dim>::update_fibre_activation (const double time)
1628 * const double u = get_neural_signal(time);
1630 * const unsigned int n_q_points_cell = qf_cell.size();
1631 * for (unsigned int cell=0; cell<triangulation.n_active_cells(); ++cell)
1633 * for (unsigned int q_point_cell=0; q_point_cell<n_q_points_cell; ++q_point_cell)
1635 * MuscleFibre<dim> &fibre = fibre_data[cell][q_point_cell];
1636 * fibre.update_alpha(u,dt);
1641 * template <int dim>
1642 * void LinearMuscleModelProblem<dim>::update_fibre_state ()
1644 * const unsigned int n_q_points_cell = qf_cell.size();
1646 * FEValues<dim> fe_values (fe, qf_cell, update_gradients);
1650 * Displacement gradient
1653 * std::vector< std::vector< Tensor<1,dim> > > u_grads (n_q_points_cell,
1654 * std::vector<Tensor<1,dim> >(dim));
1656 * unsigned int cell_no = 0;
1657 * for (typename DoFHandler<dim>::active_cell_iterator
1658 * cell = dof_handler.begin_active();
1659 * cell!=dof_handler.end(); ++cell, ++cell_no)
1661 * Assert(cell_no<fibre_data.size(), ExcMessage("Trying to access fibre data not stored for this cell index"));
1662 * fe_values.reinit(cell);
1663 * fe_values.get_function_gradients (solution, u_grads);
1665 * for (unsigned int q_point_cell=0; q_point_cell<n_q_points_cell; ++q_point_cell)
1667 * Assert(q_point_cell<fibre_data[cell_no].size(), ExcMessage("Trying to access fibre data not stored for this cell and qp index"));
1669 * const SymmetricTensor<2,dim> strain_tensor = get_small_strain (u_grads[q_point_cell]);
1670 * MuscleFibre<dim> &fibre = fibre_data[cell_no][q_point_cell];
1671 * fibre.update_state(strain_tensor, dt);
1679 * <a name="Linear_active_muscle_model.cc-LinearMuscleModelProblemsetup_system"></a>
1680 * <h4>LinearMuscleModelProblem::setup_system</h4>
1686 * template <int dim>
1687 * void LinearMuscleModelProblem<dim>::setup_system ()
1689 * dof_handler.distribute_dofs (fe);
1690 * hanging_node_constraints.clear ();
1691 * DoFTools::make_hanging_node_constraints (dof_handler,
1692 * hanging_node_constraints);
1693 * hanging_node_constraints.close ();
1694 * sparsity_pattern.reinit (dof_handler.n_dofs(),
1695 * dof_handler.n_dofs(),
1696 * dof_handler.max_couplings_between_dofs());
1697 * DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
1699 * hanging_node_constraints.condense (sparsity_pattern);
1701 * sparsity_pattern.compress();
1703 * system_matrix.reinit (sparsity_pattern);
1705 * solution.reinit (dof_handler.n_dofs());
1706 * system_rhs.reinit (dof_handler.n_dofs());
1708 * std::cout << " Number of active cells: "
1709 * << triangulation.n_active_cells()
1712 * std::cout << " Number of degrees of freedom: "
1713 * << dof_handler.n_dofs()
1720 * <a name="Linear_active_muscle_model.cc-LinearMuscleModelProblemassemble_system"></a>
1721 * <h4>LinearMuscleModelProblem::assemble_system</h4>
1727 * template <int dim>
1728 * SymmetricTensor<4,dim>
1729 * LinearMuscleModelProblem<dim>::get_stiffness_tensor (const unsigned int cell,
1730 * const unsigned int q_point_cell) const
1732 * static const SymmetricTensor<2,dim> I = unit_symmetric_tensor<dim>();
1734 * Assert(cell<fibre_data.size(), ExcMessage("Trying to access fibre data not stored for this cell index"));
1735 * Assert(q_point_cell<fibre_data[cell].size(), ExcMessage("Trying to access fibre data not stored for this cell and qp index"));
1736 * const MuscleFibre<dim> &fibre = fibre_data[cell][q_point_cell];
1743 * const double lambda = MuscleMatrix::lambda;
1744 * const double mu = MuscleMatrix::mu;
1750 * const double m_p = fibre.get_m_p();
1751 * const double m_s = fibre.get_m_s();
1752 * const double beta = fibre.get_beta(dt);
1753 * AssertThrow(beta != 0.0, ExcInternalError());
1754 * const double Cf = T0*(m_p + m_s*(1.0 - m_s/beta));
1755 * const Tensor<1,dim> &M = fibre.get_M();
1757 * SymmetricTensor<4,dim> C;
1758 * for (unsigned int i=0; i < dim; ++i)
1759 * for (unsigned int j=i; j < dim; ++j)
1760 * for (unsigned int k=0; k < dim; ++k)
1761 * for (unsigned int l=k; l < dim; ++l)
1765 * Matrix contribution
1768 * C[i][j][k][l] = lambda * I[i][j]*I[k][l]
1769 * + mu * (I[i][k]*I[j][l] + I[i][l]*I[j][k]);
1773 * Fibre contribution (Passive + active branches)
1776 * C[i][j][k][l] += Cf * M[i]*M[j]*M[k]*M[l];
1782 * template <int dim>
1783 * SymmetricTensor<2,dim>
1784 * LinearMuscleModelProblem<dim>::get_rhs_tensor (const unsigned int cell,
1785 * const unsigned int q_point_cell) const
1787 * Assert(cell<fibre_data.size(), ExcMessage("Trying to access fibre data not stored for this cell index"));
1788 * Assert(q_point_cell<fibre_data[cell].size(), ExcMessage("Trying to access fibre data not stored for this cell and qp index"));
1789 * const MuscleFibre<dim> &fibre = fibre_data[cell][q_point_cell];
1791 * const double m_s = fibre.get_m_s();
1792 * const double beta = fibre.get_beta(dt);
1793 * const double gamma = fibre.get_gamma(dt);
1794 * AssertThrow(beta != 0.0, ExcInternalError());
1795 * const double Sf = T0*(m_s*gamma/beta);
1796 * const Tensor<1,dim> &M = fibre.get_M();
1798 * SymmetricTensor<2,dim> S;
1799 * for (unsigned int i=0; i < dim; ++i)
1800 * for (unsigned int j=i; j < dim; ++j)
1804 * Fibre contribution (Active branch)
1807 * S[i][j] = Sf * M[i]*M[j];
1816 * <a name="Linear_active_muscle_model.cc-LinearMuscleModelProblemassemble_system"></a>
1817 * <h4>LinearMuscleModelProblem::assemble_system</h4>
1823 * template <int dim>
1824 * void LinearMuscleModelProblem<dim>::assemble_system (const double time)
1831 * system_matrix = 0;
1834 * FEValues<dim> fe_values (fe, qf_cell,
1835 * update_values | update_gradients |
1836 * update_quadrature_points | update_JxW_values);
1837 * FEFaceValues<dim> fe_face_values (fe, qf_face,
1839 * update_quadrature_points | update_JxW_values);
1841 * const unsigned int dofs_per_cell = fe.dofs_per_cell;
1842 * const unsigned int n_q_points_cell = qf_cell.size();
1843 * const unsigned int n_q_points_face = qf_face.size();
1845 * FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
1846 * Vector<double> cell_rhs (dofs_per_cell);
1848 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1855 * std::vector<Vector<double> > body_force_values (n_q_points_cell,
1856 * Vector<double>(dim));
1857 * std::vector<Vector<double> > traction_values (n_q_points_face,
1858 * Vector<double>(dim));
1860 * unsigned int cell_no = 0;
1861 * for (typename DoFHandler<dim>::active_cell_iterator
1862 * cell = dof_handler.begin_active();
1863 * cell!=dof_handler.end(); ++cell, ++cell_no)
1868 * fe_values.reinit (cell);
1869 * body_force.vector_value_list (fe_values.get_quadrature_points(),
1870 * body_force_values);
1872 * for (unsigned int q_point_cell=0; q_point_cell<n_q_points_cell; ++q_point_cell)
1874 * const SymmetricTensor<4,dim> C = get_stiffness_tensor (cell_no, q_point_cell);
1875 * const SymmetricTensor<2,dim> R = get_rhs_tensor(cell_no, q_point_cell);
1877 * for (unsigned int I=0; I<dofs_per_cell; ++I)
1879 * const unsigned int
1880 * component_I = fe.system_to_component_index(I).first;
1882 * for (unsigned int J=0; J<dofs_per_cell; ++J)
1884 * const unsigned int
1885 * component_J = fe.system_to_component_index(J).first;
1887 * for (unsigned int k=0; k < dim; ++k)
1888 * for (unsigned int l=0; l < dim; ++l)
1890 * += (fe_values.shape_grad(I,q_point_cell)[k] *
1891 * C[component_I][k][component_J][l] *
1892 * fe_values.shape_grad(J,q_point_cell)[l]) *
1893 * fe_values.JxW(q_point_cell);
1897 * for (unsigned int I=0; I<dofs_per_cell; ++I)
1899 * const unsigned int
1900 * component_I = fe.system_to_component_index(I).first;
1903 * += fe_values.shape_value(I,q_point_cell) *
1904 * body_force_values[q_point_cell](component_I) *
1905 * fe_values.JxW(q_point_cell);
1907 * for (unsigned int k=0; k < dim; ++k)
1909 * += (fe_values.shape_grad(I,q_point_cell)[k] *
1910 * R[component_I][k]) *
1911 * fe_values.JxW(q_point_cell);
1915 * for (unsigned int face = 0; face <GeometryInfo<dim>::faces_per_cell; ++face)
1917 * if (cell->face(face)->at_boundary() == true &&
1918 * ((parameters.problem == "IsotonicContraction" &&
1919 * cell->face(face)->boundary_id() == parameters.bid_CC_neumann) ||
1920 * (parameters.problem == "BicepsBrachii" &&
1921 * cell->face(face)->boundary_id() == parameters.bid_BB_neumann)) )
1923 * fe_face_values.reinit(cell, face);
1924 * traction.vector_value_list (fe_face_values.get_quadrature_points(),
1929 * Scale applied traction according to time
1932 * const double ramp = (time <= t_ramp_end ? time/t_ramp_end : 1.0);
1933 * Assert(ramp >= 0.0 && ramp <= 1.0, ExcMessage("Invalid force ramp"));
1934 * for (unsigned int q_point_face = 0; q_point_face < n_q_points_face; ++q_point_face)
1935 * traction_values[q_point_face] *= ramp;
1937 * for (unsigned int q_point_face = 0; q_point_face < n_q_points_face; ++q_point_face)
1939 * for (unsigned int I=0; I<dofs_per_cell; ++I)
1941 * const unsigned int
1942 * component_I = fe.system_to_component_index(I).first;
1945 * += fe_face_values.shape_value(I,q_point_face)*
1946 * traction_values[q_point_face][component_I]*
1947 * fe_face_values.JxW(q_point_face);
1953 * cell->get_dof_indices (local_dof_indices);
1954 * for (unsigned int i=0; i<dofs_per_cell; ++i)
1956 * for (unsigned int j=0; j<dofs_per_cell; ++j)
1957 * system_matrix.add (local_dof_indices[i],
1958 * local_dof_indices[j],
1959 * cell_matrix(i,j));
1961 * system_rhs(local_dof_indices[i]) += cell_rhs(i);
1965 * hanging_node_constraints.condense (system_matrix);
1966 * hanging_node_constraints.condense (system_rhs);
1969 * template <int dim>
1970 * void LinearMuscleModelProblem<dim>::apply_boundary_conditions ()
1972 * std::map<types::global_dof_index,double> boundary_values;
1974 * if (parameters.problem == "IsotonicContraction")
1978 * Symmetry condition on -X faces
1982 * ComponentMask component_mask_x (dim, false);
1983 * component_mask_x.set(0, true);
1984 * VectorTools::interpolate_boundary_values (dof_handler,
1985 * parameters.bid_CC_dirichlet_symm_X,
1986 * Functions::ZeroFunction<dim>(dim),
1988 * component_mask_x);
1992 * Symmetry condition on -Z/+Z faces
1996 * ComponentMask component_mask_z (dim, false);
1997 * component_mask_z.set(2, true);
1998 * VectorTools::interpolate_boundary_values (dof_handler,
1999 * parameters.bid_CC_dirichlet_symm_Z,
2000 * Functions::ZeroFunction<dim>(dim),
2002 * component_mask_z);
2006 * Fixed point on -X face
2010 * const Point<dim> fixed_point (-parameters.half_length_x,0.0,0.0);
2011 * std::vector<types::global_dof_index> fixed_dof_indices;
2012 * bool found_point_of_interest = false;
2014 * for (typename DoFHandler<dim>::active_cell_iterator
2015 * cell = dof_handler.begin_active(),
2016 * endc = dof_handler.end(); cell != endc; ++cell)
2018 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
2022 * We know that the fixed point is on the -X Dirichlet boundary
2025 * if (cell->face(face)->at_boundary() == true &&
2026 * cell->face(face)->boundary_id() == parameters.bid_CC_dirichlet_symm_X)
2028 * for (unsigned int face_vertex_index = 0; face_vertex_index < GeometryInfo<dim>::vertices_per_face; ++face_vertex_index)
2030 * if (cell->face(face)->vertex(face_vertex_index).distance(fixed_point) < 1e-6)
2032 * found_point_of_interest = true;
2033 * for (unsigned int index_component = 0; index_component < dim; ++index_component)
2034 * fixed_dof_indices.push_back(cell->face(face)->vertex_dof_index(face_vertex_index,
2035 * index_component));
2038 * if (found_point_of_interest == true) break;
2041 * if (found_point_of_interest == true) break;
2043 * if (found_point_of_interest == true) break;
2046 * Assert(found_point_of_interest == true, ExcMessage("Didn't find
point of interest"));
2047 * AssertThrow(fixed_dof_indices.size() == dim, ExcMessage("Didn't find the correct number of DoFs to fix"));
2049 * for (unsigned int i=0; i < fixed_dof_indices.size(); ++i)
2050 * boundary_values[fixed_dof_indices[i]] = 0.0;
2053 * else if (parameters.problem == "BicepsBrachii")
2055 * if (parameters.include_gravity == false)
2059 * Symmetry condition on -X surface
2063 * ComponentMask component_mask_x (dim, false);
2064 * component_mask_x.set(0, true);
2065 * VectorTools::interpolate_boundary_values (dof_handler,
2066 * parameters.bid_BB_dirichlet_X,
2067 * Functions::ZeroFunction<dim>(dim),
2069 * component_mask_x);
2074 * Fixed central point on -X surface
2078 * const Point<dim> fixed_point (0.0,0.0,0.0);
2079 * std::vector<types::global_dof_index> fixed_dof_indices;
2080 * bool found_point_of_interest = false;
2082 * for (typename DoFHandler<dim>::active_cell_iterator
2083 * cell = dof_handler.begin_active(),
2084 * endc = dof_handler.end(); cell != endc; ++cell)
2086 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
2090 * We know that the fixed point is on the -X Dirichlet boundary
2093 * if (cell->face(face)->at_boundary() == true &&
2094 * cell->face(face)->boundary_id() == parameters.bid_BB_dirichlet_X)
2096 * for (unsigned int face_vertex_index = 0; face_vertex_index < GeometryInfo<dim>::vertices_per_face; ++face_vertex_index)
2098 * if (cell->face(face)->vertex(face_vertex_index).distance(fixed_point) < 1e-6)
2100 * found_point_of_interest = true;
2101 * for (unsigned int index_component = 0; index_component < dim; ++index_component)
2102 * fixed_dof_indices.push_back(cell->face(face)->vertex_dof_index(face_vertex_index,
2103 * index_component));
2106 * if (found_point_of_interest == true) break;
2109 * if (found_point_of_interest == true) break;
2111 * if (found_point_of_interest == true) break;
2114 * Assert(found_point_of_interest == true, ExcMessage("Didn't find
point of interest"));
2115 * AssertThrow(fixed_dof_indices.size() == dim, ExcMessage("Didn't find the correct number of DoFs to fix"));
2117 * for (unsigned int i=0; i < fixed_dof_indices.size(); ++i)
2118 * boundary_values[fixed_dof_indices[i]] = 0.0;
2125 * When we apply gravity, some additional constraints
2126 * are required to support the load of the muscle, as
2127 * the material response is more compliant than would
2128 * be the case in reality.
2132 * Symmetry condition on -X surface
2136 * ComponentMask component_mask_x (dim, true);
2137 * VectorTools::interpolate_boundary_values (dof_handler,
2138 * parameters.bid_BB_dirichlet_X,
2139 * Functions::ZeroFunction<dim>(dim),
2141 * component_mask_x);
2145 * Symmetry condition on -X surface
2149 * ComponentMask component_mask_x (dim, false);
2150 * component_mask_x.set(1, true);
2151 * component_mask_x.set(2, true);
2152 * VectorTools::interpolate_boundary_values (dof_handler,
2153 * parameters.bid_BB_neumann,
2154 * Functions::ZeroFunction<dim>(dim),
2156 * component_mask_x);
2162 * Roller condition at central point on +X face
2166 * const Point<dim> roller_point (parameters.axial_length*parameters.scale,0.0,0.0);
2167 * std::vector<types::global_dof_index> fixed_dof_indices;
2168 * bool found_point_of_interest = false;
2170 * for (typename DoFHandler<dim>::active_cell_iterator
2171 * cell = dof_handler.begin_active(),
2172 * endc = dof_handler.end(); cell != endc; ++cell)
2174 * for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
2178 * We know that the fixed point is on the +X Neumann boundary
2181 * if (cell->face(face)->at_boundary() == true &&
2182 * cell->face(face)->boundary_id() == parameters.bid_BB_neumann)
2184 * for (unsigned int face_vertex_index = 0; face_vertex_index < GeometryInfo<dim>::vertices_per_face; ++face_vertex_index)
2186 * if (cell->face(face)->vertex(face_vertex_index).distance(roller_point) < 1e-6)
2188 * found_point_of_interest = true;
2189 * for (unsigned int index_component = 1; index_component < dim; ++index_component)
2190 * fixed_dof_indices.push_back(cell->face(face)->vertex_dof_index(face_vertex_index,
2191 * index_component));
2194 * if (found_point_of_interest == true) break;
2197 * if (found_point_of_interest == true) break;
2199 * if (found_point_of_interest == true) break;
2202 * Assert(found_point_of_interest == true, ExcMessage("Didn't find
point of interest"));
2203 * AssertThrow(fixed_dof_indices.size() == dim-1, ExcMessage("Didn't find the correct number of DoFs to fix"));
2205 * for (unsigned int i=0; i < fixed_dof_indices.size(); ++i)
2206 * boundary_values[fixed_dof_indices[i]] = 0.0;
2210 * AssertThrow(false, ExcNotImplemented());
2212 * MatrixTools::apply_boundary_values (boundary_values,
2222 * <a name="Linear_active_muscle_model.cc-LinearMuscleModelProblemsolve"></a>
2223 * <h4>LinearMuscleModelProblem::solve</h4>
2229 * template <int dim>
2230 * void LinearMuscleModelProblem<dim>::solve ()
2232 * SolverControl solver_control (system_matrix.m(), 1e-12);
2233 * SolverCG<> cg (solver_control);
2235 * PreconditionSSOR<> preconditioner;
2236 * preconditioner.initialize(system_matrix, 1.2);
2238 * cg.solve (system_matrix, solution, system_rhs,
2241 * hanging_node_constraints.distribute (solution);
2248 * <a name="Linear_active_muscle_model.cc-LinearMuscleModelProblemoutput_results"></a>
2249 * <h4>LinearMuscleModelProblem::output_results</h4>
2258 * template <int dim>
2259 * void LinearMuscleModelProblem<dim>::output_results (const unsigned int timestep,
2260 * const double time) const
2264 * Visual output: FEM results
2268 * std::string filename = "solution-";
2269 * filename += Utilities::int_to_string(timestep,4);
2270 * filename += ".vtk";
2271 * std::ofstream output (filename.c_str());
2273 * DataOut<dim> data_out;
2274 * data_out.attach_dof_handler (dof_handler);
2276 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
2277 * data_component_interpretation(dim,
2278 * DataComponentInterpretation::component_is_part_of_vector);
2279 * std::vector<std::string> solution_name(dim, "displacement");
2281 * data_out.add_data_vector (solution, solution_name,
2282 * DataOut<dim>::type_dof_data,
2283 * data_component_interpretation);
2284 * data_out.build_patches ();
2285 * data_out.write_vtk (output);
2290 * Visual output: FEM data
2294 * std::string filename = "fibres-";
2295 * filename += Utilities::int_to_string(timestep,4);
2296 * filename += ".vtk";
2297 * std::ofstream output (filename.c_str());
2300 * << "# vtk DataFile Version 3.0" << std::endl
2301 * << "# " << std::endl
2302 * << "ASCII"<< std::endl
2303 * << "DATASET POLYDATA"<< std::endl << std::endl;
2307 * Extract fibre data from quadrature points
2310 * const unsigned int n_cells = triangulation.n_active_cells();
2311 * const unsigned int n_q_points_cell = qf_cell.size();
2351 *
cell = dof_handler.begin_active();
2352 *
cell != dof_handler.end();
2355 *
fe_values.reinit (cell);
2356 *
fe_values.get_function_values (solution,
u_values);
2357 *
fe_values.get_function_gradients (solution,
u_grads);
2363 *
for (
unsigned int d=0;
d<dim; ++
d)
2393 *
<<
" float" << std::endl;
2396 *
for (
unsigned int j=0;
j < dim; ++
j)
2400 *
output << std::endl;
2408 *
output <<
"\nPOINT_DATA "
2410 *
<< std::endl << std::endl;
2418 *
<<
"VECTORS displacement float"
2422 *
for (
unsigned int j=0;
j < dim; ++
j)
2426 *
output << std::endl;
2428 *
output << std::endl;
2436 *
<<
"VECTORS direction float"
2440 *
for (
unsigned int j=0;
j < dim; ++
j)
2444 *
output << std::endl;
2446 *
output << std::endl;
2453 *
for (
unsigned int v=0; v <
n_results; ++v)
2458 *
<<
" float 1" << std::endl
2459 *
<<
"LOOKUP_TABLE default "
2465 *
output << std::endl;
2476 *
Point<dim>(parameters.half_length_x, 0.0, 0.0) :
2485 *
cell = dof_handler.begin_active(),
2486 *
endc = dof_handler.end(); cell !=
endc; ++cell)
2495 *
if (cell->face(face)->at_boundary() ==
true &&
2496 *
((parameters.problem ==
"IsotonicContraction" &&
2497 *
cell->face(face)->boundary_id() == parameters.bid_CC_neumann) ||
2498 *
(parameters.problem ==
"BicepsBrachii" &&
2499 *
cell->face(face)->boundary_id() == parameters.bid_BB_neumann)) )
2522 *
const std::string
filename =
"displacement_POI.csv";
2523 *
std::ofstream output;
2526 *
output.open(
filename.c_str(), std::ofstream::out);
2528 *
<<
"Time [s]" <<
"," <<
"X-displacement [mm]" << std::endl;
2531 *
output.open(
filename.c_str(), std::ios_base::app);
2546 * <a name=
"Linear_active_muscle_model.cc-LinearMuscleModelProblemrun"></a>
2547 * <
h4>LinearMuscleModelProblem::run</
h4>
2553 *
template <
int dim>
2565 *
double time = 0.0;
2570 *
<<
" @ time " << time
2626 * <a name=
"Linear_active_muscle_model.cc-Thecodemaincodefunction"></a>
2637 * ::deallog.depth_console (0);
2638 *
const unsigned int dim = 3;
2643 *
catch (std::exception &exc)
2645 *
std::cerr << std::endl << std::endl
2646 *
<<
"----------------------------------------------------"
2648 *
std::cerr <<
"Exception on processing: " << std::endl
2649 *
<< exc.what() << std::endl
2650 *
<<
"Aborting!" << std::endl
2651 *
<<
"----------------------------------------------------"
2658 *
std::cerr << std::endl << std::endl
2659 *
<<
"----------------------------------------------------"
2661 *
std::cerr <<
"Unknown exception!" << std::endl
2662 *
<<
"Aborting!" << std::endl
2663 *
<<
"----------------------------------------------------"
numbers::NumberTraits< Number >::real_type norm() const
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
@ update_values
Shape function values.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result, const bool copy_manifold_ids=false, const std::vector< types::manifold_id > &manifold_priorities={})
void hyper_ball(Triangulation< dim, spacedim > &tria, const Point< spacedim > ¢er={}, const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
void flatten_triangulation(const Triangulation< dim, spacedim1 > &in_tria, Triangulation< dim, spacedim2 > &out_tria)
@ matrix
Contents is actually a matrix.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * end(VectorType &V)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
long double gamma(const unsigned int n)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
constexpr types::global_dof_index invalid_dof_index
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Function &function, const unsigned int grainsize)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation