[fem] Refactor ForceDensityField

examples/multibody/deformable/point_source_force_field.cc

@@ -4,7 +4,7 @@ namespace drake {
 namespace examples {
 namespace deformable {
 
-using multibody::ForceDensityField;
+using multibody::ForceDensityFieldBase;
 using multibody::MultibodyPlant;
 using multibody::RigidBody;
 using systems::BasicVector;
@@ -38,7 +38,7 @@ Vector3<double> PointSourceForceField::DoEvaluateAt(
   return magnitude * p_QC_W / p_QC_W.norm();
 }
 
-std::unique_ptr<ForceDensityField<double>> PointSourceForceField::DoClone()
+std::unique_ptr<ForceDensityFieldBase<double>> PointSourceForceField::DoClone()
     const {
   return std::make_unique<PointSourceForceField>(
       *plant_, plant_->get_body(body_), p_BC_, falloff_distance_);

examples/multibody/deformable/point_source_force_field.h

@@ -2,8 +2,8 @@
 
 #include <memory>
 
-#include "drake/multibody/plant/force_density_field.h"
 #include "drake/multibody/plant/multibody_plant.h"
+#include "drake/multibody/tree/force_density_field.h"
 
 namespace drake {
 namespace examples {
@@ -61,7 +61,7 @@ class PointSourceForceField final
   Vector3<double> DoEvaluateAt(const systems::Context<double>& context,
                                const Vector3<double>& p_WQ) const final;
 
-  std::unique_ptr<ForceDensityField<double>> DoClone() const final;
+  std::unique_ptr<ForceDensityFieldBase<double>> DoClone() const final;
 
   void DoDeclareCacheEntries(multibody::MultibodyPlant<double>* plant) final;
 

multibody/fem/BUILD.bazel

@@ -32,6 +32,7 @@ drake_cc_package_library(
         ":fem_solver",
         ":fem_state",
         ":fem_state_system",
+        ":force_density_field_base",
         ":isoparametric_element",
         ":linear_constitutive_model",
         ":linear_constitutive_model_data",
@@ -251,9 +252,9 @@ drake_cc_library(
         "fem_plant_data.h",
     ],
     deps = [
+        ":force_density_field_base",
         "//common:default_scalars",
         "//common:essential",
-        "//multibody/plant:force_density_field",
         "//systems/framework:context",
     ],
 )
@@ -308,6 +309,21 @@ drake_cc_library(
     ],
 )
 
+drake_cc_library(
+    name = "force_density_field_base",
+    srcs = [
+        "force_density_field_base.cc",
+    ],
+    hdrs = [
+        "force_density_field_base.h",
+    ],
+    deps = [
+        "//common:default_scalars",
+        "//common:essential",
+        "//systems/framework:context",
+    ],
+)
+
 drake_cc_library(
     name = "isoparametric_element",
     srcs = [
@@ -715,6 +731,14 @@ drake_cc_googletest(
     ],
 )
 
+drake_cc_googletest(
+    name = "force_density_field_base_test",
+    deps = [
+        ":force_density_field_base",
+        "//multibody/plant",
+    ],
+)
+
 drake_cc_googletest(
     name = "isoparametric_element_test",
     deps = [

multibody/fem/fem_plant_data.h

@@ -3,7 +3,7 @@
 
 #include "drake/common/default_scalars.h"
 #include "drake/common/drake_copyable.h"
-#include "drake/multibody/plant/force_density_field.h"
+#include "drake/multibody/fem/force_density_field_base.h"
 #include "drake/systems/framework/context.h"
 
 namespace drake {
@@ -19,14 +19,14 @@ template <typename T>
 struct FemPlantData {
  public:
   DRAKE_NO_COPY_NO_MOVE_NO_ASSIGN(FemPlantData);
-  FemPlantData(
-      const systems::Context<T>& plant_context_in,
-      const std::vector<const ForceDensityField<T>*>& force_density_fields_in)
+  FemPlantData(const systems::Context<T>& plant_context_in,
+               const std::vector<const ForceDensityFieldBase<T>*>&
+                   force_density_fields_in)
       : plant_context(plant_context_in),
         force_density_fields(force_density_fields_in) {}
 
   const systems::Context<T>& plant_context;
-  std::vector<const ForceDensityField<T>*> force_density_fields;
+  std::vector<const ForceDensityFieldBase<T>*> force_density_fields;
 };
 
 }  // namespace fem

multibody/fem/force_density_field_base.cc

@@ -0,0 +1,13 @@
+#include "drake/multibody/fem/force_density_field_base.h"
+
+namespace drake {
+namespace multibody {
+
+template <typename T>
+ForceDensityFieldBase<T>::~ForceDensityFieldBase() = default;
+
+}  // namespace multibody
+}  // namespace drake
+
+DRAKE_DEFINE_CLASS_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_SCALARS(
+    class ::drake::multibody::ForceDensityFieldBase);

multibody/fem/force_density_field_base.h

@@ -0,0 +1,82 @@
+#pragma once
+
+#include <memory>
+
+#include "drake/common/default_scalars.h"
+#include "drake/common/drake_copyable.h"
+#include "drake/common/eigen_types.h"
+#include "drake/systems/framework/context.h"
+
+namespace drake {
+namespace multibody {
+
+template <typename T>
+class ForceDensityField;
+
+/** (Advanced) Enum for the type of force density in ForceDensityFieldBase. */
+enum class ForceDensityType {
+  /** ForceDensityFieldBase::EvaluateAt() returns the force per unit of
+   _current_ (deformed) configuration volume. */
+  kPerCurrentVolume,
+  /** ForceDensityFieldBase::EvaluateAt() returns the force per unit of
+   _reference_ configuration volume where the reference undeformed configuration
+   is defined by the input mesh provided by the user. */
+  kPerReferenceVolume,
+};
+
+/** The %ForceDensityFieldBase class is an abstract base class that represents a
+ force density field affecting deformable bodies in a MultibodyPlant. The
+ force field is described by the member function EvaluateAt() which takes as
+ input a position in the world frame and returns the force density from the
+ force density field at the given location, with unit [N/m³]. To create a
+ concrete %ForceDensityFieldBase class, inherit from ForceDensityField instead
+ of directly inheriting from %ForceDensityFieldBase.
+ @tparam_default_scalar */
+template <typename T>
+class ForceDensityFieldBase {
+ public:
+  virtual ~ForceDensityFieldBase() = 0;
+
+  /** Evaluates the force density [N/m³] with the given `context` of the
+   owning MultibodyPlant and a position in world, `p_WQ`. */
+  Vector3<T> EvaluateAt(const systems::Context<T>& context,
+                        const Vector3<T>& p_WQ) const {
+    return DoEvaluateAt(context, p_WQ);
+  }
+
+  /** Returns an identical copy of `this` ForceDensityFieldBase. */
+  std::unique_ptr<ForceDensityFieldBase<T>> Clone() const { return DoClone(); }
+
+  /* (Advanced) Returns the force density type of `this` %ForceDensityField.
+   */
+  ForceDensityType density_type() const { return density_type_; }
+
+ protected:
+  DRAKE_DEFAULT_COPY_AND_MOVE_AND_ASSIGN(ForceDensityFieldBase);
+
+  /** Derived classes must override this function to provide a threadsafe
+   implemention to the NVI EvaluateAt(). */
+  virtual Vector3<T> DoEvaluateAt(const systems::Context<T>& context,
+                                  const Vector3<T>& p_WQ) const = 0;
+
+  /** Derived classes must override this function to implement the NVI
+   Clone(). */
+  virtual std::unique_ptr<ForceDensityFieldBase<T>> DoClone() const = 0;
+
+ private:
+  friend class ForceDensityField<T>;
+
+  /* Private constructor exposed only to ForceDensityField. This prevents
+   users from creating concrete force densities that directly inherit from
+   ForceDensityFieldBase. */
+  explicit ForceDensityFieldBase(ForceDensityType density_type)
+      : density_type_(density_type) {}
+
+  ForceDensityType density_type_;
+};
+
+}  // namespace multibody
+}  // namespace drake
+
+DRAKE_DECLARE_CLASS_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_SCALARS(
+    class ::drake::multibody::ForceDensityFieldBase);

multibody/fem/test/dummy_element.h

@@ -151,7 +151,7 @@ class DummyElement final : public FemElement<DummyElement<is_linear>> {
                                  EigenPtr<Vector<T, kNumDofs>> result) const {
     for (int i = 0; i < this->num_nodes; ++i) {
       const auto node_q = data.element_q.template segment<3>(3 * i);
-      for (const multibody::ForceDensityField<T>* force_density :
+      for (const multibody::ForceDensityFieldBase<T>* force_density :
            plant_data.force_density_fields) {
         result->template segment<3>(3 * i) +=
             scale * force_density->EvaluateAt(plant_data.plant_context, node_q);

multibody/fem/test/fem_element_test.cc

@@ -4,6 +4,7 @@
 
 #include "drake/multibody/fem/test/dummy_element.h"
 #include "drake/multibody/plant/multibody_plant.h"
+#include "drake/multibody/tree/force_density_field.h"
 
 namespace drake {
 namespace multibody {

multibody/fem/test/fem_model_test.cc

@@ -11,6 +11,7 @@
 #include "drake/multibody/fem/volumetric_element.h"
 #include "drake/multibody/fem/volumetric_model.h"
 #include "drake/multibody/plant/multibody_plant.h"
+#include "drake/multibody/tree/force_density_field.h"
 
 namespace drake {
 namespace multibody {
@@ -134,7 +135,7 @@ GTEST_TEST(FemModelTest, CalcResidualWithContextDependentExternalForce) {
       return context.get_time() * unit_vector_;
     };
 
-    std::unique_ptr<ForceDensityField<double>> DoClone() const final {
+    std::unique_ptr<ForceDensityFieldBase<double>> DoClone() const final {
       return std::make_unique<TimeScaledForceDensityField>(*this);
     }
 

multibody/fem/test/force_density_field_base_test.cc

@@ -0,0 +1,113 @@
+#include "drake/multibody/fem/force_density_field_base.h"
+
+#include <gtest/gtest.h>
+
+#include "drake/multibody/plant/multibody_plant.h"
+#include "drake/multibody/tree/force_density_field.h"
+
+namespace drake {
+namespace multibody {
+namespace fem {
+namespace {
+
+using drake::systems::BasicVector;
+using Eigen::Vector3d;
+
+/* A concrete force density field for testing. */
+class ConstantForceDensityField final : public ForceDensityField<double> {
+ public:
+  /* Constructs an force density field that has the functional form given by
+   input `field` which is then scaled by a scalar value via input port. */
+  explicit ConstantForceDensityField(
+      std::function<Vector3<double>(const Vector3<double>&)> field)
+      : field_(std::move(field)) {}
+
+  /* Gets the double-valued input port -- the input port value scales the
+   force density multiplicatively. */
+  const systems::InputPort<double>& get_input_port() const {
+    return parent_system_or_throw().get_input_port(scale_port_index_);
+  }
+
+ private:
+  DRAKE_DEFAULT_COPY_AND_MOVE_AND_ASSIGN(ConstantForceDensityField);
+
+  Vector3<double> DoEvaluateAt(const systems::Context<double>& context,
+                               const Vector3<double>& p_WQ) const final {
+    return get_input_port().Eval(context)(0) * field_(p_WQ);
+  };
+
+  std::unique_ptr<ForceDensityFieldBase<double>> DoClone() const final {
+    return std::unique_ptr<ConstantForceDensityField>(
+        new ConstantForceDensityField(*this));
+  }
+
+  void DoDeclareInputPorts(multibody::MultibodyPlant<double>* plant) final {
+    scale_port_index_ =
+        this->DeclareVectorInputPort(plant, "on/off signal for the force field",
+                                     BasicVector<double>(1.0))
+            .get_index();
+  }
+
+  std::function<Vector3<double>(const Vector3<double>&)> field_;
+  systems::InputPortIndex scale_port_index_;
+};
+
+class ForceDensityFieldTest : public ::testing::Test {
+ protected:
+  void SetUp() override {
+    /* Define an arbtrary force density field. */
+    auto force_field = [](const Vector3d& x) {
+      return Vector3d(std::sin(x[0]), std::cos(x[1]), x[2]);
+    };
+    dut_ = std::make_unique<ConstantForceDensityField>(force_field);
+    dut_->DeclareSystemResources(&plant_);
+    plant_.Finalize();
+  }
+
+  std::unique_ptr<ConstantForceDensityField> dut_;
+  MultibodyPlant<double> plant_{0.01};
+};
+
+TEST_F(ForceDensityFieldTest, EvaluateAt) {
+  auto plant_context = plant_.CreateDefaultContext();
+  EXPECT_EQ(dut_->density_type(), ForceDensityType::kPerCurrentVolume);
+  const double scale = 2.71;
+  dut_->get_input_port().FixValue(plant_context.get(), Vector1d(scale));
+  const Vector3d p_WQ(1, 2, 3);
+  const Vector3d expected_force_density =
+      scale * Vector3d(std::sin(p_WQ[0]), std::cos(p_WQ[1]), p_WQ[2]);
+  EXPECT_EQ(dut_->EvaluateAt(*plant_context, p_WQ), expected_force_density);
+}
+
+TEST_F(ForceDensityFieldTest, Clone) {
+  auto clone = dut_->Clone();
+  ASSERT_NE(clone, nullptr);
+  EXPECT_EQ(clone->density_type(), dut_->density_type());
+
+  auto concrete_clone = dynamic_cast<ConstantForceDensityField*>(clone.get());
+  ASSERT_NE(concrete_clone, nullptr);
+  MultibodyPlant<double> plant_clone(0.01);
+  concrete_clone->DeclareSystemResources(&plant_clone);
+  plant_clone.Finalize();
+  auto plant_context = plant_.CreateDefaultContext();
+  auto clone_context = plant_clone.CreateDefaultContext();
+
+  const double scale = 2.71;
+  dut_->get_input_port().FixValue(plant_context.get(), Vector1d(scale));
+  concrete_clone->get_input_port().FixValue(clone_context.get(),
+                                            Vector1d(scale));
+  for (int i = 0; i < 3; ++i) {
+    for (int j = 0; j < 3; ++j) {
+      for (int k = 0; k < 3; ++k) {
+        const Vector3d p_WQ(i, j, k);
+        EXPECT_EQ(dut_->EvaluateAt(*plant_context, p_WQ),
+                  clone->EvaluateAt(*clone_context, p_WQ));
+      }
+    }
+  }
+}
+
+}  // namespace
+}  // namespace fem
+}  // namespace multibody
+}  // namespace drake

multibody/fem/test/volumetric_element_test.cc

@@ -11,6 +11,7 @@
 #include "drake/multibody/fem/linear_simplex_element.h"
 #include "drake/multibody/fem/simplex_gaussian_quadrature.h"
 #include "drake/multibody/plant/multibody_plant.h"
+#include "drake/multibody/tree/force_density_field.h"
 
 namespace drake {
 namespace multibody {
@@ -429,7 +430,7 @@ TEST_F(VolumetricElementTest, PerCurrentVolumeExternalForce) {
       return f_;
     };
 
-    std::unique_ptr<ForceDensityField<AD>> DoClone() const final {
+    std::unique_ptr<ForceDensityFieldBase<AD>> DoClone() const final {
       return std::make_unique<ConstantForceDensityField>(*this);
     }
 

multibody/fem/volumetric_element.h

@@ -441,7 +441,7 @@ class VolumetricElement
       const Matrix3<T>& deformation_gradient =
           data.deformation_gradient_data[q].deformation_gradient();
       Vector3<T> scaled_force = Vector3<T>::Zero();
-      for (const multibody::ForceDensityField<T>* force_density :
+      for (const multibody::ForceDensityFieldBase<T>* force_density :
            plant_data.force_density_fields) {
         DRAKE_ASSERT(force_density != nullptr);
         const T change_of_volume =