Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: 42889b664f84
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: fadba09a8378
Choose a head ref
  • 15 commits
  • 12 files changed
  • 1 contributor

Commits on Dec 1, 2019

  1. Copy the full SHA
    001ffb0 View commit details
  2. Copy the full SHA
    1c77ec9 View commit details

Commits on Dec 2, 2019

  1. Copy the full SHA
    a2cf9c2 View commit details
  2. More tensor operations.

    pleroy committed Dec 2, 2019
    Copy the full SHA
    d9f1c68 View commit details
  3. Inertia and angular momentum.

    pleroy committed Dec 2, 2019
    Copy the full SHA
    22d070d View commit details

Commits on Dec 7, 2019

  1. Copy the full SHA
    82be288 View commit details
  2. Copy the full SHA
    be0b05a View commit details
  3. Copy the full SHA
    20d0f9f View commit details

Commits on Dec 8, 2019

  1. Compilation errors.

    pleroy committed Dec 8, 2019
    Copy the full SHA
    eb0ce87 View commit details
  2. Copy the full SHA
    46f3c36 View commit details

Commits on Dec 12, 2019

  1. After egg's review.

    pleroy committed Dec 12, 2019
    Copy the full SHA
    a13ca10 View commit details

Commits on Dec 13, 2019

  1. Another round of review.

    pleroy committed Dec 13, 2019
    Copy the full SHA
    6a835a0 View commit details
  2. Compute in the RigidPileUp.

    pleroy committed Dec 13, 2019
    Copy the full SHA
    1376d4c View commit details
  3. Merge.

    pleroy committed Dec 13, 2019
    Copy the full SHA
    bf14341 View commit details

Commits on Dec 14, 2019

  1. Merge pull request #2397 from pleroy/PileUpInertia

    Propagate the dynamical elements from the parts to the pile-up
    pleroy authored Dec 14, 2019
    Copy the full SHA
    fadba09 View commit details
3 changes: 2 additions & 1 deletion geometry/symmetric_bilinear_form.hpp
Original file line number Diff line number Diff line change
@@ -196,7 +196,7 @@ SymmetricBilinearForm<Product<LScalar, RScalar>, Frame> SymmetricProduct(
Vector<RScalar, Frame> const& right);

