[PL] SD; Add material forces computation.

ProcessLib/Deformation/MaterialForces.h

+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <cassert>
+#include <vector>
+
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "ProcessLib/Deformation/GMatrix.h"
+#include "ProcessLib/Parameter/SpatialPosition.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformation
+{
+struct MaterialForcesInterface
+{
+    virtual std::vector<double> const& getMaterialForces(
+        std::vector<double> const& local_x,
+        std::vector<double>& nodal_values) = 0;
+
+    virtual ~MaterialForcesInterface() = default;
+};
+
+template <int DisplacementDim, typename ShapeFunction,
+          typename ShapeMatricesType, typename NodalForceVectorType,
+          typename NodalDisplacementVectorType, typename GradientVectorType,
+          typename GradientMatrixType, typename IPData,
+          typename IntegrationMethod>
+std::vector<double> const& getMaterialForces(
+    std::vector<double> const& local_x, std::vector<double>& nodal_values,
+    IntegrationMethod const& _integration_method, IPData const& _ip_data,
+    MeshLib::Element const& element, bool const is_axially_symmetric)
+{
+    unsigned const n_integration_points =
+        _integration_method.getNumberOfPoints();
+
+    nodal_values.clear();
+    auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
+        nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
+
+    for (unsigned ip = 0; ip < n_integration_points; ip++)
+    {
+        auto const& sigma = _ip_data[ip].sigma;
+        auto const& N = _ip_data[ip].N;
+        auto const& dNdx = _ip_data[ip].dNdx;
+
+        auto const& psi = _ip_data[ip].free_energy_density;
+
+        auto const x_coord =
+            interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+                element, _ip_data[ip].N);
+
+        // For the 2D case the 33-component is needed (and the four entries
+        // of the non-symmetric matrix); In 3d there are nine entries.
+        GradientMatrixType G(
+            DisplacementDim * DisplacementDim + (DisplacementDim == 2 ? 1 : 0),
+            ShapeFunction::NPOINTS * DisplacementDim);
+        Deformation::computeGMatrix<DisplacementDim, ShapeFunction::NPOINTS>(
+            dNdx, G, is_axially_symmetric, N, x_coord);
+
+        GradientVectorType const grad_u =
+            G * Eigen::Map<NodalForceVectorType const>(
+                    local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
+
+        GradientVectorType eshelby_stress;
+        eshelby_stress.setZero(DisplacementDim * DisplacementDim +
+                               (DisplacementDim == 2 ? 1 : 0));
+
+        if (DisplacementDim == 3)
+        {
+            eshelby_stress[0] = eshelby_stress[DisplacementDim + 1] =
+                eshelby_stress[8] = psi;
+
+            eshelby_stress[0] -= sigma[0] * grad_u[0] +
+                                 sigma[3] / std::sqrt(2) * grad_u[3] +
+                                 sigma[5] / std::sqrt(2) * grad_u[6];
+
+            eshelby_stress[1] -= sigma[3] / std::sqrt(2) * grad_u[0] +
+                                 sigma[1] * grad_u[3] +
+                                 sigma[4] / std::sqrt(2) * grad_u[6];
+
+            eshelby_stress[2] -= sigma[5] / std::sqrt(2) * grad_u[0] +
+                                 sigma[4] / std::sqrt(2) * grad_u[3] +
+                                 sigma[2] * grad_u[6];
+
+            eshelby_stress[3] -= sigma[0] * grad_u[1] +
+                                 sigma[3] / std::sqrt(2) * grad_u[4] +
+                                 sigma[5] / std::sqrt(2) * grad_u[7];
+
+            eshelby_stress[4] -= sigma[3] / std::sqrt(2) * grad_u[1] +
+                                 sigma[1] * grad_u[4] +
+                                 sigma[4] / std::sqrt(2) * grad_u[7];
+
+            eshelby_stress[5] -= sigma[5] / std::sqrt(2) * grad_u[1] +
+                                 sigma[4] / std::sqrt(2) * grad_u[4] +
+                                 sigma[2] * grad_u[7];
+
+            eshelby_stress[6] -= sigma[0] * grad_u[2] +
+                                 sigma[3] / std::sqrt(2) * grad_u[5] +
+                                 sigma[5] / std::sqrt(2) * grad_u[8];
+
+            eshelby_stress[7] -= sigma[3] / std::sqrt(2) * grad_u[2] +
+                                 sigma[1] * grad_u[5] +
+                                 sigma[4] / std::sqrt(2) * grad_u[8];
+
+            eshelby_stress[8] -= sigma[5] / std::sqrt(2) * grad_u[2] +
+                                 sigma[4] / std::sqrt(2) * grad_u[5] +
+                                 sigma[2] * grad_u[8];
+        }
+        else if (DisplacementDim == 2)
+        {
+            eshelby_stress[0] = eshelby_stress[DisplacementDim + 1] =
+                eshelby_stress[4] = psi;
+
+            eshelby_stress[0] -=
+                sigma[0] * grad_u[0] + sigma[3] / std::sqrt(2) * grad_u[2];
+
+            eshelby_stress[1] -=
+                sigma[3] / std::sqrt(2) * grad_u[0] + sigma[1] * grad_u[2];
+
+            eshelby_stress[2] -=
+                sigma[0] * grad_u[1] + sigma[3] / std::sqrt(2) * grad_u[3];
+
+            eshelby_stress[3] -=
+                sigma[3] / std::sqrt(2) * grad_u[1] + sigma[1] * grad_u[3];
+
+            // for axial-symmetric case the following is not zero in general
+            eshelby_stress[4] -= sigma[2] * grad_u[4];
+        }
+        else
+            OGS_FATAL(
+                "Material forces not implemented for displacement dimension "
+                "other than 2 and 3.");
+
+        auto const& w = _ip_data[ip].integration_weight;
+        local_b += G.transpose() * eshelby_stress * w;
+    }
+
+    return nodal_values;
+}
+
+template <typename LocalAssemblerInterface>
+void writeMaterialForces(
+    MeshLib::PropertyVector<double>& material_forces,
+    std::vector<std::unique_ptr<LocalAssemblerInterface>> const&
+        local_assemblers,
+    NumLib::LocalToGlobalIndexMap const& local_to_global_index_map,
+    GlobalVector const& x)
+{
+    DBUG("Compute material forces for small deformation process.");
+
+    // Material forces
+    std::fill(std::begin(material_forces), std::end(material_forces), 0);
+
+    GlobalExecutor::executeDereferenced(
+        [](const std::size_t mesh_item_id,
+           LocalAssemblerInterface& local_assembler,
+           const NumLib::LocalToGlobalIndexMap& dof_table,
+           GlobalVector const& x,
+           std::vector<double>& node_values) {
+            auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
+            std::vector<double> local_data;
+            auto const local_x = x.get(indices);
+
+            local_assembler.getMaterialForces(local_x, local_data);
+
+            assert(local_data.size() == indices.size());
+            for (std::size_t i = 0; i < indices.size(); ++i)
+                node_values[indices[i]] += local_data[i];
+        },
+        local_assemblers, local_to_global_index_map, x, material_forces);
+}
+
+}  // namespace SmallDeformation
+}  // namespace ProcessLib

