14#include <boost/python.hpp>
15#include <boost/python/def.hpp>
16#include <boost/python/numpy.hpp>
17namespace bp = boost::python;
18namespace np = boost::python::numpy;
34 IntegrationType::GAUSS;
85template <
int SPACE_DIM, IntegrationType I,
typename OpBase>
89template <
int SPACE_DIM, IntegrationType I,
typename OpBase>
129 boost::shared_ptr<MatrixDouble> u_ptr,
130 boost::shared_ptr<MatrixDouble> stress_ptr,
131 boost::shared_ptr<MatrixDouble> strain_ptr,
132 boost::shared_ptr<VectorDouble> o_ptr) = 0;
137 boost::shared_ptr<MatrixDouble> u_ptr,
138 boost::shared_ptr<MatrixDouble> stress_ptr,
139 boost::shared_ptr<MatrixDouble> strain_ptr,
140 boost::shared_ptr<MatrixDouble> o_ptr) = 0;
145 boost::shared_ptr<MatrixDouble> u_ptr,
146 boost::shared_ptr<MatrixDouble> stress_ptr,
147 boost::shared_ptr<MatrixDouble> strain_ptr,
148 boost::shared_ptr<MatrixDouble> o_ptr) = 0;
155 std::array<double, 3> ¢roid,
156 std::array<double, 6> &bbox,
162boost::shared_ptr<ObjectiveFunctionData>
166 const std::string block_name,
int dim);
210 Vec objective_function_gradient,
253 char objective_function_file_name[255] =
"objective_function.py";
255 PETSC_NULLPTR, PETSC_NULLPTR,
"-objective_function",
256 objective_function_file_name, 255, PETSC_NULLPTR);
259 auto file_exists = [](std::string myfile) {
260 std::ifstream file(myfile.c_str());
266 if (!file_exists(objective_function_file_name)) {
267 MOFEM_LOG(
"WORLD", Sev::error) <<
"Objective function file NOT found: "
268 << objective_function_file_name;
291 auto create_vec_modes = [&](
auto block_name) {
296 std::regex((boost::format(
"%s(.*)") % block_name).str())
300 int nb_total_modes = 0;
301 for (
auto &bc : bcs) {
302 auto id = bc->getMeshsetId();
305 nb_total_modes += nb_modes;
309 <<
"Total number of modes to apply: " << nb_total_modes;
316 auto get_modes_bounding_boxes = [&](
auto block_name) {
321 std::regex((boost::format(
"%s(.*)") % block_name).str())
325 for (
auto &bc : bcs) {
326 auto meshset = bc->getMeshset();
331 std::vector<double> x(verts.size());
332 std::vector<double> y(verts.size());
333 std::vector<double> z(verts.size());
335 std::array<double, 3> centroid = {0, 0, 0};
336 for (
int i = 0;
i != verts.size(); ++
i) {
341 MPI_Allreduce(MPI_IN_PLACE, centroid.data(), 3, MPI_DOUBLE, MPI_SUM,
343 int nb_vertex = verts.size();
344 MPI_Allreduce(MPI_IN_PLACE, &nb_vertex, 1, MPI_INT, MPI_SUM,
347 centroid[0] /= nb_vertex;
348 centroid[1] /= nb_vertex;
349 centroid[2] /= nb_vertex;
351 std::array<double, 6> bbox = {centroid[0], centroid[1], centroid[2],
352 centroid[0], centroid[1], centroid[2]};
353 for (
int i = 0;
i != verts.size(); ++
i) {
354 bbox[0] = std::min(bbox[0], x[
i]);
355 bbox[1] = std::min(bbox[1], y[
i]);
356 bbox[2] = std::min(bbox[2], z[
i]);
357 bbox[3] = std::max(bbox[3], x[
i]);
358 bbox[4] = std::max(bbox[4], y[
i]);
359 bbox[5] = std::max(bbox[5], z[
i]);
361 MPI_Allreduce(MPI_IN_PLACE, &bbox[0], 3, MPI_DOUBLE, MPI_MIN,
363 MPI_Allreduce(MPI_IN_PLACE, &bbox[3], 3, MPI_DOUBLE, MPI_MAX,
367 <<
"Block: " << bc->getName() <<
" centroid: " << centroid[0] <<
" "
368 << centroid[1] <<
" " << centroid[2];
370 <<
"Block: " << bc->getName() <<
" bbox: " << bbox[0] <<
" "
371 << bbox[1] <<
" " << bbox[2] <<
" " << bbox[3] <<
" " << bbox[4]
381 auto eval_objective_and_gradient = [](Tao tao,
Vec x, PetscReal *
f,
Vec g,
382 void *ctx) -> PetscErrorCode {
387 CHKERR TaoGetIterationNumber(tao, &iter);
389 auto set_geometry = [&](
Vec x) {
394 CHKERR VecScatterCreateToAll(x, &ctx, &x_local);
396 CHKERR VecScatterBegin(ctx, x, x_local, INSERT_VALUES, SCATTER_FORWARD);
397 CHKERR VecScatterEnd(ctx, x, x_local, INSERT_VALUES, SCATTER_FORWARD);
399 CHKERR VecScatterDestroy(&ctx);
402 CHKERR VecCopy(ex_ptr->initialGeometry, current_geometry);
404 CHKERR VecGetArrayRead(x_local, &
a);
405 const double *coeff =
a;
406 for (
auto &mode_vec : ex_ptr->
modeVecs) {
408 <<
"Adding mode with coeff: " << *coeff;
409 CHKERR VecAXPY(current_geometry, (*coeff), mode_vec);
412 CHKERR VecRestoreArrayRead(x_local, &
a);
413 CHKERR VecGhostUpdateBegin(current_geometry, INSERT_VALUES,
415 CHKERR VecGhostUpdateEnd(current_geometry, INSERT_VALUES,
418 ->setOtherLocalGhostVector(
"ADJOINT",
"ADJOINT_FIELD",
"GEOMETRY",
420 INSERT_VALUES, SCATTER_REVERSE);
422 CHKERR VecDestroy(&x_local);
427 CHKERR KSPReset(ex_ptr->kspElastic);
428 CHKERR ex_ptr->solveElastic();
431 CHKERR ex_ptr->calculateGradient(f,
g, adjoint_vector);
432 CHKERR ex_ptr->postprocessElastic(iter, adjoint_vector);
438 CHKERR create_vec_modes(
"OPTIMISE");
439 CHKERR get_modes_bounding_boxes(
"OPTIMISE");
452 CHKERR TaoSetType(tao, TAOLMVM);
464 INSERT_VALUES, SCATTER_FORWARD);
493 CHKERR TaoSetSolution(tao, x0);
494 CHKERR TaoSetFromOptions(tao);
498 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Topology optimization completed";
551 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
552 const char *list_bases[LASBASETOPT] = {
"ainsworth",
"demkowicz"};
553 PetscInt choice_base_value = AINSWORTH;
555 LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
557 switch (choice_base_value) {
561 <<
"Set AINSWORTH_LEGENDRE_BASE for displacements";
566 <<
"Set DEMKOWICZ_JACOBI_BASE for displacements";
600 auto project_ho_geometry = [&]() {
604 CHKERR project_ho_geometry();
635 auto create_adjoint_dm = [&]() {
638 auto add_field = [&]() {
643 for (
auto tt = MBEDGE; tt <= moab::CN::TypeDimensionMap[
SPACE_DIM].second;
653 auto add_adjoint_fe_impl = [&]() {
678 auto block_name =
"OPTIMISE";
682 std::regex((boost::format(
"%s(.*)") % block_name).str())
686 for (
auto bc : bcs) {
688 bc->getMeshset(),
SPACE_DIM - 1,
"ADJOINT_BOUNDARY_FE");
694 simple->getBitRefLevelMask());
699 auto set_adjoint_dm_imp = [&]() {
703 simple->getBitRefLevelMask());
705 CHKERR DMSetFromOptions(adjoint_dm);
712 CHKERR DMSetUp(adjoint_dm);
737 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_X",
739 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Y",
741 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Z",
743 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
744 "REMOVE_ALL",
"U", 0, 3);
746 simple->getProblemName(),
"U");
762 boost::make_shared<Range>(opt_ents));
764 boost::make_shared<Range>(opt_ents));
765 CHKERR DMSetUp(subset_dm_bdy);
773 CHKERR DMSetUp(subset_dm_domain);
776 auto remove_dofs = [&]() {
779 std::array<Range, 3> remove_dim_ents;
787 for (
int d = 0; d != 3; ++d) {
789 <<
"Removing topology modes on block OPT_REMOVE_" << (char)(
'X' + d)
790 <<
" with " << remove_dim_ents[d].size() <<
" entities";
798 CHKERR skin.find_skin(0, body_ents,
false, boundary_ents);
799 for (
int d = 0; d != 3; ++d) {
800 boundary_ents = subtract(boundary_ents, remove_dim_ents[d]);
802 ParallelComm *pcomm =
804 CHKERR pcomm->filter_pstatus(boundary_ents,
805 PSTATUS_SHARED | PSTATUS_MULTISHARED,
806 PSTATUS_NOT, -1,
nullptr);
807 for (
auto d =
SPACE_DIM - 2; d >= 0; --d) {
811 moab::Interface::UNION);
812 boundary_ents.merge(ents);
816 boundary_ents.merge(verts);
821 boundary_ents.merge(opt_ents);
823 "SUBSET_DOMAIN",
"ADJOINT_FIELD", boundary_ents);
824 for (
int d = 0; d != 3; ++d) {
826 "SUBSET_DOMAIN",
"ADJOINT_FIELD", remove_dim_ents[d], d, d);
841 auto get_lhs_fe = [&]() {
842 auto fe_lhs = boost::make_shared<BoundaryEle>(
mField);
843 fe_lhs->getRuleHook = [](int, int,
int p_data) {
844 return 2 * p_data + p_data - 1;
846 auto &pip = fe_lhs->getOpPtrVector();
851 pip.push_back(
new OpMass(
"ADJOINT_FIELD",
"ADJOINT_FIELD",
852 [](
double,
double,
double) {
return 1.; }));
856 auto get_rhs_fe = [&]() {
857 auto fe_rhs = boost::make_shared<BoundaryEle>(
mField);
858 fe_rhs->getRuleHook = [](int, int,
int p_data) {
859 return 2 * p_data + p_data - 1;
861 auto &pip = fe_rhs->getOpPtrVector();
868 auto block_name =
"OPTIMISE";
872 std::regex((boost::format(
"%s(.*)") % block_name).str())
881 struct OpMode :
public OP {
882 OpMode(
const std::string name,
883 boost::shared_ptr<ObjectiveFunctionData> python_ptr,
int id,
885 std::vector<std::array<double, 3>> mode_centroids,
886 std::vector<std::array<double, 6>> mode_bboxes,
int block_counter,
887 int mode_counter, boost::shared_ptr<Range> range =
nullptr)
888 :
OP(name, name, OP::OPROW, range), pythonPtr(python_ptr), iD(
id),
889 modeVecs(mode_vecs), modeCentroids(mode_centroids),
890 modeBboxes(mode_bboxes), blockCounter(block_counter),
891 modeCounter(mode_counter) {}
897 if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
905 auto nb_base_functions = data.
getN().size2();
908 CHKERR pythonPtr->blockModes(iD, OP::getCoordsAtGaussPts(),
909 modeCentroids[blockCounter],
910 modeBboxes[blockCounter], blockModes);
912 auto nb_integration_pts = getGaussPts().size2();
913 if (blockModes.size2() != 3 * nb_integration_pts) {
915 <<
"Number of modes does not match number of integration points: "
916 << blockModes.size2() <<
"!=" << 3 * nb_integration_pts;
922 int nb_modes = blockModes.size1();
923 for (
auto mode = 0; mode != nb_modes; ++mode) {
926 auto t_mode = getFTensor1FromPtr<3>(&blockModes(mode, 0));
928 const double vol = OP::getMeasure();
930 auto t_w = OP::getFTensor0IntegrationWeight();
934 for (
int gg = 0; gg != nb_integration_pts; gg++) {
937 const double alpha = t_w * vol;
939 auto t_nf = getFTensor1FromPtr<SPACE_DIM>(nf.data().data());
941 for (; rr != nb_rows /
SPACE_DIM; ++rr) {
942 t_nf(
i) += alpha * t_base * t_mode(
i);
946 for (; rr < nb_base_functions; ++rr)
951 Vec vec = modeVecs[modeCounter + mode];
953 auto *indices = data.
getIndices().data().data();
954 auto *nf_data = nf.data().data();
962 boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
964 std::vector<std::array<double, 3>> modeCentroids;
965 std::vector<std::array<double, 6>> modeBboxes;
967 std::vector<SmartPetscObj<Vec>> modeVecs;
972 auto solve_bdy = [&]() {
975 auto fe_lhs = get_lhs_fe();
976 auto fe_rhs = get_rhs_fe();
977 int block_counter = 0;
978 int mode_counter = 0;
979 for (
auto &bc : bcs) {
980 auto id = bc->getMeshsetId();
984 auto range = boost::make_shared<Range>(ents);
985 auto &pip_rhs = fe_rhs->getOpPtrVector();
988 mode_counter, range));
995 mode_counter += nb_modes;
997 <<
"Setting mode block block: " << bc->getName()
998 <<
" with ID: " << bc->getMeshsetId()
999 <<
" total modes: " << mode_counter;
1005 CHKERR VecGhostUpdateBegin(
v, ADD_VALUES, SCATTER_REVERSE);
1006 CHKERR VecGhostUpdateEnd(
v, ADD_VALUES, SCATTER_REVERSE);
1013 CHKERR MatAssemblyBegin(
M, MAT_FINAL_ASSEMBLY);
1014 CHKERR MatAssemblyEnd(
M, MAT_FINAL_ASSEMBLY);
1017 CHKERR KSPSetOperators(solver,
M,
M);
1018 CHKERR KSPSetFromOptions(solver);
1022 CHKERR KSPSolve(solver, f,
v);
1027 CHKERR VecGhostUpdateBegin(
v, INSERT_VALUES, SCATTER_FORWARD);
1028 CHKERR VecGhostUpdateEnd(
v, INSERT_VALUES, SCATTER_FORWARD);
1036 auto get_elastic_fe_lhs = [&]() {
1037 auto fe = boost::make_shared<DomainEle>(
mField);
1038 fe->getRuleHook = [](int, int,
int p_data) {
1039 return 2 * p_data + p_data - 1;
1041 auto &pip = fe->getOpPtrVector();
1044 CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, A, I, DomainEleOp>(
1045 mField, pip,
"ADJOINT_FIELD",
"MAT_ADJOINT", Sev::noisy);
1049 auto get_elastic_fe_rhs = [&]() {
1050 auto fe = boost::make_shared<DomainEle>(
mField);
1051 fe->getRuleHook = [](int, int,
int p_data) {
1052 return 2 * p_data + p_data - 1;
1054 auto &pip = fe->getOpPtrVector();
1057 CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
1058 mField, pip,
"ADJOINT_FIELD",
"MAT_ADJOINT", Sev::noisy);
1062 auto adjoint_gradient_postprocess = [&](
auto mode) {
1064 auto post_proc_mesh = boost::make_shared<moab::Core>();
1065 auto post_proc_begin =
1066 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
mField,
1068 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
1071 auto geom_vec = boost::make_shared<MatrixDouble>();
1074 boost::make_shared<PostProcEleDomain>(
mField, post_proc_mesh);
1076 post_proc_fe->getOpPtrVector(), {H1},
"GEOMETRY");
1077 post_proc_fe->getOpPtrVector().push_back(
1083 post_proc_fe->getOpPtrVector().push_back(
1087 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1091 {{
"MODE", geom_vec}},
1102 post_proc_begin->getFEMethod());
1106 post_proc_begin->getFEMethod());
1108 CHKERR post_proc_end->writeFile(
"mode_" + std::to_string(mode) +
".h5m");
1113 auto solve_domain = [&]() {
1115 auto fe_lhs = get_elastic_fe_lhs();
1116 auto fe_rhs = get_elastic_fe_rhs();
1125 CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
1126 CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
1128 auto solver =
createKSP(mField.get_comm());
1129 CHKERR KSPSetOperators(solver, M, M);
1130 CHKERR KSPSetFromOptions(solver);
1133 int mode_counter = 0;
1134 for (
auto &f : modeVecs) {
1135 CHKERR mField.getInterface<
FieldBlas>()->setField(0,
"ADJOINT_FIELD");
1143 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
1144 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
1146 CHKERR VecGhostUpdateBegin(
v, INSERT_VALUES, SCATTER_FORWARD);
1147 CHKERR VecGhostUpdateEnd(
v, INSERT_VALUES, SCATTER_FORWARD);
1161 for (
int i = 0;
i < modeVecs.size(); ++
i) {
1162 CHKERR adjoint_gradient_postprocess(
i);
1185 CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
1186 CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
1190 pip->getOpDomainLhsPipeline(), {H1},
"GEOMETRY");
1192 pip->getOpDomainRhsPipeline(), {H1},
"GEOMETRY");
1194 pip->getOpBoundaryRhsPipeline(), {NOSPACE},
"GEOMETRY");
1196 pip->getOpBoundaryLhsPipeline(), {NOSPACE},
"GEOMETRY");
1200 CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, A, I, DomainEleOp>(
1201 mField, pip->getOpDomainLhsPipeline(),
"U",
"MAT_ELASTIC", Sev::noisy);
1205 CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
1206 mField, pip->getOpDomainRhsPipeline(),
"U",
"MAT_ELASTIC", Sev::noisy);
1210 pip->getOpDomainRhsPipeline(),
mField,
"U", Sev::inform);
1216 pip->getOpBoundaryRhsPipeline(),
mField,
"U", -1, Sev::inform);
1219 pip->getOpBoundaryLhsPipeline(),
mField,
"U", Sev::noisy);
1232 CHKERR VecZeroEntries(d);
1235 auto set_essential_bc = [&]() {
1240 auto pre_proc_rhs = boost::make_shared<FEMethod>();
1241 auto post_proc_rhs = boost::make_shared<FEMethod>();
1242 auto post_proc_lhs = boost::make_shared<FEMethod>();
1244 auto get_pre_proc_hook = [&]() {
1248 pre_proc_rhs->preProcessHook = get_pre_proc_hook();
1250 auto get_post_proc_hook_rhs = [
this, post_proc_rhs]() {
1254 post_proc_rhs, 1.)();
1258 auto get_post_proc_hook_lhs = [
this, post_proc_lhs]() {
1262 post_proc_lhs, 1.)();
1266 post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
1267 post_proc_lhs->postProcessHook = get_post_proc_hook_lhs;
1269 ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
1270 ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
1271 ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
1275 auto setup_and_solve = [&](
auto solver) {
1277 BOOST_LOG_SCOPED_THREAD_ATTR(
"Timeline", attrs::timer());
1278 MOFEM_LOG(
"TIMER", Sev::noisy) <<
"KSPSetUp";
1280 MOFEM_LOG(
"TIMER", Sev::noisy) <<
"KSPSetUp <= Done";
1281 MOFEM_LOG(
"TIMER", Sev::noisy) <<
"KSPSolve";
1282 CHKERR KSPSolve(solver, f, d);
1283 MOFEM_LOG(
"TIMER", Sev::noisy) <<
"KSPSolve <= Done";
1290 CHKERR set_essential_bc();
1293 CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
1294 CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
1297 auto evaluate_field_at_the_point = [&]() {
1301 std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
1304 field_eval_coords.data(), &coords_dim,
1310 auto field_eval_data =
1314 ->buildTree<SPACE_DIM>(field_eval_data,
simple->getDomainFEName());
1316 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1317 auto no_rule = [](int, int, int) {
return -1; };
1318 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
1319 field_eval_fe_ptr->getRuleHook = no_rule;
1321 field_eval_fe_ptr->getOpPtrVector().push_back(
1325 ->evalFEAtThePoint<SPACE_DIM>(
1326 field_eval_coords.data(), 1e-12,
simple->getProblemName(),
1327 simple->getDomainFEName(), field_eval_data,
1334 MOFEM_LOG(
"FieldEvaluator", Sev::inform)
1335 <<
"U_X: " << t_disp(0) <<
" U_Y: " << t_disp(1);
1337 MOFEM_LOG(
"FieldEvaluator", Sev::inform)
1338 <<
"U_X: " << t_disp(0) <<
" U_Y: " << t_disp(1)
1339 <<
" U_Z: " << t_disp(2);
1347 CHKERR evaluate_field_at_the_point();
1359 auto det_ptr = boost::make_shared<VectorDouble>();
1360 auto jac_ptr = boost::make_shared<MatrixDouble>();
1361 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1363 pip->getDomainRhsFE().reset();
1364 pip->getDomainLhsFE().reset();
1365 pip->getBoundaryRhsFE().reset();
1366 pip->getBoundaryLhsFE().reset();
1370 auto post_proc_mesh = boost::make_shared<moab::Core>();
1371 auto post_proc_begin = boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
1373 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
1378 auto calculate_stress_ops = [&](
auto &pip) {
1383 auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEle>(
1384 mField, pip,
"U",
"MAT_ELASTIC", Sev::noisy);
1387 auto u_ptr = boost::make_shared<MatrixDouble>();
1389 auto x_ptr = boost::make_shared<MatrixDouble>();
1393 using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
1395 DataMapMat vec_map = {{
"U", u_ptr}, {
"GEOMETRY", x_ptr}};
1396 DataMapMat sym_map = {{
"STRAIN", common_ptr->getMatStrain()},
1397 {
"STRESS", common_ptr->getMatCauchyStress()}};
1399 if (adjoint_vector) {
1400 auto adjoint_ptr = boost::make_shared<MatrixDouble>();
1402 "U", adjoint_ptr, adjoint_vector));
1403 vec_map[
"ADJOINT"] = adjoint_ptr;
1407 return boost::make_tuple(u_ptr, x_ptr, common_ptr->getMatStrain(),
1408 common_ptr->getMatCauchyStress(), vec_map,
1412 auto get_tag_id_on_pmesh = [&](
bool post_proc_skin) {
1413 int def_val_int = 0;
1416 "MAT_ELASTIC", 1, MB_TYPE_INTEGER, tag_mat,
1417 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
1418 auto meshset_vec_ptr =
1420 std::regex((boost::format(
"%s(.*)") %
"MAT_ELASTIC").str()));
1423 std::unique_ptr<Skinner> skin_ptr;
1424 if (post_proc_skin) {
1426 auto boundary_meshset =
simple->getBoundaryMeshSet();
1431 for (
auto m : meshset_vec_ptr) {
1435 int const id =
m->getMeshsetId();
1436 ents_3d = ents_3d.subset_by_dimension(
SPACE_DIM);
1437 if (post_proc_skin) {
1439 CHKERR skin_ptr->find_skin(0, ents_3d,
false, skin_faces);
1440 ents_3d = intersect(skin_ents, skin_faces);
1448 auto post_proc_domain = [&](
auto post_proc_mesh) {
1450 boost::make_shared<PostProcEleDomain>(
mField, post_proc_mesh);
1453 auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr, vec_map, sym_map] =
1454 calculate_stress_ops(post_proc_fe->getOpPtrVector());
1456 post_proc_fe->getOpPtrVector().push_back(
1460 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1474 post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(
false)});
1475 return post_proc_fe;
1478 auto post_proc_boundary = [&](
auto post_proc_mesh) {
1480 boost::make_shared<PostProcEleBdy>(
mField, post_proc_mesh);
1482 post_proc_fe->getOpPtrVector(), {},
"GEOMETRY");
1486 auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr, vec_map, sym_map] =
1487 calculate_stress_ops(op_loop_side->getOpPtrVector());
1488 post_proc_fe->getOpPtrVector().push_back(op_loop_side);
1489 auto mat_traction_ptr = boost::make_shared<MatrixDouble>();
1490 post_proc_fe->getOpPtrVector().push_back(
1492 vec_map[
"T"] = mat_traction_ptr;
1496 post_proc_fe->getOpPtrVector().push_back(
1500 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1514 post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(
true)});
1515 return post_proc_fe;
1518 PetscBool post_proc_skin_only = PETSC_FALSE;
1520 post_proc_skin_only = PETSC_TRUE;
1522 &post_proc_skin_only, PETSC_NULLPTR);
1524 if (post_proc_skin_only == PETSC_FALSE) {
1525 pip->getDomainRhsFE() = post_proc_domain(post_proc_mesh);
1527 pip->getBoundaryRhsFE() = post_proc_boundary(post_proc_mesh);
1531 post_proc_begin->getFEMethod());
1532 CHKERR pip->loopFiniteElements();
1534 post_proc_end->getFEMethod());
1536 CHKERR post_proc_end->writeFile(
"out_elastic_" + std::to_string(iter) +
1559template <
int SPACE_DIM>
1566 boost::shared_ptr<HookeOps::CommonData> comm_ptr)
1585template <
int SPACE_DIM>
1599 const double vol = OP::getMeasure();
1601 auto t_w = OP::getFTensor0IntegrationWeight();
1605 auto t_cauchy_stress =
1606 getFTensor2SymmetricFromMat<SPACE_DIM>(*(commPtr->getMatCauchyStress()));
1609 for (
int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1611 const double alpha = t_w * vol;
1614 auto t_nf = OP::template getNf<SPACE_DIM>();
1617 for (; rr != OP::nbRows /
SPACE_DIM; rr++) {
1619 t_nf(
j) += alpha * t_row_grad(
i) * t_cauchy_stress(
i,
j);
1624 for (; rr < OP::nbRowBaseFunctions; ++rr) {
1634template <
int SPACE_DIM>
1641 boost::shared_ptr<HookeOps::CommonData> comm_ptr,
1642 boost::shared_ptr<MatrixDouble> jac,
1643 boost::shared_ptr<MatrixDouble> diff_jac,
1644 boost::shared_ptr<VectorDouble> cof_vals)
1646 jac(jac), diffJac(diff_jac), cofVals(cof_vals) {}
1650 boost::shared_ptr<MatrixDouble>
jac;
1665 t_diff(
i,
j,
k,
l) = 0;
1666 t_diff(0, 0, 0, 0) = 1;
1667 t_diff(1, 1, 1, 1) = 1;
1669 t_diff(1, 0, 1, 0) = 0.5;
1670 t_diff(1, 0, 0, 1) = 0.5;
1672 t_diff(0, 1, 0, 1) = 0.5;
1673 t_diff(0, 1, 1, 0) = 0.5;
1675 if constexpr (DIM == 3) {
1676 t_diff(2, 2, 2, 2) = 1;
1678 t_diff(2, 0, 2, 0) = 0.5;
1679 t_diff(2, 0, 0, 2) = 0.5;
1680 t_diff(0, 2, 0, 2) = 0.5;
1681 t_diff(0, 2, 2, 0) = 0.5;
1683 t_diff(2, 1, 2, 1) = 0.5;
1684 t_diff(2, 1, 1, 2) = 0.5;
1685 t_diff(1, 2, 1, 2) = 0.5;
1686 t_diff(1, 2, 2, 1) = 0.5;
1692template <
int SPACE_DIM>
1707 const double vol = OP::getMeasure();
1709 auto t_w = OP::getFTensor0IntegrationWeight();
1711 auto t_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(jac));
1713 auto t_diff_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(diffJac));
1720 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(commPtr->matGradPtr));
1722 auto t_cauchy_stress =
1723 getFTensor2SymmetricFromMat<SPACE_DIM>(*(commPtr->getMatCauchyStress()));
1726 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*(commPtr->matDPtr));
1728 for (
int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1730 const double alpha = t_w * vol;
1738 t_diff_inv_jac(
i,
j) =
1739 -(t_inv_jac(
i,
I) * t_diff_jac(
I,
J)) * t_inv_jac(
J,
j);
1741 t_diff_grad(
i,
j) = t_grad_u(
i,
k) * t_diff_inv_jac(
k,
j);
1745 t_diff_strain(
i,
j) = t_diff_symm(
i,
j,
k,
l) * t_diff_grad(
k,
l);
1749 t_diff_stress(
i,
j) = t_D(
i,
j,
k,
l) * t_diff_strain(
k,
l);
1752 auto t_nf = OP::template getNf<SPACE_DIM>();
1755 for (; rr != OP::nbRows /
SPACE_DIM; rr++) {
1758 t_diff_row_grad(
k) = t_row_grad(
j) * t_diff_inv_jac(
j,
k);
1761 t_nf(
j) += alpha * t_diff_row_grad(
i) * t_cauchy_stress(
i,
j);
1764 t_nf(
j) += (alpha * t_cof) * t_row_grad(
i) * t_cauchy_stress(
i,
j);
1767 t_nf(
j) += alpha * t_row_grad(
i) * t_diff_stress(
i,
j);
1772 for (; rr < OP::nbRowBaseFunctions; ++rr) {
1805 boost::shared_ptr<HookeOps::CommonData> comm_ptr,
1806 boost::shared_ptr<MatrixDouble> jac_ptr,
1807 boost::shared_ptr<MatrixDouble> diff_jac,
1808 boost::shared_ptr<VectorDouble> cof_vals,
1809 boost::shared_ptr<MatrixDouble> d_grad_ptr,
1810 boost::shared_ptr<MatrixDouble> u_ptr,
1811 boost::shared_ptr<double> glob_objective_ptr,
1812 boost::shared_ptr<double> glob_objective_grad_ptr)
1841 auto nb_gauss_pts = getGaussPts().size2();
1842 auto objective_ptr = boost::make_shared<VectorDouble>(nb_gauss_pts);
1843 auto objective_dstress =
1844 boost::make_shared<MatrixDouble>(symm_size, nb_gauss_pts);
1845 auto objective_dstrain =
1846 boost::make_shared<MatrixDouble>(symm_size, nb_gauss_pts);
1847 auto obj_grad = boost::make_shared<MatrixDouble>(
SPACE_DIM, nb_gauss_pts);
1849 auto evaluate_python = [&]() {
1851 auto &coords = OP::getCoordsAtGaussPts();
1863 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(
commPtr->matGradPtr));
1865 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*(
commPtr->matDPtr));
1866 auto t_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(
jacPtr));
1867 auto t_diff_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(
diffJacPtr));
1869 auto t_d_grad = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(
dGradPtr));
1872 auto t_obj_dstress =
1873 getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstress);
1874 auto t_obj_dstrain =
1875 getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstrain);
1877 auto vol = OP::getMeasure();
1878 auto t_w = getFTensor0IntegrationWeight();
1879 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
1886 t_diff_inv_jac(
i,
j) =
1887 -(t_inv_jac(
i,
I) * t_diff_jac(
I,
J)) * t_inv_jac(
J,
j);
1889 t_diff_grad(
i,
j) = t_grad_u(
i,
k) * t_diff_inv_jac(
k,
j);
1892 t_d_strain(
i,
j) = t_diff_symm(
i,
j,
k,
l) * (
1902 auto alpha = t_w * vol;
1904 (*globObjectivePtr) += alpha * t_obj;
1905 (*globObjectiveGradPtr) +=
1909 t_obj_dstress(
i,
j) * (t_D(
i,
j,
k,
l) * t_d_strain(
k,
l))
1913 t_obj_dstrain(
i,
j) * t_d_strain(
i,
j)
1936 CHKERR evaluate_python();
1948 boost::shared_ptr<MatrixDouble>
uPtr;
1955 Vec objective_function_gradient,
1956 Vec adjoint_vector) {
1963 auto get_essential_fe = [
this]() {
1964 auto post_proc_rhs = boost::make_shared<FEMethod>();
1965 auto get_post_proc_hook_rhs = [
this, post_proc_rhs]() {
1969 post_proc_rhs, 0)();
1972 post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
1973 return post_proc_rhs;
1976 auto get_fd_direvative_fe = [&]() {
1977 auto fe = boost::make_shared<DomainEle>(
mField);
1978 fe->getRuleHook = [](int, int,
int p_data) {
1979 return 2 * p_data + p_data - 1;
1981 auto &pip = fe->getOpPtrVector();
1984 constexpr bool debug =
false;
1985 if constexpr (
debug) {
1986 auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
1987 mField, pip,
"U",
"MAT_ELASTIC", Sev::noisy);
1991 HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
1992 mField, pip,
"U",
"MAT_ELASTIC", Sev::noisy);
1997 auto calulate_fd_residual = [&](
auto eps,
auto diff_vec,
auto fd_vec) {
2000 constexpr bool debug =
false;
2006 auto norm2_field = [&](
const double val) {
2011 MPI_Allreduce(MPI_IN_PLACE, &nrm2, 1, MPI_DOUBLE, MPI_SUM,
2013 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Geometry norm: " << sqrt(nrm2);
2017 if constexpr (
debug)
2023 initial_current_geometry, INSERT_VALUES, SCATTER_FORWARD);
2024 CHKERR VecAssemblyBegin(initial_current_geometry);
2025 CHKERR VecAssemblyEnd(initial_current_geometry);
2027 if constexpr (
debug)
2030 auto perturb_geometry = [&](
auto eps,
auto diff_vec) {
2033 CHKERR VecCopy(initial_current_geometry, current_geometry);
2034 CHKERR VecAXPY(current_geometry,
eps, diff_vec);
2037 current_geometry, INSERT_VALUES, SCATTER_REVERSE);
2041 auto fe = get_fd_direvative_fe();
2044 auto calc_impl = [&](
auto f,
auto eps) {
2046 CHKERR VecZeroEntries(f);
2050 simple->getDomainFEName(), fe);
2051 CHKERR VecAssemblyBegin(f);
2052 CHKERR VecAssemblyEnd(f);
2053 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
2054 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
2055 auto post_proc_rhs = get_essential_fe();
2056 post_proc_rhs->f = f;
2058 post_proc_rhs.get());
2063 CHKERR VecWAXPY(fd_vec, -1.0, fm, fp);
2064 CHKERR VecScale(fd_vec, 1.0 / (2.0 *
eps));
2068 initial_current_geometry, INSERT_VALUES, SCATTER_REVERSE);
2070 if constexpr (
debug)
2076 auto get_direvative_fe = [&](
auto diff_vec) {
2077 auto fe_adjoint = boost::make_shared<DomainEle>(
mField);
2078 fe_adjoint->getRuleHook = [](int, int,
int p_data) {
2079 return 2 * p_data + p_data - 1;
2081 auto &pip = fe_adjoint->getOpPtrVector();
2083 auto jac_ptr = boost::make_shared<MatrixDouble>();
2084 auto det_ptr = boost::make_shared<VectorDouble>();
2085 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
2086 auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
2087 auto cof_ptr = boost::make_shared<VectorDouble>();
2095 "GEOMETRY", jac_ptr));
2101 pip.push_back(
new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
2104 auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
2105 mField, pip,
"U",
"MAT_ELASTIC", Sev::noisy);
2107 "U", common_ptr, jac_ptr, diff_jac_ptr, cof_ptr));
2112 auto get_objective_fe = [&](
auto diff_vec,
auto grad_vec,
2113 auto glob_objective_ptr,
2114 auto glob_objective_grad_ptr) {
2115 auto fe_adjoint = boost::make_shared<DomainEle>(
mField);
2116 fe_adjoint->getRuleHook = [](int, int,
int p_data) {
2117 return 2 * p_data + p_data - 1;
2119 auto &pip = fe_adjoint->getOpPtrVector();
2122 auto jac_ptr = boost::make_shared<MatrixDouble>();
2123 auto det_ptr = boost::make_shared<VectorDouble>();
2124 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
2125 auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
2126 auto cof_ptr = boost::make_shared<VectorDouble>();
2127 auto d_grad_ptr = boost::make_shared<MatrixDouble>();
2128 auto u_ptr = boost::make_shared<MatrixDouble>();
2133 "GEOMETRY", jac_ptr));
2138 pip.push_back(
new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
2140 "U", d_grad_ptr, grad_vec));
2143 auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
2144 mField, pip,
"U",
"MAT_ELASTIC", Sev::noisy);
2146 pythonPtr, common_ptr, jac_ptr, diff_jac_ptr, cof_ptr, d_grad_ptr,
2147 u_ptr, glob_objective_ptr, glob_objective_grad_ptr));
2152 auto dm =
simple->getDM();
2157 auto adjoint_fe = get_direvative_fe(dm_diff_vec);
2158 auto glob_objective_ptr = boost::make_shared<double>(0.0);
2159 auto glob_objective_grad_ptr = boost::make_shared<double>(0.0);
2160 auto objective_fe = get_objective_fe(dm_diff_vec, d, glob_objective_ptr,
2161 glob_objective_grad_ptr);
2163 auto set_variance_of_geometry = [&](
auto mode,
auto mod_vec) {
2169 dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
2170 CHKERR VecGhostUpdateBegin(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
2171 CHKERR VecGhostUpdateEnd(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
2175 auto calculate_variance_internal_forces = [&](
auto mode,
auto mod_vec) {
2177 CHKERR VecZeroEntries(f);
2178 CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
2179 CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
2182 CHKERR VecAssemblyBegin(f);
2183 CHKERR VecAssemblyEnd(f);
2184 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
2185 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
2186 auto post_proc_rhs = get_essential_fe();
2187 post_proc_rhs->f = f;
2189 post_proc_rhs.get());
2190 CHKERR VecScale(f, -1.0);
2193 constexpr bool debug =
false;
2194 if constexpr (
debug) {
2196 CHKERR VecNorm(f, NORM_2, &norm0);
2199 CHKERR calulate_fd_residual(
eps, dm_diff_vec, fd_check);
2201 CHKERR VecAXPY(fd_check, -1.0, f);
2202 CHKERR VecNorm(fd_check, NORM_2, &nrm);
2204 <<
" FD check for internal forces [ " << mode <<
" ]: " << nrm
2205 <<
" / " << norm0 <<
" ( " << (nrm / norm0) <<
" )";
2212 auto calulate_variance_of_displacement = [&](
auto mode,
auto mod_vec) {
2215 CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
2216 CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
2220 auto calculate_variance_of_objective_function = [&](
auto mode,
auto mod_vec) {
2222 *glob_objective_ptr = 0.0;
2223 *glob_objective_grad_ptr = 0.0;
2226 CHKERR VecSetValue(objective_function_gradient, mode,
2227 *glob_objective_grad_ptr, ADD_VALUES);
2228 std::array<double, 2> array = {*glob_objective_ptr,
2229 *glob_objective_grad_ptr};
2230 MPI_Allreduce(MPI_IN_PLACE, array.data(), 2, MPI_DOUBLE, MPI_SUM,
2232 *glob_objective_ptr = array[0];
2233 *glob_objective_grad_ptr = array[1];
2234 (*objective_function_value) += *glob_objective_ptr;
2236 MOFEM_LOG(
"WORLD", Sev::verbose) <<
" Objective gradient [ " << mode
2237 <<
" ]: " << *glob_objective_grad_ptr;
2241 CHKERR VecZeroEntries(objective_function_gradient);
2242 CHKERR VecZeroEntries(adjoint_vector);
2243 *objective_function_value = 0.0;
2247 CHKERR set_variance_of_geometry(mode, mod_vec);
2248 CHKERR calculate_variance_internal_forces(mode, mod_vec);
2249 CHKERR calulate_variance_of_displacement(mode, mod_vec);
2250 CHKERR calculate_variance_of_objective_function(mode, mod_vec);
2252 CHKERR VecAXPY(adjoint_vector, *glob_objective_grad_ptr, dm_diff_vec);
2256 (*objective_function_value) /=
modeVecs.size();
2258 <<
"Objective function: " << *glob_objective_ptr;
2260 CHKERR VecAssemblyBegin(objective_function_gradient);
2261 CHKERR VecAssemblyEnd(objective_function_gradient);
2263 CHKERR VecAssemblyBegin(adjoint_vector);
2264 CHKERR VecAssemblyEnd(adjoint_vector);
2265 CHKERR VecGhostUpdateBegin(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
2266 CHKERR VecGhostUpdateEnd(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
2311 const char param_file[] =
"param_file.petsc";
2314 auto core_log = logging::core::get();
2326 DMType dm_name =
"DMMOFEM";
2328 DMType dm_name_mg =
"DMMOFEM_MG";
2333 moab::Core mb_instance;
2334 moab::Interface &moab = mb_instance;
2351 if (Py_FinalizeEx() < 0) {
2407 boost::shared_ptr<MatrixDouble> u_ptr,
2408 boost::shared_ptr<MatrixDouble> stress_ptr,
2409 boost::shared_ptr<MatrixDouble> strain_ptr,
2410 boost::shared_ptr<VectorDouble> o_ptr);
2428 boost::shared_ptr<MatrixDouble> u_ptr,
2429 boost::shared_ptr<MatrixDouble> stress_ptr,
2430 boost::shared_ptr<MatrixDouble> strain_ptr,
2431 boost::shared_ptr<MatrixDouble> o_ptr);
2449 boost::shared_ptr<MatrixDouble> u_ptr,
2450 boost::shared_ptr<MatrixDouble> stress_ptr,
2451 boost::shared_ptr<MatrixDouble> strain_ptr,
2452 boost::shared_ptr<MatrixDouble> o_ptr);
2470 boost::shared_ptr<MatrixDouble> u_ptr,
2471 boost::shared_ptr<MatrixDouble> stress_ptr,
2472 boost::shared_ptr<MatrixDouble> strain_ptr,
2473 boost::shared_ptr<MatrixDouble> o_ptr);
2493 std::array<double, 3> ¢roid,
2522 np::ndarray coords, np::ndarray u,
2524 np::ndarray stress, np::ndarray strain, np::ndarray &o
2543 np::ndarray coords, np::ndarray u,
2545 np::ndarray stress, np::ndarray strain, np::ndarray &o
2564 np::ndarray coords, np::ndarray u,
2566 np::ndarray stress, np::ndarray strain, np::ndarray &o
2586 np::ndarray coords, np::ndarray u,
2588 np::ndarray stress, np::ndarray strain, np::ndarray &o
2608 int block_id, np::ndarray coords, np::ndarray centroid, np::ndarray bbodx,
2666boost::shared_ptr<ObjectiveFunctionData>
2668 auto ptr = boost::make_shared<ObjectiveFunctionDataImpl>();
2702 auto main_module = bp::import(
"__main__");
2716 }
catch (bp::error_already_set
const &) {
2729 auto t_f = getFTensor2FromMat<3, 3>(f);
2730 auto t_s = getFTensor2SymmetricFromMat<SPACE_DIM>(s);
2731 for (
int ii = 0; ii != s.size2(); ++ii) {
2732 t_f(
i,
j) = t_s(
i,
j);
2744 ptr + 0 * s.size2(), ptr + 1 * s.size2(),
2745 ptr + 2 * s.size2(), ptr + 3 * s.size2(),
2746 ptr + 4 * s.size2(), ptr + 5 * s.size2(),
2747 ptr + 6 * s.size2(), ptr + 7 * s.size2(),
2751 auto t_s = getFTensor2SymmetricFromMat<SPACE_DIM>(s);
2752 for (
int ii = 0; ii != s.size2(); ++ii) {
2753 t_s(
i,
j) = (t_f(
i,
j) || t_f(
j,
i)) / 2.0;
2787 MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
2788 boost::shared_ptr<MatrixDouble> stress_ptr,
2789 boost::shared_ptr<MatrixDouble> strain_ptr,
2790 boost::shared_ptr<VectorDouble> o_ptr) {
2798 auto np_u =
convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
2802 auto full_stress =
copyToFull(*(stress_ptr));
2803 auto full_strain =
copyToFull(*(strain_ptr));
2806 auto np_stress =
convertToNumPy(full_stress.data(), full_stress.size1(),
2807 full_stress.size2());
2808 auto np_strain =
convertToNumPy(full_strain.data(), full_strain.size1(),
2809 full_strain.size2());
2812 np::ndarray np_output = np::empty(bp::make_tuple(strain_ptr->size2()),
2813 np::dtype::get_builtin<double>());
2820 o_ptr->resize(stress_ptr->size2(),
false);
2821 double *val_ptr =
reinterpret_cast<double *
>(np_output.get_data());
2822 std::copy(val_ptr, val_ptr + strain_ptr->size2(), o_ptr->data().begin());
2824 }
catch (bp::error_already_set
const &) {
2863 MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
2864 boost::shared_ptr<MatrixDouble> stress_ptr,
2865 boost::shared_ptr<MatrixDouble> strain_ptr,
2866 boost::shared_ptr<MatrixDouble> o_ptr) {
2873 auto np_u =
convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
2876 auto full_stress =
copyToFull(*(stress_ptr));
2877 auto full_strain =
copyToFull(*(strain_ptr));
2880 auto np_stress =
convertToNumPy(full_stress.data(), full_stress.size1(),
2881 full_stress.size2());
2882 auto np_strain =
convertToNumPy(full_strain.data(), full_strain.size1(),
2883 full_strain.size2());
2886 np::ndarray np_output =
2887 np::empty(bp::make_tuple(full_strain.size1(), full_strain.size2()),
2888 np::dtype::get_builtin<double>());
2895 o_ptr->resize(stress_ptr->size1(), stress_ptr->size2(),
false);
2896 double *val_ptr =
reinterpret_cast<double *
>(np_output.get_data());
2900 }
catch (bp::error_already_set
const &) {
2940 MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
2941 boost::shared_ptr<MatrixDouble> stress_ptr,
2942 boost::shared_ptr<MatrixDouble> strain_ptr,
2943 boost::shared_ptr<MatrixDouble> o_ptr) {
2950 auto np_u =
convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
2953 auto full_stress =
copyToFull(*(stress_ptr));
2954 auto full_strain =
copyToFull(*(strain_ptr));
2955 auto np_stress =
convertToNumPy(full_stress.data(), full_stress.size1(),
2956 full_stress.size2());
2957 auto np_strain =
convertToNumPy(full_strain.data(), full_strain.size1(),
2958 full_strain.size2());
2961 np::ndarray np_output =
2962 np::empty(bp::make_tuple(full_strain.size1(), full_strain.size2()),
2963 np::dtype::get_builtin<double>());
2968 o_ptr->resize(strain_ptr->size1(), strain_ptr->size2(),
false);
2969 double *val_ptr =
reinterpret_cast<double *
>(np_output.get_data());
2972 }
catch (bp::error_already_set
const &) {
3014 MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
3015 boost::shared_ptr<MatrixDouble> stress_ptr,
3016 boost::shared_ptr<MatrixDouble> strain_ptr,
3017 boost::shared_ptr<MatrixDouble> o_ptr) {
3024 auto np_u =
convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
3027 auto full_stress =
copyToFull(*(stress_ptr));
3028 auto full_strain =
copyToFull(*(strain_ptr));
3029 auto np_stress =
convertToNumPy(full_stress.data(), full_stress.size1(),
3030 full_stress.size2());
3031 auto np_strain =
convertToNumPy(full_strain.data(), full_strain.size1(),
3032 full_strain.size2());
3035 np::ndarray np_output =
3036 np::empty(bp::make_tuple(u_ptr->size1(), u_ptr->size2()),
3037 np::dtype::get_builtin<double>());
3045 o_ptr->resize(u_ptr->size1(), u_ptr->size2(),
false);
3046 double *val_ptr =
reinterpret_cast<double *
>(np_output.get_data());
3047 std::copy(val_ptr, val_ptr + u_ptr->size1() * u_ptr->size2(),
3048 o_ptr->data().begin());
3050 }
catch (bp::error_already_set
const &) {
3094 int block_id,
MatrixDouble &coords, std::array<double, 3> ¢roid,
3111 np::ndarray np_output =
3112 np::empty(bp::make_tuple(coords.size1(), nb_modes, 3),
3113 np::dtype::get_builtin<double>());
3120 o_ptr.resize(nb_modes, coords.size1() * coords.size2(),
false);
3121 double *val_ptr =
reinterpret_cast<double *
>(np_output.get_data());
3123 std::copy(val_ptr, val_ptr + coords.size1() * coords.size2() * nb_modes,
3124 o_ptr.data().begin());
3126 }
catch (bp::error_already_set
const &) {
3136 np::ndarray coords, np::ndarray u,
3138 np::ndarray stress, np::ndarray strain, np::ndarray &o
3147 }
catch (bp::error_already_set
const &) {
3157 np::ndarray coords, np::ndarray u,
3159 np::ndarray stress, np::ndarray strain, np::ndarray &o
3166 o = bp::extract<np::ndarray>(
3169 }
catch (bp::error_already_set
const &) {
3179 np::ndarray coords, np::ndarray u,
3181 np::ndarray stress, np::ndarray strain, np::ndarray &o
3188 o = bp::extract<np::ndarray>(
3191 }
catch (bp::error_already_set
const &) {
3201 np::ndarray coords, np::ndarray u,
3203 np::ndarray stress, np::ndarray strain, np::ndarray &o
3212 }
catch (bp::error_already_set
const &) {
3227 }
catch (bp::error_already_set
const &) {
3237 np::ndarray centroid,
3243 o = bp::extract<np::ndarray>(
3245 }
catch (bp::error_already_set
const &) {
3276 auto dtype = np::dtype::get_builtin<double>();
3277 auto size = bp::make_tuple(rows, nb_gauss_pts);
3278 auto stride = bp::make_tuple(nb_gauss_pts *
sizeof(
double),
sizeof(
double));
3279 return (np::from_data(data.data(), dtype, size, stride, bp::object()));
3284 auto dtype = np::dtype::get_builtin<double>();
3285 auto size = bp::make_tuple(s);
3286 auto stride = bp::make_tuple(
sizeof(
double));
3287 return (np::from_data(ptr, dtype, size, stride, bp::object()));
3291 const std::string block_name,
int dim) {
3297 std::regex((boost::format(
"%s(.*)") % block_name).str())
3301 for (
auto bc : bcs) {
3305 "get meshset ents");
3309 for (
auto dd = dim - 1; dd >= 0; --dd) {
3313 moab::Interface::UNION),
3333 CHKERR moab.add_entities(*out_meshset, r);
3334 CHKERR moab.write_file(name.c_str(),
"VTK",
"", out_meshset->get_ptr(), 1);
Calculate traction for linear problem.
Implementation of elastic spring bc.
Natural boundary condition applying pressure from fluid.
Implementation of Hookes operator Hookes for linear elastic problems in MoFEM.
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Implementation of natural boundary conditions.
Boundary conditions in domain, i.e. body forces.
#define FTENSOR_INDEX(DIM, I)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
auto diff_symmetrize(FTensor::Number< DIM >)
static char help[]
[calculateGradient]
PetscBool is_plane_strain
constexpr int SPACE_DIM
[Define dimension]
constexpr double poisson_ratio
Poisson's ratio ν
constexpr int BASE_DIM
[Constants and material properties]
constexpr double shear_modulus_G
Shear modulus G = E/(2(1+ν))
constexpr IntegrationType I
Use Gauss quadrature for integration.
constexpr double bulk_modulus_K
Bulk modulus K = E/(3(1-2ν))
boost::shared_ptr< ObjectiveFunctionData > create_python_objective_function(std::string py_file)
Factory function to create Python-integrated objective function interface.
constexpr AssemblyType A
[Define dimension]
Range get_range_from_block(MoFEM::Interface &m_field, const std::string block_name, int dim)
constexpr double young_modulus
[Material properties for linear elasticity]
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
auto createDMMatrix(DM dm)
Get smart matrix from DM.
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
auto createKSP(MPI_Comm comm)
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
PetscErrorCode TaoSetObjectiveAndGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
Sets the objective function value and gradient for a TAO optimization solver.
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
auto createTao(MPI_Comm comm)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
PetscBool do_eval_field
Evaluate field.
ElementsAndOps< SPACE_DIM >::SideEle SideEle
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
FTensor::Index< 'm', 3 > m
Boundary conditions marker.
MoFEMErrorCode boundaryCondition()
Apply essential boundary conditions.
MoFEMErrorCode assembleSystem()
Setup operators in finite element pipeline.
MoFEMErrorCode readMesh()
Read mesh from file and setup meshsets.
SmartPetscObj< KSP > kspElastic
Linear solver for elastic problem.
std::vector< SmartPetscObj< Vec > > modeVecs
Topology mode vectors (design variables)
FieldApproximationBase base
Choice of finite element basis functions
std::vector< std::array< double, 3 > > modeCentroids
Centroids of optimization blocks
MoFEMErrorCode topologyModes()
Compute topology optimization modes.
SmartPetscObj< Vec > initialGeometry
Initial geometry field.
int fieldOrder
Polynomial order for approximation.
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
Main driver function for the optimization process.
MoFEMErrorCode calculateGradient(PetscReal *objective_function_value, Vec objective_function_gradient, Vec adjoint_vector)
Calculate objective function gradient using adjoint method.
MoFEMErrorCode setupAdJoint()
Setup adjoint fields and finite elements
MoFEM::Interface & mField
Reference to MoFEM interface.
std::vector< std::array< double, 6 > > modeBBoxes
Bounding boxes of optimization blocks.
boost::shared_ptr< ObjectiveFunctionData > pythonPtr
Interface to Python objective function.
SmartPetscObj< DM > adjointDM
Data manager for adjoint problem.
MoFEMErrorCode setupProblem()
Setup fields, approximation spaces and DOFs.
MoFEMErrorCode postprocessElastic(int iter, SmartPetscObj< Vec > adjoint_vector=nullptr)
Post-process and output results.
MoFEMErrorCode solveElastic()
Solve forward elastic problem.
boost::shared_ptr< MatrixDouble > vectorFieldPtr
Field values at evaluation points.
Add operators pushing bases from local to physical configuration.
Boundary condition manager for finite element problem setup.
MoFEMErrorCode synchroniseEntities(Range &ent, std::map< int, Range > *received_ents, int verb=DEFAULT_VERBOSITY)
synchronize entity range on processors (collective)
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Definition of the displacement bc data structure.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
Class (Function) to enforce essential constrains on the left hand side diagonal.
Class (Function) to enforce essential constrains on the right hand side diagonal.
Class (Function) to enforce essential constrains.
MoFEMErrorCode fieldLambdaOnValues(OneFieldFunctionOnValues lambda, const std::string field_name, Range *ents_ptr=nullptr)
field lambda
Field evaluator interface.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Interface for managing meshsets containing materials and boundary conditions.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Specialization for double precision vector field values calculation.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
Template struct for dimension-specific finite element types.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Set integration rule for domain right-hand side finite element.
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Implementation of ObjectiveFunctionData interface using Python integration.
MoFEMErrorCode objectiveGradientStrainImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
Internal implementation for strain gradient computation.
MoFEMErrorCode numberOfModes(int block_id, int &modes)
Return number of topology optimization modes for given material block.
MoFEMErrorCode blockModesImpl(int block_id, np::ndarray coords, np::ndarray centroid, np::ndarray bbodx, np::ndarray &o_ptr)
Internal implementation for topology mode generation.
bp::object objectNumberOfModes
Python function: numberOfModes(block_id) -> int.
MoFEMErrorCode evalObjectiveGradientStrain(MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< MatrixDouble > o_ptr)
Compute gradient of objective function with respect to strain.
bp::object objectiveGradientStress
Python function: ∂f/∂σ(coords,u,stress,strain) -> gradient
ObjectiveFunctionDataImpl()=default
MoFEMErrorCode initPython(const std::string py_file)
Initialize Python interpreter and load objective function script.
MoFEMErrorCode blockModes(int block_id, MatrixDouble &coords, std::array< double, 3 > ¢roid, std::array< double, 6 > &bbodx, MatrixDouble &o_ptr)
Define spatial topology modes for design optimization.
bp::object objectiveGradientStrain
Python function: ∂f/∂ε(coords,u,stress,strain) -> gradient.
MoFEMErrorCode evalObjectiveGradientU(MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< MatrixDouble > o_ptr)
Compute gradient of objective function with respect to displacement.
virtual ~ObjectiveFunctionDataImpl()=default
MoFEMErrorCode evalObjectiveGradientStress(MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< MatrixDouble > o_ptr)
Compute gradient of objective function with respect to stress.
bp::object objectiveFunction
Python function: f(coords,u,stress,strain) -> objective.
MoFEMErrorCode evalObjectiveFunction(MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< VectorDouble > o_ptr)
Evaluate objective function at current state.
bp::object objectiveGradientU
Python function: ∂f/∂u(coords,u,stress,strain) -> gradient.
MatrixDouble copyToFull(MatrixDouble &s)
Convert symmetric tensor storage to full matrix format.
np::ndarray convertToNumPy(std::vector< double > &data, int rows, int nb_gauss_pts)
Convert std::vector to NumPy array for Python interface.
MoFEMErrorCode objectiveGradientUImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
Internal implementation for displacement gradient computation.
MoFEMErrorCode objectiveGradientStressImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
Internal implementation for stress gradient computation.
void copyToSymmetric(double *ptr, MatrixDouble &s)
Convert full matrix to symmetric tensor storage format
bp::object objectBlockModes
Python function: blockModes(block_id,coords,centroid,bbox) -> modes.
bp::object mainNamespace
Main Python namespace for script execution.
MoFEMErrorCode objectiveFunctionImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
Internal implementation for objective function evaluation.
Abstract interface for Python-defined objective functions.
virtual ~ObjectiveFunctionData()=default
virtual MoFEMErrorCode evalObjectiveGradientStress(MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< MatrixDouble > o_ptr)=0
Evaluate gradient of objective function w.r.t. stress: ∂f/∂σ
virtual MoFEMErrorCode evalObjectiveGradientStrain(MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< MatrixDouble > o_ptr)=0
Evaluate gradient of objective function w.r.t. strain: ∂f/∂ε
virtual MoFEMErrorCode numberOfModes(int block_id, int &modes)=0
Return number of topology optimization modes for given block.
virtual MoFEMErrorCode blockModes(int block_id, MatrixDouble &coords, std::array< double, 3 > ¢roid, std::array< double, 6 > &bbox, MatrixDouble &o_ptr)=0
Define spatial modes for topology optimization.
virtual MoFEMErrorCode evalObjectiveFunction(MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< VectorDouble > o_ptr)=0
Evaluate objective function f(coords, u, stress, strain)
boost::shared_ptr< MatrixDouble > diffJac
boost::shared_ptr< MatrixDouble > jac
OpAdJointGradTimesSymTensor(const std::string field_name, boost::shared_ptr< HookeOps::CommonData > comm_ptr, boost::shared_ptr< MatrixDouble > jac, boost::shared_ptr< MatrixDouble > diff_jac, boost::shared_ptr< VectorDouble > cof_vals)
boost::shared_ptr< HookeOps::CommonData > commPtr
boost::shared_ptr< VectorDouble > cofVals
Forward declaration of operator for gradient times symmetric tensor operations.
Operator for computing objective function and gradients in topology optimization.
boost::shared_ptr< MatrixDouble > dGradPtr
boost::shared_ptr< double > globObjectiveGradPtr
boost::shared_ptr< ObjectiveFunctionData > pythonPtr
boost::shared_ptr< double > globObjectivePtr
ForcesAndSourcesCore::UserDataOperator OP
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Compute objective function contributions at element level.
boost::shared_ptr< HookeOps::CommonData > commPtr
boost::shared_ptr< VectorDouble > cofVals
boost::shared_ptr< MatrixDouble > uPtr
OpAdJointObjective(boost::shared_ptr< ObjectiveFunctionData > python_ptr, boost::shared_ptr< HookeOps::CommonData > comm_ptr, boost::shared_ptr< MatrixDouble > jac_ptr, boost::shared_ptr< MatrixDouble > diff_jac, boost::shared_ptr< VectorDouble > cof_vals, boost::shared_ptr< MatrixDouble > d_grad_ptr, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< double > glob_objective_ptr, boost::shared_ptr< double > glob_objective_grad_ptr)
boost::shared_ptr< MatrixDouble > diffJacPtr
boost::shared_ptr< MatrixDouble > jacPtr
OpAdJointTestOp(const std::string field_name, boost::shared_ptr< HookeOps::CommonData > comm_ptr)
boost::shared_ptr< HookeOps::CommonData > commPtr
PipelineManager::ElementsAndOpsByDim< 2 >::FaceSideEle SideEle
PipelineManager::ElementsAndOpsByDim< 3 >::FaceSideEle SideEle