131 {
133
134
135 auto get_ents_on_mesh_skin = [&]() {
136 std::map<std::string, Range> boundary_entities_vec;
137
138 auto get_block_vec_impl = [&](auto block_name) {
140 std::regex(
141
142 (boost::format("%s(.*)") % block_name).str()
143
144 ));
145 };
146
147 auto add_blockset_ents_impl = [&](auto &&vec) {
149 for (auto it : vec) {
150 Range boundary_entities;
152 boundary_entities, true);
153 std::string meshset_name = it->getName();
154 boundary_entities_vec[meshset_name] = boundary_entities;
155 boundary_entities.clear();
156 }
158 };
159
160 for (auto block_name :
161 {"SPATIAL_DISP", "FIX", "CONTACT", "SPATIAL_ROTATION",
162 "NORMAL_DISPLACEMENT", "ANALYTICAL_DISPLACEMENT", "SPRING_BC"}) {
164 "add blockset entities");
165 }
166
167 return boundary_entities_vec;
168 };
169
170 auto boundary_entities_vec = get_ents_on_mesh_skin();
171
172 std::vector<std::tuple<std::string, Range, std::array<double, 6>>>
173 reactionForces;
174 for (const auto &pair : boundary_entities_vec) {
175 reactionForces.push_back(
176 std::make_tuple(pair.first, pair.second,
177 std::array<double, 6>{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}));
178 }
179
182 };
183 auto face_fe = boost::make_shared<FaceElementForcesAndSourcesCore>(
eP.
mField);
184 face_fe->getRuleHook = integration_rule_face;
186 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
189
190 auto op_side =
192 face_fe->getOpPtrVector().push_back(op_side);
193 auto side_fe_ptr = op_side->getSideFEPtr();
194 auto base_ptr =
195 boost::make_shared<EshelbianPlasticity::CGGUserPolynomialBase>();
196 side_fe_ptr->getUserPolynomialBase() = base_ptr;
198 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
201 auto piola_scale_ptr = boost::make_shared<double>(1.0);
202 side_fe_ptr->getOpPtrVector().push_back(
203 new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
205 constexpr bool add_bubble_field = false;
206 if (add_bubble_field) {
207 side_fe_ptr->getOpPtrVector().push_back(
208 new OpCalculateHTensorTensorField<3, 3>(
210 }
211 side_fe_ptr->getOpPtrVector().push_back(
213 side_fe_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
215 for (
auto &[
name, ents, reaction_vec] : reactionForces) {
216 face_fe->getOpPtrVector().push_back(
218 }
219
221 *face_fe);
222
223 for (
auto &[
name, ents, reaction_vec] : reactionForces) {
224 std::array<double, 6> block_reaction_force{0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
225
226 MPI_Allreduce(reaction_vec.data(), &block_reaction_force, 6, MPI_DOUBLE,
228
229 for (auto &force : block_reaction_force) {
230 if (std::abs(force) < 1e-12) {
231 force = 0.0;
232 }
233 }
234
237 "Step %d time %3.4g Block %s Reaction force [%3.6e, "
238 "%3.6e, %3.6e]",
239 ts_step, ts_t,
name.c_str(), block_reaction_force[0],
240 block_reaction_force[1], block_reaction_force[2]);
242 "Step %d time %3.4g Block %s Moment [%3.6e, %3.6e, %3.6e]",
243 ts_step, ts_t,
name.c_str(), block_reaction_force[3],
244 block_reaction_force[4], block_reaction_force[5]);
245 }
247};
static constexpr int approx_order
boost::shared_ptr< Range > frontAdjEdges
const std::string spatialL2Disp
const std::string materialH1Positions
const std::string piolaStress
const std::string bubbleField
const std::string skinElement
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
std::map< std::string, std::array< double, 6 > > reactionForcesMap
