Newer
Older
/**
* \file
*
* \copyright
* Copyright (c) 2012-2023, 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 "MathLib/LinAlg/LinAlg.h"
#include "MathLib/LinAlg/MatrixVectorTraits.h"
#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Fem/InitShapeMatrices.h"
#include "ProcessLib/Deformation/GMatrix.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 StressData,
typename OutputData, 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,
StressData const& stress_data, OutputData const& output_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 = stress_data[ip].stress_data.sigma;
auto const& N = _ip_data[ip].N;
auto const& dNdx = _ip_data[ip].dNdx;
auto const& psi =
output_data[ip].free_energy_density_data.free_energy_density;
NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
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
"Material forces not implemented for displacement "
"dimension "
auto const& w = _ip_data[ip].integration_weight;
local_b += G.transpose() * eshelby_stress * w;
}
return nodal_values;
}
template <typename LocalAssemblerInterface>
void writeMaterialForces(
std::unique_ptr<GlobalVector>& 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.");
// Prepare local storage and fetch values.
MathLib::LinAlg::setLocalAccessibleVector(x);
if (!material_forces)
{
material_forces =
MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
}
MathLib::LinAlg::set(*material_forces, 0);
GlobalExecutor::executeDereferenced(
[](const std::size_t mesh_item_id,
LocalAssemblerInterface& local_assembler,
const NumLib::LocalToGlobalIndexMap& dof_table,
GlobalVector const& x,
GlobalVector& 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());
node_values.add(indices, local_data);
local_assemblers, local_to_global_index_map, x, *material_forces);
MathLib::LinAlg::finalizeAssembly(*material_forces);
}
} // namespace SmallDeformation
} // namespace ProcessLib