Skip to content

Commit

Permalink
Address feature review comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@xuchenhan-tri
xuchenhan-tri committed Dec 13, 2024
1 parent aff1947 commit 33f2d23
Show file tree
Hide file tree
Showing 3 changed files with 130 additions and 86 deletions.
73 changes: 42 additions & 31 deletions multibody/mpm/particle_data.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ ConstitutiveModelVariant<T> MakeConstitutiveModel(
}

template <typename T>
DeformationGradientDataVariant<T> MakeStrainData(
DeformationGradientDataVariant<T> MakeDeformationGradientData(
const fem::DeformableBodyConfig<double>& config) {
switch (config.material_model()) {
case fem::MaterialModel::kCorotated:
Expand All @@ -53,12 +53,14 @@ T ParticleData<T>::ComputeTotalEnergy(const std::vector<Matrix3<T>>& F,
[&, this](auto&& model) {
const Matrix3<T>& F_p = F[p];
const Matrix3<T>& F0_p = F0[p];
using StrainDataType = typename std::decay_t<decltype(model)>::Data;
StrainDataType& strain_data_p =
std::get<StrainDataType>(strain_data_[p]);
strain_data_p.UpdateData(F_p, F0_p);
using DeformationGradientDataType =
typename std::remove_reference_t<decltype(model)>::Data;
DeformationGradientDataType& deformation_gradient_data_p =
std::get<DeformationGradientDataType>(
deformation_gradient_data_[p]);
deformation_gradient_data_p.UpdateData(F_p, F0_p);
T Psi;
model.CalcElasticEnergyDensity(strain_data_p, &Psi);
model.CalcElasticEnergyDensity(deformation_gradient_data_p, &Psi);
result += Psi * volume_[p];
},
constitutive_models_[p]);
Expand All @@ -80,13 +82,16 @@ void ParticleData<T>::ComputeKirchhoffStress(
[&, this](auto& model) {
const Matrix3<T>& F_p = F[p];
const Matrix3<T>& F0_p = F0[p];
using StrainDataType = typename std::decay_t<decltype(model)>::Data;
StrainDataType& strain_data_p =
std::get<StrainDataType>(strain_data_[p]);
strain_data_p.UpdateData(F_p, F0_p);
using DeformationGradientDataType =
typename std::decay_t<decltype(model)>::Data;
DeformationGradientDataType& deformation_gradient_data_p =
std::get<DeformationGradientDataType>(
deformation_gradient_data_[p]);
deformation_gradient_data_p.UpdateData(F_p, F0_p);
auto& tau_volume_p = (*volume_scaled_stress)[p];
const Matrix3<T>& particle_F = F_[p];
model.CalcFirstPiolaStress(strain_data_p, &tau_volume_p);
model.CalcFirstPiolaStress(deformation_gradient_data_p,
&tau_volume_p);
tau_volume_p *= volume_[p] * particle_F.transpose();
},
constitutive_models_[p]);
Expand All @@ -107,44 +112,50 @@ void ParticleData<T>::ComputePK1StressDerivatives(
[&, this](auto& model) {
const Matrix3<T>& F_p = F[p];
const Matrix3<T>& F0_p = F0[p];
using StrainDataType = typename std::decay_t<decltype(model)>::Data;
StrainDataType& strain_data_p =
std::get<StrainDataType>(strain_data_[p]);
strain_data_p.UpdateData(F_p, F0_p);
using DeformationGradientDataType =
typename std::decay_t<decltype(model)>::Data;
DeformationGradientDataType& deformation_gradient_data_p =
std::get<DeformationGradientDataType>(
deformation_gradient_data_[p]);
deformation_gradient_data_p.UpdateData(F_p, F0_p);
auto& dPdF_volume_p = (*dPdF_volume)[p];
model.CalcFirstPiolaStressDerivative(strain_data_p, &dPdF_volume_p);
model.CalcFirstPiolaStressDerivative(deformation_gradient_data_p,
&dPdF_volume_p);
dPdF_volume_p.mutable_data() *= volume_[p];
},
constitutive_models_[p]);
}
}

template <typename T>
void ParticleData<T>::Sample(const std::vector<Vector3<double>>& positions,
double total_volume,
const fem::DeformableBodyConfig<double>& config) {
void ParticleData<T>::AddParticles(
const std::vector<Vector3<double>>& positions, double total_volume,
const fem::DeformableBodyConfig<double>& config) {
DRAKE_THROW_UNLESS(total_volume > 0);
DRAKE_THROW_UNLESS(!positions.empty());
const int num_new_particles = ssize(positions);
const double mass_density = config.mass_density();
const double volume_per_particle = total_volume / num_new_particles;
const ConstitutiveModelVariant<T> constitutive_model =
MakeConstitutiveModel<T>(config);
const DeformationGradientDataVariant<T> strain_data =
MakeStrainData<T>(config);
const DeformationGradientDataVariant<T> deformation_gradient_data =
MakeDeformationGradientData<T>(config);
for (int i = 0; i < num_new_particles; ++i) {
m_.push_back(mass_density * volume_per_particle);
x_.emplace_back(positions[i].cast<T>());
v_.emplace_back(Vector3<T>::Zero());
F_.emplace_back(Matrix3<T>::Identity());
F0_.emplace_back(Matrix3<T>::Identity());
C_.emplace_back(Matrix3<T>::Zero());
volume_.push_back(volume_per_particle);
constitutive_models_.push_back(constitutive_model);
tau_volume_.emplace_back(Matrix3<T>::Zero());
in_constraint_.push_back(false);
strain_data_.push_back(strain_data);
}
m_.insert(m_.end(), num_new_particles, mass_density * volume_per_particle);
v_.insert(v_.end(), num_new_particles, Vector3<T>::Zero());
F_.insert(F_.end(), num_new_particles, Matrix3<T>::Identity());
F0_.insert(F0_.end(), num_new_particles, Matrix3<T>::Identity());
C_.insert(C_.end(), num_new_particles, Matrix3<T>::Zero());
volume_.insert(volume_.end(), num_new_particles, volume_per_particle);
constitutive_models_.insert(constitutive_models_.end(), num_new_particles,
constitutive_model);
tau_volume_.insert(tau_volume_.end(), num_new_particles, Matrix3<T>::Zero());
in_constraint_.insert(in_constraint_.end(), num_new_particles, false);
deformation_gradient_data_.insert(deformation_gradient_data_.end(),
num_new_particles,
deformation_gradient_data);
}

} // namespace internal
Expand Down
38 changes: 20 additions & 18 deletions multibody/mpm/particle_data.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,9 +29,9 @@ using DeformationGradientDataVariant =
fem::internal::LinearCorotatedModelData<T>,
fem::internal::LinearConstitutiveModelData<T>>;

/* The collection of all physical attributes we care about for all particles.
All quantities are measured and expressed in the world frame (when
applicable).
/* The collection of all physical attributes we care about for all particles in
a full MPM model. All quantities are measured and expressed in the world frame
(when applicable).
@tparam T The scalar type, can be a double, float, or AutoDiffXd. */
template <typename T>
class ParticleData {
Expand Down Expand Up @@ -59,17 +59,12 @@ class ParticleData {
const std::vector<bool>& in_constraint() const { return in_constraint_; }
const std::vector<Matrix3<T>>& tau_volume() const { return tau_volume_; }

/* Mutators */
std::vector<T>& mutable_m() { return m_; }
/* Mutable accessors */
std::vector<Vector3<T>>& mutable_x() { return x_; }
std::vector<Vector3<T>>& mutable_v() { return v_; }
std::vector<Matrix3<T>>& mutable_F() { return F_; }
std::vector<Matrix3<T>>& mutable_F0() { return F0_; }
std::vector<Matrix3<T>>& mutable_C() { return C_; }
std::vector<T>& mutable_volume() { return volume_; }
std::vector<ConstitutiveModelVariant<T>>& mutable_constitutive_models() {
return constitutive_models_;
}
std::vector<bool>& mutable_in_constraint() { return in_constraint_; }
std::vector<Matrix3<T>>& mutable_tau_volume() { return tau_volume_; }

Expand All @@ -89,7 +84,7 @@ class ParticleData {
@param[in] F The deformation gradients of the particles.
@param[in] F0 The deformation gradients at the previous time step.
@param[out] tau_volume The Kirchhoff stress of each particle scaled by the
volume of the particle.
reference volume of the particle.
@param[in] parallelism Specifies the degree of parallelism to use.
@pre F and F0 have the same size and ordering as this ParticleData.
@pre volume_scaled_stress != nullptr.
Expand Down Expand Up @@ -125,35 +120,42 @@ class ParticleData {
// TODO(xuchenhan-tri): Support Rayleigh damping for MPM.
/* Appends default particle data to `this` ParticleData using the given
parameters. The velocities of the particles are set to zeros. The deformation
gradients are set to identities.
gradients are set to identities. This function can be called repeatedly on a
single ParticleData object to add particles in multiple passes.
@param positions The positions of the new particles in the world frame.
@param total_volume The per particle volume of the new particles is given by
total_volume divided by the number of new particles.
This is using the assumption that the new particles are
evenly distributed across the domain.
@param config Provides the constitutive models and the mass densities
of the new particles.
@pre total_volume > 0. */
void Sample(const std::vector<Vector3<double>>& positions,
double total_volume,
const fem::DeformableBodyConfig<double>& config);
void AddParticles(const std::vector<Vector3<double>>& positions,
double total_volume,
const fem::DeformableBodyConfig<double>& config);

/* Per particle state and data. All of the following fields have the same
size and ordering. */
std::vector<T> m_; // mass
/* State */
std::vector<Vector3<T>> x_; // positions
std::vector<Vector3<T>> v_; // velocity
std::vector<Matrix3<T>> F_; // deformation gradient
std::vector<Matrix3<T>>
F0_; // deformation gradient at the previous time step
std::vector<Matrix3<T>> C_; // affine velocity field
std::vector<T> volume_; // reference volume
/* Constant data. */
std::vector<T> m_; // mass
std::vector<T> volume_; // reference volume
std::vector<ConstitutiveModelVariant<T>> constitutive_models_;
/* State dependent data. */
std::vector<bool> in_constraint_; // whether the particle is participating
// in a constraint
std::vector<Matrix3<T>>
tau_volume_; // Kirchhoff stress scaled by reference volume
mutable std::vector<DeformationGradientDataVariant<T>>
strain_data_; // Deformation gradient dependent data that is used to
// calculate the energy density and its derivatives.
deformation_gradient_data_; // Deformation gradient dependent data that
// is used to calculate the energy density
// and its derivatives.
};

} // namespace internal
Expand Down
Loading

0 comments on commit 33f2d23

Please sign in to comment.