ProcessLib/SmallDeformation/LocalAssemblerInterface.h

@@ -13,6 +13,7 @@
 
 #include "MaterialLib/SolidModels/MechanicsBase.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "ProcessLib/Deformation/MaterialForces.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
 
 namespace ProcessLib
@@ -22,6 +23,7 @@ namespace SmallDeformation
 template <int DisplacementDim>
 struct SmallDeformationLocalAssemblerInterface
     : public ProcessLib::LocalAssemblerInterface,
+      public MaterialForcesInterface,
       public NumLib::ExtrapolatableElement
 {
     virtual std::vector<double> const& getIntPtFreeEnergyDensity(

ProcessLib/SmallDeformation/SmallDeformationFEM.h

@@ -18,6 +18,7 @@
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Deformation/GMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
@@ -93,6 +94,10 @@ public:
     using NodalDisplacementVectorType =
         typename BMatricesType::NodalForceVectorType;
 
+    using GMatricesType = GMatrixPolicyType<ShapeFunction, DisplacementDim>;
+    using GradientVectorType = typename GMatricesType::GradientVectorType;
+    using GradientMatrixType = typename GMatricesType::GradientMatrixType;
+
     SmallDeformationLocalAssembler(SmallDeformationLocalAssembler const&) =
         delete;
     SmallDeformationLocalAssembler(SmallDeformationLocalAssembler&&) = delete;
@@ -249,6 +254,35 @@ public:
         }
     }
 
