[fem] Refactor ForceDensityField

The goal of the commit is the remove the dependency of the fem library
from the tree library so that later we create deformable multibody
element (in the tree library) that depends on fem.

Author: xuchenhan-tri

examples/multibody/deformable/point_source_force_field.h

@@ -2,7 +2,7 @@
 
 #include <memory>
 
-#include "drake/multibody/plant/force_density_field.h"
+#include "drake/multibody/tree/force_density_field_impl.h"
 #include "drake/multibody/plant/multibody_plant.h"
 
 namespace drake {
@@ -16,7 +16,7 @@ namespace deformable {
  a double-valued input port (see maximum_force_density_input_port()). If the
  port is unconnected, it reads as zero. */
 class PointSourceForceField final
-    : public multibody::ForceDensityField<double> {
+    : public multibody::ForceDensityFieldImpl<double> {
  public:
   DRAKE_NO_COPY_NO_MOVE_NO_ASSIGN(PointSourceForceField);
 

multibody/fem/BUILD.bazel

@@ -32,6 +32,7 @@ drake_cc_package_library(
         ":fem_solver",
         ":fem_state",
         ":fem_state_system",
+        ":force_density_field",
         ":isoparametric_element",
         ":linear_constitutive_model",
         ":linear_constitutive_model_data",
@@ -251,9 +252,9 @@ drake_cc_library(
         "fem_plant_data.h",
     ],
     deps = [
+        ":force_density_field",
         "//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",
+    srcs = [
+        "force_density_field.cc",
+    ],
+    hdrs = [
+        "force_density_field.h",
+    ],
+    deps = [
+        "//common:default_scalars",
+        "//common:essential",
+        "//systems/framework:context",
+    ],
+)
+
 drake_cc_library(
     name = "isoparametric_element",
     srcs = [

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.h"
 #include "drake/systems/framework/context.h"
 
 namespace drake {

multibody/fem/force_density_field.cc

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

multibody/fem/force_density_field.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 ForceDensityFieldImpl;
+
+/** (Advanced) Enum for the type of force density in ForceDensityField. */
+enum class ForceDensityType {
+  /** ForceDensityField::EvaluateAt() returns the force per unit of _current_
+   (deformed) configuration volume. */
+  kPerCurrentVolume,
+  /** ForceDensityField::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 %ForceDensityField 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 %ForceDensityField class, inherit from ForceDensityFieldImpl instead
+ of directly inheriting from %ForceDensityField.
+ @tparam_default_scalar */
+template <typename T>
+class ForceDensityField {
+ public:
+  virtual ~ForceDensityField() = 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` ForceDensityField. */
+  std::unique_ptr<ForceDensityField<T>> Clone() const { return DoClone(); }
+
+  /* (Advanced) Returns the force density type of `this` %ForceDensityFieldImpl.
+   */
+  ForceDensityType density_type() const { return density_type_; }
+
+ protected:
+  DRAKE_DEFAULT_COPY_AND_MOVE_AND_ASSIGN(ForceDensityField);
+
+  /** 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<ForceDensityField<T>> DoClone() const = 0;
+
+ private:
+  friend class ForceDensityFieldImpl<T>;
+
+  /* Private constructor exposed only to ForceDensityFieldImpl. This prevents
+   users from creating concrete force densities that directly inherit from
+   ForceDensityField. */
+  explicit ForceDensityField(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::ForceDensityField);

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_impl.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_impl.h"
 
 namespace drake {
 namespace multibody {
@@ -118,7 +119,8 @@ GTEST_TEST(FemModelTest, CalcResidualWithContextDependentExternalForce) {
   unique_ptr<FemState<double>> fem_state = model.MakeFemState();
 
   /* A force field where the magnitude of the force density depends on time. */
-  class TimeScaledForceDensityField final : public ForceDensityField<double> {
+  class TimeScaledForceDensityField final
+      : public ForceDensityFieldImpl<double> {
    public:
     DRAKE_DEFAULT_COPY_AND_MOVE_AND_ASSIGN(TimeScaledForceDensityField);
 

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_impl.h"
 
 namespace drake {
 namespace multibody {
@@ -419,7 +420,7 @@ TEST_F(VolumetricElementTest, PerCurrentVolumeExternalForce) {
 
   const Vector3<AD> force_per_current_volume(4, 5, 6);
   /* A constant per current volume force density field. */
-  class ConstantForceDensityField final : public ForceDensityField<AD> {
+  class ConstantForceDensityField final : public ForceDensityFieldImpl<AD> {
    public:
     explicit ConstantForceDensityField(const Vector3<AD>& f) : f_(f) {}
 

multibody/plant/BUILD.bazel

@@ -26,7 +26,6 @@ drake_cc_package_library(
         ":distance_constraint_params",
         ":externally_applied_spatial_force",
         ":externally_applied_spatial_force_multiplexer",
-        ":force_density_field",
         ":internal_geometry_names",
         ":multibody_plant_config",
         ":multibody_plant_config_functions",
@@ -92,7 +91,7 @@ drake_cc_library(
         "discrete_step_memory.cc",
         "discrete_update_manager.cc",
         "dummy_physical_model.cc",
-        "force_density_field_declare_system_resources.cc",
+        "force_density_field_impl_declare_system_resources.cc",
         "geometry_contact_data.cc",
         "make_discrete_update_manager.cc",
         "multibody_element_getter.cc",
@@ -238,22 +237,6 @@ drake_cc_library(
     ],
 )
 
-drake_cc_library(
-    name = "force_density_field",
-    srcs = [
-        "force_density_field.cc",
-    ],
-    hdrs = [
-        "force_density_field.h",
-    ],
-    deps = [
-        "//common:default_scalars",
-        "//common:essential",
-        "//multibody/tree",
-        "//systems/framework:context",
-    ],
-)
-
 drake_cc_library(
     name = "internal_geometry_names",
     srcs = ["internal_geometry_names.cc"],

multibody/plant/deformable_driver.cc

@@ -17,7 +17,6 @@
 #include "drake/multibody/fem/fem_model.h"
 #include "drake/multibody/plant/contact_properties.h"
 #include "drake/multibody/plant/discrete_update_manager.h"
-#include "drake/multibody/plant/force_density_field.h"
 #include "drake/multibody/plant/hydroelastic_quadrature_point_data.h"
 #include "drake/multibody/plant/multibody_plant.h"
 #include "drake/systems/framework/context.h"

multibody/plant/deformable_model.cc

@@ -14,6 +14,7 @@
 #include "drake/multibody/fem/velocity_newmark_scheme.h"
 #include "drake/multibody/fem/volumetric_model.h"
 #include "drake/multibody/plant/multibody_plant.h"
+#include "drake/multibody/tree/force_density_field_impl.h"
 
 namespace drake {
 namespace multibody {
@@ -583,7 +584,10 @@ void DeformableModel<T>::DoDeclareSystemResources() {
    them. */
   for (std::unique_ptr<ForceDensityField<T>>& force_density :
        force_densities_) {
-    force_density->DeclareSystemResources(this->mutable_plant());
+    /* We know that the static cast is safe because all concrete force density
+     field is derived from ForceDensityFieldImpl. */
+    static_cast<ForceDensityFieldImpl<T>*>(force_density.get())
+        ->DeclareSystemResources(this->mutable_plant());
   }
 }
 

multibody/plant/deformable_model.h

@@ -13,9 +13,9 @@
 #include "drake/multibody/fem/deformable_body_config.h"
 #include "drake/multibody/fem/discrete_time_integrator.h"
 #include "drake/multibody/fem/fem_model.h"
+#include "drake/multibody/fem/force_density_field.h"
 #include "drake/multibody/plant/constraint_specs.h"
 #include "drake/multibody/plant/deformable_ids.h"
-#include "drake/multibody/plant/force_density_field.h"
 #include "drake/multibody/plant/physical_model.h"
 #include "drake/multibody/tree/rigid_body.h"
 

multibody/plant/force_density_field_declare_system_resources.cc renamed to multibody/plant/force_density_field_impl_declare_system_resources.cc

@@ -1,6 +1,6 @@
 #include "drake/common/default_scalars.h"
-#include "drake/multibody/plant/force_density_field.h"
 #include "drake/multibody/plant/multibody_plant.h"
+#include "drake/multibody/tree/force_density_field_impl.h"
 
 namespace drake {
 namespace multibody {
@@ -10,7 +10,7 @@ namespace multibody {
 // isn't a no-op.
 template <typename T>
 template <typename>
-void ForceDensityField<T>::DeclareSystemResources(
+void ForceDensityFieldImpl<T>::DeclareSystemResources(
     internal::MultibodyTreeSystem<T>* tree_system) {
   DRAKE_DEMAND(tree_system != nullptr);
   /* `this` force field isn't already associated with a system. */
@@ -29,7 +29,7 @@ void ForceDensityField<T>::DeclareSystemResources(
 
 // clang-format off
 DRAKE_DEFINE_FUNCTION_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_SCALARS((
-    &ForceDensityField<T>::template DeclareSystemResources<void>
+    &ForceDensityFieldImpl<T>::template DeclareSystemResources<void>
 ));
 // clang-format on
 

multibody/plant/test/deformable_model_test.cc

@@ -6,6 +6,7 @@
 #include "drake/common/test_utilities/expect_throws_message.h"
 #include "drake/multibody/plant/multibody_plant.h"
 #include "drake/multibody/plant/multibody_plant_config_functions.h"
+#include "drake/multibody/tree/force_density_field_impl.h"
 #include "drake/systems/framework/diagram_builder.h"
 
 namespace drake {
@@ -348,7 +349,7 @@ TEST_F(DeformableModelTest, AddFixedConstraint) {
 
 TEST_F(DeformableModelTest, ExternalForces) {
   /* A user defined force density field. */
-  class ConstantForceDensityField final : public ForceDensityField<double> {
+  class ConstantForceDensityField final : public ForceDensityFieldImpl<double> {
    public:
     DRAKE_DEFAULT_COPY_AND_MOVE_AND_ASSIGN(ConstantForceDensityField);
 

multibody/tree/BUILD.bazel

@@ -92,6 +92,7 @@ drake_cc_library(
         "door_hinge.cc",
         "element_collection.cc",
         "fixed_offset_frame.cc",
+        "force_density_field_impl.cc",
         "force_element.cc",
         "frame.cc",
         "joint.cc",
@@ -137,6 +138,7 @@ drake_cc_library(
         "door_hinge.h",
         "element_collection.h",
         "fixed_offset_frame.h",
+        "force_density_field_impl.h",
         "force_element.h",
         "frame.h",
         "joint.h",
@@ -189,6 +191,7 @@ drake_cc_library(
         "//common:unused",
         "//common/trajectories:piecewise_constant_curvature_trajectory",
         "//math:geometric_transform",
+        "//multibody/fem:force_density_field",
         "//multibody/topology",
         "//systems/framework:leaf_system",
     ],

multibody/plant/force_density_field.cc renamed to multibody/tree/force_density_field_impl.cc

@@ -1,13 +1,13 @@
-#include "drake/multibody/plant/force_density_field.h"
+#include "drake/multibody/tree/force_density_field_impl.h"
 
 namespace drake {
 namespace multibody {
 
 template <typename T>
-ForceDensityField<T>::~ForceDensityField() = default;
+ForceDensityFieldImpl<T>::~ForceDensityFieldImpl() = default;
 
 template <typename T>
-systems::CacheEntry& ForceDensityField<T>::DeclareCacheEntry(
+systems::CacheEntry& ForceDensityFieldImpl<T>::DeclareCacheEntry(
     internal::MultibodyTreeSystem<T>* plant, std::string description,
     systems::ValueProducer value_producer,
     std::set<systems::DependencyTicket> prerequisites_of_calc) {
@@ -17,15 +17,15 @@ systems::CacheEntry& ForceDensityField<T>::DeclareCacheEntry(
 }
 
 template <typename T>
-systems::InputPort<T>& ForceDensityField<T>::DeclareAbstractInputPort(
+systems::InputPort<T>& ForceDensityFieldImpl<T>::DeclareAbstractInputPort(
     internal::MultibodyTreeSystem<T>* plant, std::string name,
     const AbstractValue& model_value) {
   return internal::MultibodyTreeSystemElementAttorney<
       T>::DeclareAbstractInputPort(plant, std::move(name), model_value);
 }
 
 template <typename T>
-systems::InputPort<T>& ForceDensityField<T>::DeclareVectorInputPort(
+systems::InputPort<T>& ForceDensityFieldImpl<T>::DeclareVectorInputPort(
     internal::MultibodyTreeSystem<T>* plant, std::string name,
     const systems::BasicVector<T>& model_vector) {
   return internal::MultibodyTreeSystemElementAttorney<
@@ -39,7 +39,7 @@ GravityForceField<T>::~GravityForceField() = default;
 }  // namespace drake
 
 DRAKE_DEFINE_CLASS_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_SCALARS(
-    class ::drake::multibody::ForceDensityField);
+    class ::drake::multibody::ForceDensityFieldImpl);
 
 DRAKE_DEFINE_CLASS_TEMPLATE_INSTANTIATIONS_ON_DEFAULT_SCALARS(
     class ::drake::multibody::GravityForceField);