41 std::vector<double>
const& local_x, std::vector<double>& nodal_values,
42 IntegrationMethod
const& _integration_method, IPData
const& _ip_data,
43 StressData const& stress_data, OutputData
const& output_data,
46 unsigned const n_integration_points =
47 _integration_method.getNumberOfPoints();
51 nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
53 for (
unsigned ip = 0; ip < n_integration_points; ip++)
55 auto const& sigma = stress_data[ip].stress_data.
sigma;
56 auto const& N = _ip_data[ip].N_u;
57 auto const& dNdx = _ip_data[ip].dNdx_u;
60 output_data[ip].free_energy_density_data.free_energy_density;
64 element, _ip_data[ip].N_u);
69 DisplacementDim * DisplacementDim + (DisplacementDim == 2 ? 1 : 0),
70 ShapeFunction::NPOINTS * DisplacementDim);
72 dNdx, G, is_axially_symmetric, N, x_coord);
74 GradientVectorType
const grad_u =
75 G * Eigen::Map<NodalForceVectorType const>(
76 local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
78 GradientVectorType eshelby_stress;
79 eshelby_stress.setZero(DisplacementDim * DisplacementDim +
80 (DisplacementDim == 2 ? 1 : 0));
82 if (DisplacementDim == 3)
84 eshelby_stress[0] = eshelby_stress[DisplacementDim + 1] =
85 eshelby_stress[8] = psi;
87 eshelby_stress[0] -= sigma[0] * grad_u[0] +
88 sigma[3] / std::sqrt(2) * grad_u[3] +
89 sigma[5] / std::sqrt(2) * grad_u[6];
91 eshelby_stress[1] -= sigma[3] / std::sqrt(2) * grad_u[0] +
92 sigma[1] * grad_u[3] +
93 sigma[4] / std::sqrt(2) * grad_u[6];
95 eshelby_stress[2] -= sigma[5] / std::sqrt(2) * grad_u[0] +
96 sigma[4] / std::sqrt(2) * grad_u[3] +
99 eshelby_stress[3] -= sigma[0] * grad_u[1] +
100 sigma[3] / std::sqrt(2) * grad_u[4] +
101 sigma[5] / std::sqrt(2) * grad_u[7];
103 eshelby_stress[4] -= sigma[3] / std::sqrt(2) * grad_u[1] +
104 sigma[1] * grad_u[4] +
105 sigma[4] / std::sqrt(2) * grad_u[7];
107 eshelby_stress[5] -= sigma[5] / std::sqrt(2) * grad_u[1] +
108 sigma[4] / std::sqrt(2) * grad_u[4] +
109 sigma[2] * grad_u[7];
111 eshelby_stress[6] -= sigma[0] * grad_u[2] +
112 sigma[3] / std::sqrt(2) * grad_u[5] +
113 sigma[5] / std::sqrt(2) * grad_u[8];
115 eshelby_stress[7] -= sigma[3] / std::sqrt(2) * grad_u[2] +
116 sigma[1] * grad_u[5] +
117 sigma[4] / std::sqrt(2) * grad_u[8];
119 eshelby_stress[8] -= sigma[5] / std::sqrt(2) * grad_u[2] +
120 sigma[4] / std::sqrt(2) * grad_u[5] +
121 sigma[2] * grad_u[8];
123 else if (DisplacementDim == 2)
125 eshelby_stress[0] = eshelby_stress[DisplacementDim + 1] =
126 eshelby_stress[4] = psi;
129 sigma[0] * grad_u[0] + sigma[3] / std::sqrt(2) * grad_u[2];
132 sigma[3] / std::sqrt(2) * grad_u[0] + sigma[1] * grad_u[2];
135 sigma[0] * grad_u[1] + sigma[3] / std::sqrt(2) * grad_u[3];
138 sigma[3] / std::sqrt(2) * grad_u[1] + sigma[1] * grad_u[3];
141 eshelby_stress[4] -= sigma[2] * grad_u[4];
146 "Material forces not implemented for displacement "
148 "other than 2 and 3.");
151 auto const& w = _ip_data[ip].integration_weight;
152 local_b += G.transpose() * eshelby_stress * w;
160 std::unique_ptr<GlobalVector>& material_forces,
161 std::vector<std::unique_ptr<LocalAssemblerInterface>>
const&
166 DBUG(
"Compute material forces for small deformation process.");
171 if (!material_forces)
180 [](
const std::size_t mesh_item_id,
187 std::vector<double> local_data;
188 auto const local_x = x.
get(indices);
190 local_assembler.getMaterialForces(local_x, local_data);
192 assert(local_data.size() == indices.size());
193 node_values.
add(indices, local_data);
195 local_assemblers, local_to_global_index_map, x, *material_forces);