Skip to content

Commit 08e245e

Browse files
MFratersMFraters
committed
Add direction dilation.
1 parent 25755ff commit 08e245e

File tree

5 files changed

+66
-32
lines changed

5 files changed

+66
-32
lines changed

‎include/aspect/material_model/interface.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1232,7 +1232,7 @@ namespace aspect
12321232
* A scalar value per evaluation point that specifies the prescribed dilation
12331233
* in that point.
12341234
*/
1235-
std::vector<double> dilation;
1235+
std::vector<std::vector<double>> dilation;
12361236
};
12371237

12381238

‎source/material_model/interface.cc‎

Lines changed: 13 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1000,8 +1000,8 @@ namespace aspect
10001000

10011001
template <int dim>
10021002
PrescribedPlasticDilation<dim>::PrescribedPlasticDilation (const unsigned int n_points)
1003-
: NamedAdditionalMaterialOutputs<dim>(std::vector<std::string>(1, "prescribed_dilation")),
1004-
dilation(n_points, numbers::signaling_nan<double>())
1003+
: NamedAdditionalMaterialOutputs<dim>(std::vector<std::string>(3, "prescribed_dilation")),
1004+
dilation(dim, std::vector<double>(n_points, numbers::signaling_nan<double>()))
10051005
{}
10061006

10071007

@@ -1010,8 +1010,17 @@ namespace aspect
10101010
std::vector<double> PrescribedPlasticDilation<dim>::get_nth_output(const unsigned int idx) const
10111011
{
10121012
(void)idx;
1013-
Assert(idx==0, ExcInternalError());
1014-
return dilation;
1013+
Assert(idx<3, ExcInternalError());
1014+
switch (idx)
1015+
{
1016+
case 0:
1017+
return dilation[0];
1018+
case 1:
1019+
return dilation[1];
1020+
case 2:
1021+
return dilation[2];
1022+
}
1023+
return std::vector<double>();
10151024
}
10161025

10171026

‎source/simulator/assemblers/newton_stokes.cc‎

Lines changed: 23 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -437,22 +437,32 @@ namespace aspect
437437
* JxW;
438438

439439
if (enable_prescribed_dilation)
440-
data.local_rhs(i) += (
441-
// RHS of - (div u,q) = - (R,q)
442-
- pressure_scaling
443-
* prescribed_dilation->dilation[q]
444-
* scratch.phi_p[i]
445-
) * JxW;
440+
{
441+
const unsigned int index_direction=fe.system_to_component_index(i).first;
442+
// Only want the velocity components and not the pressure one (which is the last one), so add 1
443+
if (introspection.is_stokes_component(index_direction+1))
444+
data.local_rhs(i) += (
445+
// RHS of - (div u,q) = - (R,q)
446+
- pressure_scaling
447+
* prescribed_dilation->dilation[index_direction][q]
448+
* scratch.phi_p[i]
449+
) * JxW;
450+
}
446451

447452
// Only assemble this term if we are running incompressible, otherwise this term
448453
// is already included on the LHS of the equation.
449454
if (enable_prescribed_dilation && !material_model_is_compressible)
450-
data.local_rhs(i) += (
451-
// RHS of momentum eqn: - \int 2/3 eta R, div v
452-
- 2.0 / 3.0 * eta
453-
* prescribed_dilation->dilation[q]
454-
* scratch.div_phi_u[i]
455-
) * JxW;
455+
{
456+
const unsigned int index_direction=fe.system_to_component_index(i).first;
457+
// Only want the velocity components and not the pressure one (which is the last one), so add 1
458+
if (introspection.is_stokes_component(index_direction+1))
459+
data.local_rhs(i) += (
460+
// RHS of momentum eqn: - \int 2/3 eta R, div v
461+
- 2.0 / 3.0 * eta
462+
* prescribed_dilation->dilation[index_direction][q]
463+
* scratch.div_phi_u[i]
464+
) * JxW;
465+
}
456466
}
457467

458468
// and then the matrix, if necessary
@@ -623,7 +633,7 @@ namespace aspect
623633

624634
Assert(!this->get_parameters().enable_prescribed_dilation
625635
||
626-
outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation.size()
636+
outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation[0].size()
627637
== n_points, ExcInternalError());
628638

629639
if (this->get_newton_handler().parameters.newton_derivative_scaling_factor != 0)

‎source/simulator/assemblers/stokes.cc‎

Lines changed: 23 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -418,22 +418,32 @@ namespace aspect
418418
* JxW;
419419

420420
if (prescribed_dilation != nullptr)
421-
data.local_rhs(i) += (
422-
// RHS of - (div u,q) = - (R,q)
423-
- pressure_scaling
424-
* prescribed_dilation->dilation[q]
425-
* scratch.phi_p[i]
426-
) * JxW;
421+
{
422+
const unsigned int index_direction=fe.system_to_component_index(i).first;
423+
// Only want the velocity components and not the pressure one (which is the last one), so add 1
424+
if (introspection.is_stokes_component(index_direction+1))
425+
data.local_rhs(i) += (
426+
// RHS of - (div u,q) = - (R,q)
427+
- pressure_scaling
428+
* prescribed_dilation->dilation[index_direction][q]
429+
* scratch.phi_p[i]
430+
) * JxW;
431+
}
427432

428433
// Only assemble this term if we are running incompressible, otherwise this term
429434
// is already included on the LHS of the equation.
430435
if (prescribed_dilation != nullptr && !material_model_is_compressible)
431-
data.local_rhs(i) += (
432-
// RHS of momentum eqn: - \int 2/3 eta R, div v
433-
- 2.0 / 3.0 * eta
434-
* prescribed_dilation->dilation[q]
435-
* scratch.div_phi_u[i]
436-
) * JxW;
436+
{
437+
const unsigned int index_direction=fe.system_to_component_index(i).first;
438+
// Only want the velocity components and not the pressure one (which is the last one), so add 1
439+
if (introspection.is_stokes_component(index_direction+1))
440+
data.local_rhs(i) += (
441+
// RHS of momentum eqn: - \int 2/3 eta R, div v
442+
- 2.0 / 3.0 * eta
443+
* prescribed_dilation->dilation[index_direction][q]
444+
* scratch.div_phi_u[i]
445+
) * JxW;
446+
}
437447

438448
if (scratch.rebuild_stokes_matrix)
439449
for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
@@ -501,7 +511,7 @@ namespace aspect
501511

502512
Assert(!this->get_parameters().enable_prescribed_dilation
503513
||
504-
outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation.size()
514+
outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation[0].size()
505515
== n_points, ExcInternalError());
506516

507517
// Elasticity:

‎tests/prescribed_dilation.cc‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -208,7 +208,12 @@ namespace aspect
208208
}
209209
if (prescribed_dilation)
210210
{
211-
prescribed_dilation->dilation[i] = x;
211+
prescribed_dilation->dilation[0][i] = x;
212+
prescribed_dilation->dilation[1][i] = x;
213+
if (dim == 3)
214+
{
215+
prescribed_dilation->dilation[2][i] = x;
216+
}
212217
}
213218

214219
}

0 commit comments

Comments
 (0)