// Symmetric bilinear forms act on bivectors through this function.
// |Anticommutator(F, B)| is (tr(F)𝟙 - F)B in ℝ³ representation . In matrix
// |Anticommutator(F, B)| is (tr(F)𝟙 - F)B in ℝ³ representation. In matrix
// representation it is FB + BF = {F, B}.
template<typename LScalar, typename RScalar, typename Frame>
Bivector<Product<LScalar, RScalar>, Frame> Anticommutator(
@@ -219,6 +219,7 @@ std::ostream& operator<<(std::ostream& out,

} // namespace internal_symmetric_bilinear_form

using internal_symmetric_bilinear_form::Anticommutator;
using internal_symmetric_bilinear_form::InnerProductForm;
using internal_symmetric_bilinear_form::SymmetricBilinearForm;
using internal_symmetric_bilinear_form::SymmetricProduct;
10 changes: 10 additions & 0 deletions ksp_plugin/frames.hpp
Original file line number Diff line number Diff line change
@@ -121,6 +121,15 @@ using MainBodyCentred = Frame<serialization::Frame::PluginTag,
serialization::Frame::MAIN_BODY_CENTRED,
NonInertial>;

// The |PileUp| is seen as a rigid body; the degrees of freedom of the parts in
// the frame of that body can be set, however their motion is not integrated;
// this is simply applied as an offset from the rigid body motion of the
// |PileUp|. The origin of |RigidPileUp| is the centre of mass of the pile up.
// Its axes are those of Barycentric.
using RigidPileUp = Frame<serialization::Frame::PluginTag,
serialization::Frame::RIGID_PILE_UP,
NonInertial>;

// Convenient instances of types from |physics| for the above frames.
using NavigationFrame = DynamicFrame<Barycentric, Navigation>;
using NavigationManœuvre = Manœuvre<Barycentric, Navigation>;
@@ -145,6 +154,7 @@ using internal_frames::Navigation;
using internal_frames::NavigationFrame;
using internal_frames::NavigationManœuvre;
using internal_frames::RigidPart;
using internal_frames::RigidPileUp;
using internal_frames::World;
using internal_frames::WorldSun;
using internal_frames::sun_looking_glass;
9 changes: 2 additions & 7 deletions ksp_plugin/part_subsets.cpp
Original file line number Diff line number Diff line change
@@ -20,9 +20,7 @@ using physics::Ephemeris;

namespace base {

Subset<Part>::Properties::Properties(not_null<ksp_plugin::Part*> const part)
: total_mass_(part->inertia_tensor().mass()),
total_intrinsic_force_(part->intrinsic_force()) {
Subset<Part>::Properties::Properties(not_null<ksp_plugin::Part*> const part) {
if (part->is_piled_up()) {
missing_ = part->containing_pile_up()->parts().size() - 1;
}
@@ -45,8 +43,6 @@ void Subset<Part>::Properties::MergeWith(Properties& other) {
}
}
parts_.splice(parts_.end(), other.parts_);
total_mass_ += other.total_mass_;
total_intrinsic_force_ += other.total_intrinsic_force_;
grounded_ |= other.grounded_;
}

@@ -71,8 +67,7 @@ void Subset<Part>::Properties::Collect(
collected_ = true;
if (EqualsExistingPileUp()) {
PileUp& pile_up = *parts_.front()->containing_pile_up();
pile_up.set_mass(total_mass_);
pile_up.set_intrinsic_force(total_intrinsic_force_);
pile_up.RecomputeFromParts();
} else {
if (StrictSubsetOfExistingPileUp()) {
for (auto const part : parts_) {
5 changes: 0 additions & 5 deletions ksp_plugin/part_subsets.hpp
Original file line number Diff line number Diff line change
@@ -93,11 +93,6 @@ class Subset<ksp_plugin::Part>::Properties final {
int missing_;
// The list of parts in this subset.
std::list<not_null<ksp_plugin::Part*>> parts_;
// The sum of the masses of the |parts_|.
quantities::Mass total_mass_;
// The sum of the |intrinsic_force|s on the |parts_|.
geometry::Vector<quantities::Force, ksp_plugin::Barycentric>
total_intrinsic_force_;
// Whether the subset touches the ground.
bool grounded_ = false;
};
144 changes: 98 additions & 46 deletions ksp_plugin/pile_up.cpp
Original file line number Diff line number Diff line change
@@ -6,11 +6,13 @@
#include <map>

#include "base/map_util.hpp"
#include "geometry/grassmann.hpp"
#include "geometry/identity.hpp"
#include "geometry/named_quantities.hpp"
#include "ksp_plugin/integrators.hpp"
#include "ksp_plugin/part.hpp"
#include "physics/rigid_motion.hpp"
#include "quantities/named_quantities.hpp"

namespace principia {
namespace ksp_plugin {
@@ -21,13 +23,18 @@ using base::FindOrDie;
using base::make_not_null_unique;
using geometry::AngularVelocity;
using geometry::BarycentreCalculator;
using geometry::Bivector;
using geometry::Identity;
using geometry::OrthogonalMap;
using geometry::Position;
using geometry::RigidTransformation;
using geometry::Velocity;
using geometry::Wedge;
using physics::Anticommutator;
using physics::DegreesOfFreedom;
using physics::RigidMotion;
using quantities::AngularMomentum;
using quantities::si::Radian;
using ::std::placeholders::_1;
using ::std::placeholders::_2;
using ::std::placeholders::_3;
@@ -48,15 +55,7 @@ PileUp::PileUp(
history_(make_not_null_unique<DiscreteTrajectory<Barycentric>>()),
deletion_callback_(std::move(deletion_callback)) {
LOG(INFO) << "Constructing pile up at " << this;
BarycentreCalculator<DegreesOfFreedom<Barycentric>, Mass> calculator;
Vector<Force, Barycentric> total_intrinsic_force;
for (not_null<Part*> const part : parts_) {
total_intrinsic_force += part->intrinsic_force();
calculator.Add(part->degrees_of_freedom(), part->inertia_tensor().mass());
}
mass_ = calculator.weight();
intrinsic_force_ = total_intrinsic_force;
DegreesOfFreedom<Barycentric> const barycentre = calculator.Get();
DegreesOfFreedom<Barycentric> const barycentre = RecomputeFromParts(parts_);
history_->Append(t, barycentre);

RigidMotion<Barycentric, RigidPileUp> const barycentric_to_pile_up{
@@ -81,15 +80,6 @@ PileUp::~PileUp() {
}
}

void PileUp::set_mass(Mass const& mass) {
mass_ = mass;
}

void PileUp::set_intrinsic_force(
Vector<Force, Barycentric> const& intrinsic_force) {
intrinsic_force_ = intrinsic_force;
}

std::list<not_null<Part*>> const& PileUp::parts() const {
return parts_;
}
@@ -103,28 +93,6 @@ void PileUp::SetPartApparentDegreesOfFreedom(
<< degrees_of_freedom;
}

void PileUp::NudgeParts() const {
auto const actual_centre_of_mass =
psychohistory_->back().degrees_of_freedom;

RigidMotion<Barycentric, RigidPileUp> const barycentric_to_pile_up{
RigidTransformation<Barycentric, RigidPileUp>{
actual_centre_of_mass.position(),
RigidPileUp::origin,
Identity<Barycentric, RigidPileUp>().Forget()},
AngularVelocity<Barycentric>(),
actual_centre_of_mass.velocity()};
auto const pile_up_to_barycentric = barycentric_to_pile_up.Inverse();
for (not_null<Part*> const part : parts_) {
DegreesOfFreedom<Barycentric> const actual_part_degrees_of_freedom =
pile_up_to_barycentric(
FindOrDie(actual_part_degrees_of_freedom_, part));
part->set_rigid_motion(
RigidMotion<RigidPart, Barycentric>::MakeNonRotatingMotion(
actual_part_degrees_of_freedom));
}
}

Status PileUp::DeformAndAdvanceTime(Instant const& t) {
absl::MutexLock l(lock_.get());
Status status;
@@ -136,12 +104,14 @@ Status PileUp::DeformAndAdvanceTime(Instant const& t) {
return status;
}

void PileUp::RecomputeFromParts() {
RecomputeFromParts(parts_);
}

void PileUp::WriteToMessage(not_null<serialization::PileUp*> message) const {
for (not_null<Part*> const part : parts_) {
message->add_part_id(part->part_id());
}
mass_.WriteToMessage(message->mutable_mass());
intrinsic_force_.WriteToMessage(message->mutable_intrinsic_force());
history_->WriteToMessage(message->mutable_history(),
/*forks=*/{psychohistory_});
for (auto const& pair : actual_part_degrees_of_freedom_) {
@@ -233,9 +203,6 @@ not_null<std::unique_ptr<PileUp>> PileUp::ReadFromMessage(
std::move(deletion_callback)));
}

pile_up->mass_ = Mass::ReadFromMessage(message.mass());
pile_up->intrinsic_force_ =
Vector<Force, Barycentric>::ReadFromMessage(message.intrinsic_force());
for (auto const& pair : message.actual_part_degrees_of_freedom()) {
std::uint32_t const part_id = pair.first;
serialization::Pair const& degrees_of_freedom = pair.second;
@@ -250,6 +217,7 @@ not_null<std::unique_ptr<PileUp>> PileUp::ReadFromMessage(
part_id_to_part(part_id),
DegreesOfFreedom<ApparentBubble>::ReadFromMessage(degrees_of_freedom));
}
pile_up->RecomputeFromParts();
return check_not_null(std::move(pile_up));
}

@@ -362,7 +330,7 @@ Status PileUp::AdvanceTime(Instant const& t) {
}
history_->DeleteFork(psychohistory_);

auto const a = intrinsic_force_ / mass_;
auto const a = intrinsic_force_ / inertia_tensor_.mass();
// NOTE(phl): |a| used to be captured by copy below, which is the logical
// thing to do. However, since it contains an |R3Element|, it must be
// aligned on a 16-byte boundary. Unfortunately, VS2015 gets confused and
@@ -397,6 +365,28 @@ Status PileUp::AdvanceTime(Instant const& t) {
return status;
}

void PileUp::NudgeParts() const {
auto const actual_centre_of_mass =
psychohistory_->back().degrees_of_freedom;

RigidMotion<Barycentric, RigidPileUp> const barycentric_to_pile_up{
RigidTransformation<Barycentric, RigidPileUp>{
actual_centre_of_mass.position(),
RigidPileUp::origin,
Identity<Barycentric, RigidPileUp>().Forget()},
AngularVelocity<Barycentric>(),
actual_centre_of_mass.velocity()};
auto const pile_up_to_barycentric = barycentric_to_pile_up.Inverse();
for (not_null<Part*> const part : parts_) {
DegreesOfFreedom<Barycentric> const actual_part_degrees_of_freedom =
pile_up_to_barycentric(
FindOrDie(actual_part_degrees_of_freedom_, part));
part->set_rigid_motion(
RigidMotion<RigidPart, Barycentric>::MakeNonRotatingMotion(
actual_part_degrees_of_freedom));
}
}

template<PileUp::AppendToPartTrajectory append_to_part_trajectory>
void PileUp::AppendToPart(DiscreteTrajectory<Barycentric>::Iterator it) const {
auto const& pile_up_dof = it->degrees_of_freedom;
@@ -416,6 +406,68 @@ void PileUp::AppendToPart(DiscreteTrajectory<Barycentric>::Iterator it) const {
}
}

DegreesOfFreedom<Barycentric> PileUp::RecomputeFromParts(
std::list<not_null<Part*>> const& parts) {
// First compute the overall mass and centre of mass of the pile-up. Also
// compute the overall force applied to it.
// TODO(phl): We assume that forces are applied at the centre of mass of the
// pile-up, but they are really applied at some unknown point of the parts, so
// this introduces a torque.
Mass pile_up_mass;
Vector<Force, Barycentric> pile_up_intrinsic_force;
BarycentreCalculator<DegreesOfFreedom<Barycentric>, Mass> calculator;
for (not_null<Part*> const part : parts) {
pile_up_intrinsic_force += part->intrinsic_force();
calculator.Add(part->degrees_of_freedom(), part->inertia_tensor().mass());
}
pile_up_mass = calculator.weight();
DegreesOfFreedom<Barycentric> const pile_up_barycentre = calculator.Get();

// Then compute the inertia tensor and the angular momentum of the pile-up.
InertiaTensor<RigidPileUp> pile_up_inertia_tensor;
Bivector<AngularMomentum, RigidPileUp> pile_up_angular_momentum;
RigidMotion<RigidPileUp, Barycentric> const pile_up_to_barycentric(
RigidTransformation<RigidPileUp, Barycentric>(
RigidPileUp::origin,
pile_up_barycentre.position(),
Identity<RigidPileUp, Barycentric>().Forget()),
AngularVelocity<Barycentric>{},
pile_up_barycentre.velocity());
RigidMotion<Barycentric, RigidPileUp> const barycentric_to_pile_up =
pile_up_to_barycentric.Inverse();
for (not_null<Part*> const part : parts) {
RigidMotion<RigidPart, Barycentric> const& part_to_barycentric =
part->rigid_motion();
RigidMotion<RigidPart, RigidPileUp> const part_to_pile_up =
barycentric_to_pile_up * part_to_barycentric;

// NOTE(phl): The following call reads like the tensor "transforms" the
// rigid motion. Improve this API.
InertiaTensor<RigidPileUp> const& part_inertia_tensor =
part->inertia_tensor().Transform(
part_to_pile_up.rigid_transformation());
pile_up_inertia_tensor += part_inertia_tensor;

Bivector<AngularMomentum, RigidPileUp> const part_angular_momentum =
Anticommutator(
part_inertia_tensor,
part_to_pile_up.Inverse().angular_velocity_of_to_frame());
DegreesOfFreedom<RigidPileUp> const part_degrees_of_freedom =
part_to_pile_up({RigidPart::origin, Velocity<RigidPart>{}});
pile_up_angular_momentum +=
part_angular_momentum +
Wedge(part_degrees_of_freedom.position() - RigidPileUp::origin,
part->inertia_tensor().mass() *
part_degrees_of_freedom.velocity()) *
Radian;
}

angular_momentum_ = pile_up_angular_momentum;
inertia_tensor_ = pile_up_inertia_tensor;
intrinsic_force_ = pile_up_intrinsic_force;
return pile_up_barycentre;
}

PileUpFuture::PileUpFuture(not_null<PileUp const*> const pile_up,
std::future<Status> future)
: pile_up(pile_up),
Loading