25 namespace SmallDeformation
30 std::vector<double>
const& local_x,
31 std::vector<double>& nodal_values) = 0;
36 template <
int DisplacementDim,
typename ShapeFunction,
37 typename ShapeMatricesType,
typename NodalForceVectorType,
38 typename NodalDisplacementVectorType,
typename GradientVectorType,
39 typename GradientMatrixType,
typename IPData,
40 typename IntegrationMethod>
42 std::vector<double>
const& local_x, std::vector<double>& nodal_values,
43 IntegrationMethod
const& _integration_method, IPData
const& _ip_data,
46 unsigned const n_integration_points =
47 _integration_method.getNumberOfPoints();
50 auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
51 nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
53 for (
unsigned ip = 0; ip < n_integration_points; ip++)
55 auto const& sigma = _ip_data[ip].sigma;
56 auto const& N = _ip_data[ip].N;
57 auto const& dNdx = _ip_data[ip].dNdx;
59 auto const& psi = _ip_data[ip].free_energy_density;
62 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
63 element, _ip_data[ip].N);
68 DisplacementDim * DisplacementDim + (DisplacementDim == 2 ? 1 : 0),
69 ShapeFunction::NPOINTS * DisplacementDim);
70 Deformation::computeGMatrix<DisplacementDim, ShapeFunction::NPOINTS>(
71 dNdx, G, is_axially_symmetric, N, x_coord);
73 GradientVectorType
const grad_u =
74 G * Eigen::Map<NodalForceVectorType const>(
75 local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
77 GradientVectorType eshelby_stress;
78 eshelby_stress.setZero(DisplacementDim * DisplacementDim +
79 (DisplacementDim == 2 ? 1 : 0));
81 if (DisplacementDim == 3)
83 eshelby_stress[0] = eshelby_stress[DisplacementDim + 1] =
84 eshelby_stress[8] = psi;
86 eshelby_stress[0] -= sigma[0] * grad_u[0] +
87 sigma[3] / std::sqrt(2) * grad_u[3] +
88 sigma[5] / std::sqrt(2) * grad_u[6];
90 eshelby_stress[1] -= sigma[3] / std::sqrt(2) * grad_u[0] +
91 sigma[1] * grad_u[3] +
92 sigma[4] / std::sqrt(2) * grad_u[6];
94 eshelby_stress[2] -= sigma[5] / std::sqrt(2) * grad_u[0] +
95 sigma[4] / std::sqrt(2) * grad_u[3] +
98 eshelby_stress[3] -= sigma[0] * grad_u[1] +
99 sigma[3] / std::sqrt(2) * grad_u[4] +
100 sigma[5] / std::sqrt(2) * grad_u[7];
102 eshelby_stress[4] -= sigma[3] / std::sqrt(2) * grad_u[1] +
103 sigma[1] * grad_u[4] +
104 sigma[4] / std::sqrt(2) * grad_u[7];
106 eshelby_stress[5] -= sigma[5] / std::sqrt(2) * grad_u[1] +
107 sigma[4] / std::sqrt(2) * grad_u[4] +
108 sigma[2] * grad_u[7];
110 eshelby_stress[6] -= sigma[0] * grad_u[2] +
111 sigma[3] / std::sqrt(2) * grad_u[5] +
112 sigma[5] / std::sqrt(2) * grad_u[8];
114 eshelby_stress[7] -= sigma[3] / std::sqrt(2) * grad_u[2] +
115 sigma[1] * grad_u[5] +
116 sigma[4] / std::sqrt(2) * grad_u[8];
118 eshelby_stress[8] -= sigma[5] / std::sqrt(2) * grad_u[2] +
119 sigma[4] / std::sqrt(2) * grad_u[5] +
120 sigma[2] * grad_u[8];
122 else if (DisplacementDim == 2)
124 eshelby_stress[0] = eshelby_stress[DisplacementDim + 1] =
125 eshelby_stress[4] = psi;
128 sigma[0] * grad_u[0] + sigma[3] / std::sqrt(2) * grad_u[2];
131 sigma[3] / std::sqrt(2) * grad_u[0] + sigma[1] * grad_u[2];
134 sigma[0] * grad_u[1] + sigma[3] / std::sqrt(2) * grad_u[3];
137 sigma[3] / std::sqrt(2) * grad_u[1] + sigma[1] * grad_u[3];
140 eshelby_stress[4] -= sigma[2] * grad_u[4];
145 "Material forces not implemented for displacement "
147 "other than 2 and 3.");
150 auto const& w = _ip_data[ip].integration_weight;
151 local_b += G.transpose() * eshelby_stress * w;
157 template <
typename LocalAssemblerInterface>
159 std::unique_ptr<GlobalVector>& material_forces,
160 std::vector<std::unique_ptr<LocalAssemblerInterface>>
const&
165 DBUG(
"Compute material forces for small deformation process.");
170 if (!material_forces)
179 [](
const std::size_t mesh_item_id,
185 std::vector<double> local_data;
186 auto const local_x = x.
get(indices);
188 local_assembler.getMaterialForces(local_x, local_data);
190 assert(local_data.size() == indices.size());
191 node_values.
add(indices, local_data);
193 local_assemblers, local_to_global_index_map, x, *material_forces);
void DBUG(char const *fmt, Args const &... args)
Global vector based on Eigen vector.
void add(IndexType rowId, double v)
add entry
double get(IndexType rowId) const
get entry
void finalizeAssembly(PETScMatrix &A)
void setLocalAccessibleVector(PETScVector const &x)
void set(PETScVector &x, double const a)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
static void executeDereferenced(F const &f, C const &c, Args_ &&... args)