Skip to content

Commit

Permalink
Merge pull request #5641 from gassmoeller/use_unroll_for_tensor
Browse files Browse the repository at this point in the history
Use constructor and unroll function for tensors
  • Loading branch information
bangerth authored May 30, 2024
2 parents 48bbd7a + b4e9f78 commit d22bb44
Show file tree
Hide file tree
Showing 6 changed files with 24 additions and 31 deletions.
11 changes: 5 additions & 6 deletions cookbooks/finite_strain/doc/finite_strain.cc
Original file line number Diff line number Diff line change
Expand Up @@ -18,17 +18,16 @@
<http://www.gnu.org/licenses/>.
*/

for (unsigned int q=0; q < in.position.size(); ++q)
for (unsigned int q=0; q < in.n_evaluation_points(); ++q)
{
// Convert the compositional fields into the tensor quantity they represent.
Tensor<2,dim> strain;
for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
strain[Tensor<2,dim>::unrolled_to_component_indices(i)] = in.composition[q][i];
const Tensor<2,dim> strain(make_array_view(&in.composition[q][0],
&in.composition[q][0] + Tensor<2,dim>::n_independent_components));

// Compute the strain accumulated in this timestep.
const Tensor<2,dim> strain_increment = this->get_timestep() * (velocity_gradients[q] * strain);

// Output the strain increment component-wise to its respective compositional field's reaction terms.
for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
out.reaction_terms[q][i] = strain_increment[Tensor<2,dim>::unrolled_to_component_indices(i)];
strain_increment.unroll(&out.reaction_terms[q][0],
&out.reaction_terms[q][0] + Tensor<2,dim>::n_independent_components);
}
9 changes: 4 additions & 5 deletions cookbooks/finite_strain/finite_strain.cc
Original file line number Diff line number Diff line change
Expand Up @@ -80,16 +80,15 @@ namespace aspect
for (unsigned int q=0; q < in.n_evaluation_points(); ++q)
{
// Convert the compositional fields into the tensor quantity they represent.
Tensor<2,dim> strain;
for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
strain[Tensor<2,dim>::unrolled_to_component_indices(i)] = in.composition[q][i];
const Tensor<2,dim> strain(make_array_view(&in.composition[q][0],
&in.composition[q][0] + Tensor<2,dim>::n_independent_components));

// Compute the strain accumulated in this timestep.
const Tensor<2,dim> strain_increment = this->get_timestep() * (velocity_gradients[q] * strain);

// Output the strain increment component-wise to its respective compositional field's reaction terms.
for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
out.reaction_terms[q][i] = strain_increment[Tensor<2,dim>::unrolled_to_component_indices(i)];
strain_increment.unroll(&out.reaction_terms[q][0],
&out.reaction_terms[q][0] + Tensor<2,dim>::n_independent_components);
}
}
}
Expand Down
5 changes: 2 additions & 3 deletions source/material_model/rheology/strain_dependent.cc
Original file line number Diff line number Diff line change
Expand Up @@ -395,9 +395,8 @@ namespace aspect
case finite_strain_tensor:
{
// Calculate second invariant of left stretching tensor "L"
Tensor<2,dim> strain;
for (unsigned int q = 0; q < Tensor<2,dim>::n_independent_components ; ++q)
strain[Tensor<2,dim>::unrolled_to_component_indices(q)] = composition[q];
const Tensor<2,dim> strain(make_array_view(&composition[0],
&composition[0] + Tensor<2,dim>::n_independent_components));
const SymmetricTensor<2,dim> L = symmetrize( strain * transpose(strain) );

const double strain_ii = std::fabs(second_invariant(L));
Expand Down
12 changes: 5 additions & 7 deletions source/particle/property/integrated_strain.cc
Original file line number Diff line number Diff line change
Expand Up @@ -43,11 +43,8 @@ namespace aspect
const std::vector<Tensor<1,dim>> &gradients,
typename ParticleHandler<dim>::particle_iterator &particle) const
{
const auto data = particle->get_properties();

Tensor<2,dim> old_strain;
for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
old_strain[Tensor<2,dim>::unrolled_to_component_indices(i)] = data[data_position + i];
const Tensor<2,dim> old_strain(make_array_view(&particle->get_properties()[data_position],
&particle->get_properties()[data_position] + Tensor<2,dim>::n_independent_components));

Tensor<2,dim> grad_u;
for (unsigned int d=0; d<dim; ++d)
Expand Down Expand Up @@ -75,8 +72,9 @@ namespace aspect
// strain of the current time step
new_strain = old_strain + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;

for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
data[data_position + i] = new_strain[Tensor<2,dim>::unrolled_to_component_indices(i)];
// unroll and store the new strain
new_strain.unroll(&particle->get_properties()[data_position],
&particle->get_properties()[data_position] + Tensor<2,dim>::n_independent_components);
}

template <int dim>
Expand Down
9 changes: 4 additions & 5 deletions tests/simple_shear.cc
Original file line number Diff line number Diff line change
Expand Up @@ -80,16 +80,15 @@ namespace aspect
for (unsigned int q=0; q < in.n_evaluation_points(); ++q)
{
// Convert the compositional fields into the tensor quantity they represent.
Tensor<2,dim> strain;
for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
strain[Tensor<2,dim>::unrolled_to_component_indices(i)] = in.composition[q][i];
const Tensor<2,dim> strain(make_array_view(&in.composition[q][0],
&in.composition[q][0] + Tensor<2,dim>::n_independent_components));

// Compute the strain accumulated in this timestep.
const Tensor<2,dim> strain_increment = this->get_timestep() * (velocity_gradients[q] * strain);

// Output the strain increment component-wise to its respective compositional field's reaction terms.
for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
out.reaction_terms[q][i] = strain_increment[Tensor<2,dim>::unrolled_to_component_indices(i)];
strain_increment.unroll(&out.reaction_terms[q][0],
&out.reaction_terms[q][0] + Tensor<2,dim>::n_independent_components);
}
}
}
Expand Down
9 changes: 4 additions & 5 deletions tests/simple_shear_output_the_mobility.cc
Original file line number Diff line number Diff line change
Expand Up @@ -80,16 +80,15 @@ namespace aspect
for (unsigned int q=0; q < in.n_evaluation_points(); ++q)
{
// Convert the compositional fields into the tensor quantity they represent.
Tensor<2,dim> strain;
for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
strain[Tensor<2,dim>::unrolled_to_component_indices(i)] = in.composition[q][i];
const Tensor<2,dim> strain(make_array_view(&in.composition[q][0],
&in.composition[q][0] + Tensor<2,dim>::n_independent_components));

// Compute the strain accumulated in this timestep.
const Tensor<2,dim> strain_increment = this->get_timestep() * (velocity_gradients[q] * strain);

// Output the strain increment component-wise to its respective compositional field's reaction terms.
for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
out.reaction_terms[q][i] = strain_increment[Tensor<2,dim>::unrolled_to_component_indices(i)];
strain_increment.unroll(&out.reaction_terms[q][0],
&out.reaction_terms[q][0] + Tensor<2,dim>::n_independent_components);
}
}
}
Expand Down

0 comments on commit d22bb44

Please sign in to comment.