[PL] PF: Add energy computation at postTimestep.

MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h

@@ -94,7 +94,8 @@ public:
                                  KelvinMatrix& C_tensile,
                                  KelvinMatrix& C_compressive,
                                  KelvinVector& sigma_real,
-                                 double const degradation) const override
+                                 double const degradation,
+                                 double& elastic_energy) const override
     {
         using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
         // calculation of deviatoric parts
@@ -120,6 +121,7 @@ public:
             sigma_compressive.noalias() = KelvinVector::Zero();
             C_tensile.template topLeftCorner<3, 3>().setConstant(K);
             C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
+            elastic_energy = degradation * strain_energy_tensile;
         }
         else
         {
@@ -129,6 +131,8 @@ public:
                 K * eps_curr_trace * Invariants::identity2;
             C_tensile.noalias() = 2 * mu * P_dev * KelvinMatrix::Identity();
             C_compressive.template topLeftCorner<3, 3>().setConstant(K);
+            elastic_energy = K / 2 * eps_curr_trace * eps_curr_trace +
+                             mu * epsd_curr.transpose() * epsd_curr;
         }
 
         sigma_real.noalias() = degradation * sigma_tensile + sigma_compressive;

MaterialLib/SolidModels/PhaseFieldExtension.h

@@ -32,7 +32,8 @@ struct PhaseFieldExtension : public MechanicsBase<DisplacementDim>
                                          KelvinMatrix& C_tensile,
                                          KelvinMatrix& C_compressive,
                                          KelvinVector& sigma_real,
-                                         double const degradation) const = 0;
+                                         double const degradation,
+                                         double& elastic_energy) const = 0;
 
     /// Dynamic size Kelvin vector and matrix wrapper for the polymorphic
     /// constitutive relation compute function.
@@ -48,7 +49,8 @@ struct PhaseFieldExtension : public MechanicsBase<DisplacementDim>
         Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
             C_compressive,
         Eigen::Matrix<double, Eigen::Dynamic, 1>& sigma_real,
-        double const degradation) const
+        double const degradation,
+        double& elastic_energy) const
     {
         // TODO Avoid copies of data:
         // Using MatrixBase<Derived> not possible because template functions
@@ -71,7 +73,8 @@ struct PhaseFieldExtension : public MechanicsBase<DisplacementDim>
                                                     C_tensile_,
                                                     C_compressive_,
                                                     sigma_real_,
-                                                    degradation);
+                                                    degradation,
+                                                    elastic_energy);
 
         sigma_tensile = sigma_tensile_;
         sigma_compressive = sigma_compressive_;

ProcessLib/PhaseField/LocalAssemblerInterface.h

@@ -102,6 +102,16 @@ struct PhaseFieldLocalAssemblerInterface
         GlobalVector const& x, double const t, double& crack_volume,
         bool const use_monolithic_scheme,
         CoupledSolutionsForStaggeredScheme const* const cpl_xs) = 0;
+
+    virtual void computeEnergy(
+        std::size_t mesh_item_id,
+        std::vector<
+            std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
+            dof_tables,
+        GlobalVector const& x, double const t, double& elastic_energy,
+        double& surface_energy, double& pressure_work,
+        bool const use_monolithic_scheme,
+        CoupledSolutionsForStaggeredScheme const* const cpl_xs) = 0;
 };
 
 }  // namespace PhaseField

ProcessLib/PhaseField/PhaseFieldFEM-impl.h

@@ -502,5 +502,84 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
         crack_volume += (N_u * u).dot(dNdx * d) * w;
     }
 }
+
+template <typename ShapeFunction, typename IntegrationMethod,
+          int DisplacementDim>
+void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
+                              DisplacementDim>::
+    computeEnergy(std::size_t mesh_item_id,
+                  std::vector<std::reference_wrapper<
+                      NumLib::LocalToGlobalIndexMap>> const& dof_tables,
+                  GlobalVector const& /*x*/, double const t,
+                  double& elastic_energy, double& surface_energy,
+                  double& pressure_work, bool const /*use_monolithic_scheme*/,
+                  CoupledSolutionsForStaggeredScheme const* const cpl_xs)
+{
+    assert(cpl_xs != nullptr);
+
+    std::vector<std::vector<GlobalIndexType>> indices_of_processes;
+    indices_of_processes.reserve(dof_tables.size());
+    std::transform(dof_tables.begin(), dof_tables.end(),
+                   std::back_inserter(indices_of_processes),
+                   [&](NumLib::LocalToGlobalIndexMap const& dof_table) {
+                       return NumLib::getIndices(mesh_item_id, dof_table);
+                   });
+
+    auto local_coupled_xs =
+        getCurrentLocalSolutions(*cpl_xs, indices_of_processes);
+    assert(local_coupled_xs.size() == 2);
+
+    auto const& local_u = local_coupled_xs[0];
+    auto const& local_d = local_coupled_xs[1];
+
+    assert(local_u.size() == displacement_size);
+    assert(local_d.size() == phasefield_size);
+
+    auto d = Eigen::Map<
+        typename ShapeMatricesType::template VectorType<phasefield_size> const>(
+        local_d.data(), phasefield_size);
+
+    auto u = Eigen::Map<typename ShapeMatricesType::template VectorType<
+        displacement_size> const>(local_u.data(), displacement_size);
+
+    SpatialPosition x_position;
+    x_position.setElementID(_element.getID());
+
+    int const n_integration_points = _integration_method.getNumberOfPoints();
+    for (int ip = 0; ip < n_integration_points; ip++)
+    {
+        x_position.setIntegrationPoint(ip);
+        auto const& w = _ip_data[ip].integration_weight;
+        auto const& N = _ip_data[ip].N;
+        auto const& dNdx = _ip_data[ip].dNdx;
+        double const d_ip = N.dot(d);
+        auto pressure_ip = _process_data.pressure;
+        auto u_corrected = pressure_ip * u;
+        double const gc = _process_data.crack_resistance(t, x_position)[0];
+        double const ls = _process_data.crack_length_scale(t, x_position)[0];
+
+        typename ShapeMatricesType::template MatrixType<DisplacementDim,
+                                                        displacement_size>
+            N_u = ShapeMatricesType::template MatrixType<
+                DisplacementDim, displacement_size>::Zero(DisplacementDim,
+                                                          displacement_size);
+
+        for (int i = 0; i < DisplacementDim; ++i)
+            N_u.template block<1, displacement_size / DisplacementDim>(
+                   i, i * displacement_size / DisplacementDim)
+                .noalias() = N;
+
+        elastic_energy += _ip_data[ip].elastic_energy * w;
+
+        surface_energy +=
+            0.5 * gc *
+            ((1 - d_ip) * (1 - d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) *
+            w;
+
+        if (_process_data.crack_pressure)
+            pressure_work +=
+                pressure_ip * (N_u * u_corrected).dot(dNdx * d) * w;
+    }
+}
 }  // namespace PhaseField
 }  // namespace ProcessLib

