From 8efb45a191ea9bae079f235522400f3b26751c0e Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Wed, 5 Apr 2017 13:13:35 +0200
Subject: [PATCH] [PL] SD; Add material forces computation.
---
ProcessLib/Deformation/MaterialForces.h | 187 ++++++++++++++++++
.../LocalAssemblerInterface.h | 2 +
.../SmallDeformation/SmallDeformationFEM.h | 37 ++++
.../SmallDeformationProcess-impl.h | 23 +++
.../SmallDeformationProcess.h | 3 +
5 files changed, 252 insertions(+)
create mode 100644 ProcessLib/Deformation/MaterialForces.h
diff --git a/ProcessLib/Deformation/MaterialForces.h b/ProcessLib/Deformation/MaterialForces.h
new file mode 100644
index 00000000000..7545e2a18ef
--- /dev/null
+++ b/ProcessLib/Deformation/MaterialForces.h
@@ -0,0 +1,187 @@
+/**
+ * \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
diff --git a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
index cfc6e7b4f71..fd8347879d7 100644
--- a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
+++ b/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(
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index 48c31de3f9b..e454243f6b7 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/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();
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h b/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h
index b962917fbfb..2c495666b10 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h
+++ b/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
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
index a1bb0e86a61..534ac993ee0 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.h
+++ b/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>;
--
GitLab