+    void postTimestepConcrete(std::vector<double> const& /*local_x*/) override
+    {
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            auto& ip_data = _ip_data[ip];
+
+            // Update free energy density needed for material forces.
+            ip_data.free_energy_density =
+                ip_data.solid_material.computeFreeEnergyDensity(
+                    ip_data.eps, ip_data.sigma,
+                    *ip_data.material_state_variables);
+        }
+    }
+
+    std::vector<double> const& getMaterialForces(
+        std::vector<double> const& local_x,
+        std::vector<double>& nodal_values) override
+    {
+        return ProcessLib::SmallDeformation::getMaterialForces<
+            DisplacementDim, ShapeFunction, ShapeMatricesType,
+            typename BMatricesType::NodalForceVectorType,
+            NodalDisplacementVectorType, GradientVectorType,
+            GradientMatrixType>(local_x, nodal_values, _integration_method,
+                                _ip_data, _element, _is_axially_symmetric);
+    }
+
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
         const unsigned integration_point) const override
     {
@@ -259,6 +293,9 @@ public:
     }
 
     std::vector<double> const& getIntPtFreeEnergyDensity(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
         std::vector<double>& cache) const override
     {
         cache.clear();

ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h

@@ -38,6 +38,8 @@ SmallDeformationProcess<DisplacementDim>::SmallDeformationProcess(
 {
     _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
         mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
+    _material_forces = MeshLib::getOrCreateMeshProperty<double>(
+        mesh, "MaterialForces", MeshLib::MeshItemType::Node, DisplacementDim);
 }
 
 template <int DisplacementDim>
@@ -69,6 +71,13 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
             NumLib::ComponentOrder::BY_LOCATION);
     _nodal_forces->resize(DisplacementDim * mesh.getNumberOfNodes());
 
+    _material_forces->resize(DisplacementDim * mesh.getNumberOfNodes());
+
+    Base::_secondary_variables.addSecondaryVariable(
+        "free_energy_density",
+        makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+                         &LocalAssemblerInterface::getIntPtFreeEnergyDensity));
+
     Base::_secondary_variables.addSecondaryVariable(
         "sigma",
         makeExtrapolator(
@@ -248,5 +257,19 @@ void SmallDeformationProcess<DisplacementDim>::preTimestepConcreteProcess(
         *_local_to_global_index_map, x, t, dt);
 }
 
+template <int DisplacementDim>
+void SmallDeformationProcess<DisplacementDim>::postTimestepConcreteProcess(
+    GlobalVector const& x)
+{
+    DBUG("PostTimestep SmallDeformationProcess.");
+
+    GlobalExecutor::executeMemberOnDereferenced(
+        &LocalAssemblerInterface::postTimestep, _local_assemblers,
+        *_local_to_global_index_map, x);
+
+    ProcessLib::SmallDeformation::writeMaterialForces(
+        *_material_forces, _local_assemblers, *_local_to_global_index_map, x);
+}
+
 }  // namespace SmallDeformation
 }  // namespace ProcessLib

ProcessLib/SmallDeformation/SmallDeformationProcess.h

@@ -61,6 +61,8 @@ private:
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
                                     double const dt) override;
 
+    void postTimestepConcreteProcess(GlobalVector const& x) override;
+
 private:
     SmallDeformationProcessData<DisplacementDim> _process_data;
 
@@ -69,6 +71,7 @@ private:
     std::unique_ptr<NumLib::LocalToGlobalIndexMap>
         _local_to_global_index_map_single_component;
     MeshLib::PropertyVector<double>* _nodal_forces = nullptr;
+    MeshLib::PropertyVector<double>* _material_forces = nullptr;
 };
 
 extern template class SmallDeformationProcess<2>;