ProcessLib/PhaseField/PhaseFieldFEM.h

@@ -46,7 +46,7 @@ struct IntegrationPointData final
 
     typename BMatricesType::KelvinVectorType sigma_tensile, sigma_compressive,
         sigma_real_prev, sigma_real;
-    double strain_energy_tensile;
+    double strain_energy_tensile, elastic_energy;
 
     MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material;
     std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
@@ -81,7 +81,7 @@ struct IntegrationPointData final
             .calculateDegradedStress(t, x_position, eps, strain_energy_tensile,
                                      sigma_tensile, sigma_compressive,
                                      C_tensile, C_compressive, sigma_real,
-                                     degradation);
+                                     degradation, elastic_energy);
     }
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
 };
@@ -161,6 +161,7 @@ public:
                 _process_data.history_field(0, x_position)[0];
             ip_data.sigma_real.setZero(kelvin_vector_size);
             ip_data.strain_energy_tensile = 0.0;
+            ip_data.elastic_energy = 0.0;
 
             ip_data.N = shape_matrices[ip].N;
             ip_data.dNdx = shape_matrices[ip].dNdx;
@@ -217,6 +218,16 @@ public:
         bool const use_monolithic_scheme,
         CoupledSolutionsForStaggeredScheme const* const cpl_xs) override;
 
+    void computeEnergy(
+        std::size_t mesh_item_id,
+        std::vector<
+            std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
+            dof_tables,
+        GlobalVector const& x, double const t, double& elastic_energy,
+        double& surface_energy, double& pressure_work,
+        bool const use_monolithic_scheme,
+        CoupledSolutionsForStaggeredScheme const* const cpl_xs) override;
+
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
         const unsigned integration_point) const override
     {

ProcessLib/PhaseField/PhaseFieldProcess-impl.h

@@ -360,11 +360,30 @@ template <int DisplacementDim>
 void PhaseFieldProcess<DisplacementDim>::postTimestepConcreteProcess(
     GlobalVector const& x, int const process_id)
 {
-    DBUG("PostTimestep PhaseFieldProcess %d.", process_id);
+    if (isPhaseFieldProcess(process_id))
+    {
+        DBUG("PostTimestep PhaseFieldProcess.");
 
-    GlobalExecutor::executeMemberOnDereferenced(
-        &LocalAssemblerInterface::postTimestep, _local_assemblers,
-        getDOFTable(process_id), x);
+        _process_data.elastic_energy = 0.0;
+        _process_data.surface_energy = 0.0;
+        _process_data.pressure_work = 0.0;
+
+        std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+            dof_tables;
+
+        dof_tables.emplace_back(*_local_to_global_index_map);
+        dof_tables.emplace_back(*_local_to_global_index_map_single_component);
+
+        GlobalExecutor::executeMemberOnDereferenced(
+            &LocalAssemblerInterface::computeEnergy, _local_assemblers,
+            dof_tables, x, _process_data.t, _process_data.elastic_energy,
+            _process_data.surface_energy, _process_data.pressure_work,
+            _use_monolithic_scheme, _coupled_solutions);
+
+        INFO("Elastic energy: %g Surface energy: %g Pressure work: %g ",
+             _process_data.elastic_energy, _process_data.surface_energy,
+             _process_data.pressure_work);
+    }
 }
 
 template <int DisplacementDim>

ProcessLib/PhaseField/PhaseFieldProcessData.h

@@ -98,6 +98,9 @@ struct PhaseFieldProcessData
     double pressure_error = 0.0;
     double injected_volume = 0.0;
     double crack_volume = 0.0;
+    double elastic_energy = 0.0;
+    double surface_energy = 0.0;
+    double pressure_work = 0.0;
     bool propagating_crack = false;
     bool crack_pressure = false;
 };