Skip to content

Commit

Permalink
Restricted mean muscle velocity and strain rate in output_gearing_inf…
Browse files Browse the repository at this point in the history
…o to slab in the middle
  • Loading branch information
javieralmonacid committed Apr 24, 2024
1 parent ecd4b22 commit 04a26f9
Showing 1 changed file with 37 additions and 23 deletions.
60 changes: 37 additions & 23 deletions dynamic-muscle.cc
Original file line number Diff line number Diff line change
Expand Up @@ -4974,40 +4974,52 @@ namespace Flexodeal
if (parameters.type_of_simulation != "dynamic")
return void();

double mean_muscle_velocity = 0.0, mean_strain_rate = 0.0;
double mean_muscle_velocity = 0.0, mean_strain_rate = 0.0, volume_slab = 0.0;

FEValues<dim> fe_values(fe, qf_cell,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);

for (const auto &cell : triangulation.active_cell_iterators())
{
fe_values.reinit(cell);

const std::vector<std::shared_ptr<const PointHistory<dim> > > lqph =
quadrature_point_history.get_data(cell);
Assert(lqph.size() == n_q_points, ExcInternalError());

for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
// We restrict the computation of these quantities to a slab in the
// middle of the domain to avoid averaging with outliers located
// at the ends of the block.
if (cell->center()[0] >= 3.0*parameters.length/8.0 &&
cell->center()[0] <= 5.0*parameters.length/8.0)
{
// As in output_energies(), get_velocity_previous returns the current velocity.
const Tensor<1, dim> muscle_velocity = lqph[q_point]->get_velocity_previous();
// We now retrieve information related to the fibre velocity.
const double strain_rate = lqph[q_point]->get_strain_rate();
// The rest of the quantities are related to the integrals themselves, as usual.
const double det_F = lqph[q_point]->get_det_F();
const double JxW = fe_values.JxW(q_point);
fe_values.reinit(cell);

const std::vector<std::shared_ptr<const PointHistory<dim> > > lqph =
quadrature_point_history.get_data(cell);
Assert(lqph.size() == n_q_points, ExcInternalError());

mean_muscle_velocity += muscle_velocity.norm() * det_F * JxW;
mean_strain_rate += strain_rate * det_F * JxW;
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
// As in output_energies(), get_velocity_previous returns the current velocity.
const Tensor<1, dim> muscle_velocity = lqph[q_point]->get_velocity_previous();
// We now retrieve information related to the fibre velocity.
const double strain_rate = lqph[q_point]->get_strain_rate();
// The rest of the quantities are related to the integrals themselves, as usual.
const double det_F = lqph[q_point]->get_det_F();
const double JxW = fe_values.JxW(q_point);

// We also decide to output the muscle velocity only in the x component,
// since this is the default line of action. If the line of action is changed,
// then this quantity must be updated properly.
mean_muscle_velocity += muscle_velocity[0] * det_F * JxW;
mean_strain_rate += strain_rate * det_F * JxW;
volume_slab += det_F * JxW;
}
}
}

const double current_volume = compute_vol_current();
mean_muscle_velocity = mean_muscle_velocity / current_volume;
mean_strain_rate = mean_strain_rate / current_volume;
mean_muscle_velocity = mean_muscle_velocity / volume_slab;
mean_strain_rate = mean_strain_rate / volume_slab;

const double initial_fibre_length = parameters.length / parameters.muscle_fibre_orientation_x;
// We define the initial fibre length as
// $L_f^0 = \dfrac{H}{\sin (\beta_0)} = \dfrac{H}{a_{0,z}}$
const double initial_fibre_length = parameters.height / parameters.muscle_fibre_orientation_z;
const double strain_rate_naught = parameters.max_strain_rate;

// Output time series:
Expand Down Expand Up @@ -5035,7 +5047,8 @@ namespace Flexodeal
<< "," << "Mean muscle velocity [m/s]"
<< "," << "Mean fibre strain rate (non-dim)"
<< "," << "Initial fibre length [m]"
<< "," << "Maximum strain rate [1/s]" << "\n";
<< "," << "Maximum strain rate [1/s]"
<< "," << "Volume slab [m^3]" << "\n";
}
else
output.open(filename.str(), std::ios_base::app);
Expand All @@ -5045,7 +5058,8 @@ namespace Flexodeal
<< "," << mean_muscle_velocity
<< "," << mean_strain_rate
<< "," << initial_fibre_length
<< "," << strain_rate_naught << "\n";
<< "," << strain_rate_naught
<< "," << volume_slab << "\n";
}

template <int dim>
Expand Down

0 comments on commit 04a26f9

Please sign in to comment.