#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
IntegrationType::GAUSS;
};
};
template <int SPACE_DIM, IntegrationType I, typename OpBase>
template <int SPACE_DIM, IntegrationType I, typename OpBase>
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;
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;
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;
std::array<double, 3> ¢roid,
std::array<double, 6> &bbox,
MatrixDouble &o_ptr) = 0;
};
boost::shared_ptr<ObjectiveFunctionData>
const std::string block_name, int dim);
private:
Vec objective_function_gradient,
Vec adjoint_vector);
boost::shared_ptr<ObjectiveFunctionData>
pythonPtr;
std::vector<SmartPetscObj<Vec>>
modeVecs;
};
char objective_function_file_name[255] = "objective_function.py";
PETSC_NULLPTR, PETSC_NULLPTR, "-objective_function",
objective_function_file_name, 255, PETSC_NULLPTR);
auto file_exists = [](std::string myfile) {
std::ifstream file(myfile.c_str());
if (file) {
return true;
}
return false;
};
if (!file_exists(objective_function_file_name)) {
MOFEM_LOG(
"WORLD", Sev::error) <<
"Objective function file NOT found: "
<< objective_function_file_name;
}
auto create_vec_modes = [&](auto block_name) {
std::regex((boost::format("%s(.*)") % block_name).str())
);
int nb_total_modes = 0;
for (auto &bc : bcs) {
auto id = bc->getMeshsetId();
int nb_modes;
nb_total_modes += nb_modes;
}
<< "Total number of modes to apply: " << nb_total_modes;
};
auto get_modes_bounding_boxes = [&](auto block_name) {
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto &bc : bcs) {
auto meshset = bc->getMeshset();
std::vector<double> x(verts.size());
std::vector<double> y(verts.size());
std::vector<double> z(verts.size());
std::array<double, 3> centroid = {0, 0, 0};
for (
int i = 0;
i != verts.size(); ++
i) {
}
MPI_Allreduce(MPI_IN_PLACE, centroid.data(), 3, MPI_DOUBLE, MPI_SUM,
int nb_vertex = verts.size();
MPI_Allreduce(MPI_IN_PLACE, &nb_vertex, 1, MPI_INT, MPI_SUM,
if (nb_vertex) {
centroid[0] /= nb_vertex;
centroid[1] /= nb_vertex;
centroid[2] /= nb_vertex;
}
std::array<double, 6> bbox = {centroid[0], centroid[1], centroid[2],
centroid[0], centroid[1], centroid[2]};
for (
int i = 0;
i != verts.size(); ++
i) {
bbox[0] = std::min(bbox[0], x[
i]);
bbox[1] = std::min(bbox[1], y[
i]);
bbox[2] = std::min(bbox[2], z[
i]);
bbox[3] = std::max(bbox[3], x[
i]);
bbox[4] = std::max(bbox[4], y[
i]);
bbox[5] = std::max(bbox[5], z[
i]);
}
MPI_Allreduce(MPI_IN_PLACE, &bbox[0], 3, MPI_DOUBLE, MPI_MIN,
MPI_Allreduce(MPI_IN_PLACE, &bbox[3], 3, MPI_DOUBLE, MPI_MAX,
<< "Block: " << bc->getName() << " centroid: " << centroid[0] << " "
<< centroid[1] << " " << centroid[2];
<< "Block: " << bc->getName() << " bbox: " << bbox[0] << " "
<< bbox[1] << " " << bbox[2] << " " << bbox[3] << " " << bbox[4]
<< " " << bbox[5];
}
};
auto eval_objective_and_gradient = [](Tao tao, Vec x, PetscReal *
f, Vec
g,
void *ctx) -> PetscErrorCode {
int iter;
CHKERR TaoGetIterationNumber(tao, &iter);
auto set_geometry = [&](Vec x) {
VecScatter ctx;
Vec x_local;
CHKERR VecScatterCreateToAll(x, &ctx, &x_local);
CHKERR VecScatterBegin(ctx, x, x_local, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecScatterEnd(ctx, x, x_local, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecScatterDestroy(&ctx);
CHKERR VecCopy(ex_ptr->initialGeometry, current_geometry);
CHKERR VecGetArrayRead(x_local, &
a);
for (auto &mode_vec : ex_ptr->modeVecs) {
<< "Adding mode with coeff: " << *coeff;
CHKERR VecAXPY(current_geometry, (*coeff), mode_vec);
++coeff;
}
CHKERR VecRestoreArrayRead(x_local, &
a);
CHKERR VecGhostUpdateBegin(current_geometry, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(current_geometry, INSERT_VALUES,
SCATTER_FORWARD);
->setOtherLocalGhostVector("ADJOINT", "ADJOINT_FIELD", "GEOMETRY",
INSERT_VALUES, SCATTER_REVERSE);
};
CHKERR KSPReset(ex_ptr->kspElastic);
CHKERR ex_ptr->solveElastic();
CHKERR ex_ptr->calculateGradient(f,
g, adjoint_vector);
CHKERR ex_ptr->postprocessElastic(iter, adjoint_vector);
};
CHKERR create_vec_modes(
"OPTIMISE");
CHKERR get_modes_bounding_boxes(
"OPTIMISE");
CHKERR TaoSetType(tao, TAOLMVM);
PETSC_DECIDE);
(void *)this);
INSERT_VALUES, SCATTER_FORWARD);
CHKERR TaoSetSolution(tao, x0);
CHKERR TaoSetFromOptions(tao);
}
}
enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
PetscInt choice_base_value = AINSWORTH;
LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
switch (choice_base_value) {
case AINSWORTH:
<< "Set AINSWORTH_LEGENDRE_BASE for displacements";
break;
case DEMKOWICZ:
<< "Set DEMKOWICZ_JACOBI_BASE for displacements";
break;
default:
break;
}
PETSC_NULLPTR);
auto project_ho_geometry = [&]() {
};
}
auto create_adjoint_dm = [&]() {
auto add_field = [&]() {
"ADJOINT_FIELD");
for (
auto tt = MBEDGE; tt <= moab::CN::TypeDimensionMap[
SPACE_DIM].second;
++tt)
"ADJOINT_FIELD", 1);
};
auto add_adjoint_fe_impl = [&]() {
"ADJOINT_FIELD");
"ADJOINT_FIELD");
"ADJOINT_FIELD");
"GEOMETRY");
"ADJOINT_FIELD");
"ADJOINT_FIELD");
"ADJOINT_FIELD");
"GEOMETRY");
auto block_name = "OPTIMISE";
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto bc : bcs) {
bc->getMeshset(),
SPACE_DIM - 1,
"ADJOINT_BOUNDARY_FE");
}
simple->getBitRefLevelMask());
};
auto set_adjoint_dm_imp = [&]() {
simple->getBitRefLevelMask());
CHKERR DMSetFromOptions(adjoint_dm);
PETSC_TRUE;
PETSC_FALSE;
};
return adjoint_dm;
};
}
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_X",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Y",
"U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Z",
"U", 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_ALL", "U", 0, 3);
simple->getProblemName(),
"U");
}
boost::make_shared<Range>(opt_ents));
boost::make_shared<Range>(opt_ents));
CHKERR DMSetUp(subset_dm_bdy);
CHKERR DMSetUp(subset_dm_domain);
auto remove_dofs = [&]() {
std::array<Range, 3> remove_dim_ents;
remove_dim_ents[0] =
remove_dim_ents[1] =
remove_dim_ents[2] =
for (
int d = 0;
d != 3; ++
d) {
<< "Removing topology modes on block OPT_REMOVE_" << (char)('X' + d)
<<
" with " << remove_dim_ents[
d].size() <<
" entities";
}
true);
CHKERR skin.find_skin(0, body_ents,
false, boundary_ents);
for (
int d = 0;
d != 3; ++
d) {
boundary_ents = subtract(boundary_ents, remove_dim_ents[d]);
}
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(boundary_ents,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr);
if (d >= 0) {
moab::Interface::UNION);
boundary_ents.merge(ents);
} else {
boundary_ents.merge(verts);
}
boundary_ents);
}
boundary_ents.merge(opt_ents);
"SUBSET_DOMAIN", "ADJOINT_FIELD", boundary_ents);
for (
int d = 0;
d != 3; ++
d) {
"SUBSET_DOMAIN", "ADJOINT_FIELD", remove_dim_ents[d], d, d);
}
boundary_ents);
}
};
auto get_lhs_fe = [&]() {
auto fe_lhs = boost::make_shared<BoundaryEle>(
mField);
fe_lhs->getRuleHook = [](int, int, int p_data) {
return 2 * p_data + p_data - 1;
};
auto &pip = fe_lhs->getOpPtrVector();
"GEOMETRY");
pip.push_back(
new OpMass(
"ADJOINT_FIELD",
"ADJOINT_FIELD",
[](double, double, double) { return 1.; }));
return fe_lhs;
};
auto get_rhs_fe = [&]() {
auto fe_rhs = boost::make_shared<BoundaryEle>(
mField);
fe_rhs->getRuleHook = [](int, int, int p_data) {
return 2 * p_data + p_data - 1;
};
auto &pip = fe_rhs->getOpPtrVector();
"GEOMETRY");
return fe_rhs;
};
auto block_name = "OPTIMISE";
std::regex((boost::format("%s(.*)") % block_name).str())
);
}
struct OpMode :
public OP {
OpMode(const std::string name,
boost::shared_ptr<ObjectiveFunctionData> python_ptr, int id,
std::vector<std::array<double, 3>> mode_centroids,
std::vector<std::array<double, 6>> mode_bboxes, int block_counter,
int mode_counter, boost::shared_ptr<Range> range = nullptr)
:
OP(name, name,
OP::OPROW, range), pythonPtr(python_ptr), iD(id),
modeVecs(mode_vecs), modeCentroids(mode_centroids),
modeBboxes(mode_bboxes), blockCounter(block_counter),
modeCounter(mode_counter) {}
if (OP::entsPtr) {
if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
}
if (!nb_rows) {
}
auto nb_base_functions = data.
getN().size2();
CHKERR pythonPtr->blockModes(iD, OP::getCoordsAtGaussPts(),
modeCentroids[blockCounter],
modeBboxes[blockCounter], blockModes);
auto nb_integration_pts = getGaussPts().size2();
if (blockModes.size2() != 3 * nb_integration_pts) {
<< "Number of modes does not match number of integration points: "
<< blockModes.size2() << "!=" << 3 * nb_integration_pts;
}
VectorDouble nf(nb_rows);
int nb_modes = blockModes.size1();
for (auto mode = 0; mode != nb_modes; ++mode) {
nf.clear();
auto t_mode = getFTensor1FromPtr<3>(&blockModes(mode, 0));
const double vol = OP::getMeasure();
auto t_w = OP::getFTensor0IntegrationWeight();
for (int gg = 0; gg != nb_integration_pts; gg++) {
const double alpha = t_w * vol;
auto t_nf = getFTensor1FromPtr<SPACE_DIM>(nf.data().data());
int rr = 0;
t_nf(
i) += alpha * t_base * t_mode(
i);
++t_base;
++t_nf;
}
for (; rr < nb_base_functions; ++rr)
++t_base;
++t_w;
++t_mode;
}
Vec vec = modeVecs[modeCounter + mode];
auto *nf_data = nf.data().data();
}
}
private:
boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
MatrixDouble blockModes;
std::vector<std::array<double, 3>> modeCentroids;
std::vector<std::array<double, 6>> modeBboxes;
int iD;
std::vector<SmartPetscObj<Vec>> modeVecs;
int blockCounter;
int modeCounter;
};
auto solve_bdy = [&]() {
auto fe_lhs = get_lhs_fe();
auto fe_rhs = get_rhs_fe();
int block_counter = 0;
int mode_counter = 0;
for (auto &bc : bcs) {
auto id = bc->getMeshsetId();
true);
auto range = boost::make_shared<Range>(ents);
auto &pip_rhs = fe_rhs->getOpPtrVector();
mode_counter, range));
fe_rhs);
pip_rhs.pop_back();
int nb_modes;
++block_counter;
mode_counter += nb_modes;
<< "Setting mode block block: " << bc->getName()
<< " with ID: " << bc->getMeshsetId()
<< " total modes: " << mode_counter;
}
CHKERR VecGhostUpdateBegin(
v, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(
v, ADD_VALUES, SCATTER_REVERSE);
}
fe_lhs);
CHKERR MatAssemblyBegin(
M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
M, MAT_FINAL_ASSEMBLY);
CHKERR KSPSetFromOptions(solver);
}
CHKERR VecGhostUpdateBegin(
v, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
v, INSERT_VALUES, SCATTER_FORWARD);
}
};
auto get_elastic_fe_lhs = [&]() {
auto fe = boost::make_shared<DomainEle>(
mField);
fe->getRuleHook = [](int, int, int p_data) {
return 2 * p_data + p_data - 1;
};
auto &pip = fe->getOpPtrVector();
"GEOMETRY");
mField, pip,
"ADJOINT_FIELD",
"MAT_ADJOINT", Sev::noisy);
return fe;
};
auto get_elastic_fe_rhs = [&]() {
auto fe = boost::make_shared<DomainEle>(
mField);
fe->getRuleHook = [](int, int, int p_data) {
return 2 * p_data + p_data - 1;
};
auto &pip = fe->getOpPtrVector();
"GEOMETRY");
mField, pip,
"ADJOINT_FIELD",
"MAT_ADJOINT", Sev::noisy);
return fe;
};
auto adjoint_gradient_postprocess = [&](auto mode) {
auto post_proc_mesh = boost::make_shared<moab::Core>();
auto post_proc_begin =
boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
mField,
post_proc_mesh);
auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
auto geom_vec = boost::make_shared<MatrixDouble>();
auto post_proc_fe =
boost::make_shared<PostProcEleDomain>(
mField, post_proc_mesh);
post_proc_fe->getOpPtrVector(), {H1}, "GEOMETRY");
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{},
{{"MODE", geom_vec}},
{},
{}
)
);
post_proc_begin->getFEMethod());
post_proc_fe);
post_proc_begin->getFEMethod());
CHKERR post_proc_end->writeFile(
"mode_" + std::to_string(mode) +
".h5m");
};
auto solve_domain = [&]() {
auto fe_lhs = get_elastic_fe_lhs();
auto fe_rhs = get_elastic_fe_rhs();
fe_lhs);
CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
CHKERR KSPSetOperators(solver, M, M);
CHKERR KSPSetFromOptions(solver);
int mode_counter = 0;
for (auto &f : modeVecs) {
SCATTER_REVERSE);
fe_rhs);
CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateBegin(
v, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
v, INSERT_VALUES, SCATTER_FORWARD);
SCATTER_REVERSE);
SCATTER_FORWARD);
++mode_counter;
}
};
for (
int i = 0;
i < modeVecs.size(); ++
i) {
CHKERR adjoint_gradient_postprocess(
i);
}
}
};
};
CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
pip->getOpDomainLhsPipeline(), {H1}, "GEOMETRY");
pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
pip->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
pip->getOpBoundaryLhsPipeline(), {NOSPACE}, "GEOMETRY");
mField, pip->getOpDomainLhsPipeline(),
"U",
"MAT_ELASTIC", Sev::noisy);
mField, pip->getOpDomainRhsPipeline(),
"U",
"MAT_ELASTIC", Sev::noisy);
CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
pip->getOpDomainRhsPipeline(),
mField,
"U", Sev::inform);
CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
pip->getOpBoundaryRhsPipeline(),
mField,
"U", -1, Sev::inform);
CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
pip->getOpBoundaryLhsPipeline(),
mField,
"U", Sev::noisy);
}
auto set_essential_bc = [&]() {
auto pre_proc_rhs = boost::make_shared<FEMethod>();
auto post_proc_rhs = boost::make_shared<FEMethod>();
auto post_proc_lhs = boost::make_shared<FEMethod>();
auto get_pre_proc_hook = [&]() {
{});
};
pre_proc_rhs->preProcessHook = get_pre_proc_hook();
auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
post_proc_rhs, 1.)();
};
auto get_post_proc_hook_lhs = [this, post_proc_lhs]() {
post_proc_lhs, 1.)();
};
post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
post_proc_lhs->postProcessHook = get_post_proc_hook_lhs;
ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
};
auto setup_and_solve = [&](auto solver) {
BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
MOFEM_LOG(
"TIMER", Sev::noisy) <<
"KSPSetUp";
MOFEM_LOG(
"TIMER", Sev::noisy) <<
"KSPSetUp <= Done";
MOFEM_LOG(
"TIMER", Sev::noisy) <<
"KSPSolve";
CHKERR KSPSolve(solver, f, d);
MOFEM_LOG(
"TIMER", Sev::noisy) <<
"KSPSolve <= Done";
};
CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
auto evaluate_field_at_the_point = [&]() {
int coords_dim = 3;
std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
field_eval_coords.data(), &coords_dim,
auto field_eval_data =
->buildTree<SPACE_DIM>(field_eval_data,
simple->getDomainFEName());
field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
auto no_rule = [](int, int, int) { return -1; };
auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
field_eval_fe_ptr->getRuleHook = no_rule;
field_eval_fe_ptr->getOpPtrVector().push_back(
->evalFEAtThePoint<SPACE_DIM>(
field_eval_coords.data(), 1e-12,
simple->getProblemName(),
simple->getDomainFEName(), field_eval_data,
<< "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
else
<< "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
<< " U_Z: " << t_disp(2);
}
}
};
CHKERR evaluate_field_at_the_point();
}
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pip->getDomainRhsFE().reset();
pip->getDomainLhsFE().reset();
pip->getBoundaryRhsFE().reset();
pip->getBoundaryLhsFE().reset();
auto post_proc_mesh = boost::make_shared<moab::Core>();
auto post_proc_begin = boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
auto calculate_stress_ops = [&](auto &pip) {
mField, pip,
"U",
"MAT_ELASTIC", Sev::noisy);
auto u_ptr = boost::make_shared<MatrixDouble>();
auto x_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
DataMapMat vec_map = {{"U", u_ptr}, {"GEOMETRY", x_ptr}};
DataMapMat sym_map = {{"STRAIN", common_ptr->getMatStrain()},
{"STRESS", common_ptr->getMatCauchyStress()}};
if (adjoint_vector) {
auto adjoint_ptr = boost::make_shared<MatrixDouble>();
"U", adjoint_ptr, adjoint_vector));
vec_map["ADJOINT"] = adjoint_ptr;
}
return boost::make_tuple(u_ptr, x_ptr, common_ptr->getMatStrain(),
common_ptr->getMatCauchyStress(), vec_map,
sym_map);
};
auto get_tag_id_on_pmesh = [&](bool post_proc_skin) {
int def_val_int = 0;
"MAT_ELASTIC", 1, MB_TYPE_INTEGER, tag_mat,
MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
auto meshset_vec_ptr =
std::regex((boost::format("%s(.*)") % "MAT_ELASTIC").str()));
std::unique_ptr<Skinner> skin_ptr;
if (post_proc_skin) {
auto boundary_meshset =
simple->getBoundaryMeshSet();
skin_ents, true);
}
for (
auto m : meshset_vec_ptr) {
true);
int const id =
m->getMeshsetId();
ents_3d = ents_3d.subset_by_dimension(
SPACE_DIM);
if (post_proc_skin) {
CHKERR skin_ptr->find_skin(0, ents_3d,
false, skin_faces);
ents_3d = intersect(skin_ents, skin_faces);
}
}
return tag_mat;
};
auto post_proc_domain = [&](auto post_proc_mesh) {
auto post_proc_fe =
boost::make_shared<PostProcEleDomain>(
mField, post_proc_mesh);
auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr, vec_map, sym_map] =
calculate_stress_ops(post_proc_fe->getOpPtrVector());
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{},
vec_map,
{},
sym_map
)
);
post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(false)});
return post_proc_fe;
};
auto post_proc_boundary = [&](auto post_proc_mesh) {
auto post_proc_fe =
boost::make_shared<PostProcEleBdy>(
mField, post_proc_mesh);
post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
auto op_loop_side =
auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr, vec_map, sym_map] =
calculate_stress_ops(op_loop_side->getOpPtrVector());
post_proc_fe->getOpPtrVector().push_back(op_loop_side);
auto mat_traction_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
mat_traction_ptr));
vec_map["T"] = mat_traction_ptr;
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{},
vec_map,
{},
sym_map
)
);
post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(true)});
return post_proc_fe;
};
PetscBool post_proc_skin_only = PETSC_FALSE;
post_proc_skin_only = PETSC_TRUE;
&post_proc_skin_only, PETSC_NULLPTR);
}
if (post_proc_skin_only == PETSC_FALSE) {
pip->getDomainRhsFE() = post_proc_domain(post_proc_mesh);
} else {
pip->getBoundaryRhsFE() = post_proc_boundary(post_proc_mesh);
}
post_proc_begin->getFEMethod());
CHKERR pip->loopFiniteElements();
post_proc_end->getFEMethod());
CHKERR post_proc_end->writeFile(
"out_elastic_" + std::to_string(iter) +
".h5m");
}
template <int SPACE_DIM>
boost::shared_ptr<HookeOps::CommonData> comm_ptr)
protected:
boost::shared_ptr<HookeOps::CommonData> commPtr;
};
template <int SPACE_DIM>
const double vol = OP::getMeasure();
auto t_w = OP::getFTensor0IntegrationWeight();
auto t_cauchy_stress =
getFTensor2SymmetricFromMat<SPACE_DIM>(*(commPtr->getMatCauchyStress()));
for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
const double alpha = t_w * vol;
auto t_nf = OP::template getNf<SPACE_DIM>();
int rr = 0;
t_nf(
j) += alpha * t_row_grad(
i) * t_cauchy_stress(
i,
j);
++t_row_grad;
++t_nf;
}
for (; rr < OP::nbRowBaseFunctions; ++rr) {
++t_row_grad;
}
++t_cauchy_stress;
++t_w;
}
}
template <int SPACE_DIM>
boost::shared_ptr<HookeOps::CommonData> comm_ptr,
boost::shared_ptr<MatrixDouble> jac,
boost::shared_ptr<MatrixDouble> diff_jac,
boost::shared_ptr<VectorDouble> cof_vals)
jac(jac), diffJac(diff_jac), cofVals(cof_vals) {}
protected:
boost::shared_ptr<HookeOps::CommonData> commPtr;
boost::shared_ptr<MatrixDouble> jac;
boost::shared_ptr<MatrixDouble> diffJac;
boost::shared_ptr<VectorDouble> cofVals;
};
t_diff(0, 0, 0, 0) = 1;
t_diff(1, 1, 1, 1) = 1;
t_diff(1, 0, 1, 0) = 0.5;
t_diff(1, 0, 0, 1) = 0.5;
t_diff(0, 1, 0, 1) = 0.5;
t_diff(0, 1, 1, 0) = 0.5;
if constexpr (DIM == 3) {
t_diff(2, 2, 2, 2) = 1;
t_diff(2, 0, 2, 0) = 0.5;
t_diff(2, 0, 0, 2) = 0.5;
t_diff(0, 2, 0, 2) = 0.5;
t_diff(0, 2, 2, 0) = 0.5;
t_diff(2, 1, 2, 1) = 0.5;
t_diff(2, 1, 1, 2) = 0.5;
t_diff(1, 2, 1, 2) = 0.5;
t_diff(1, 2, 2, 1) = 0.5;
}
return t_diff;
};
template <int SPACE_DIM>
const double vol = OP::getMeasure();
auto t_w = OP::getFTensor0IntegrationWeight();
auto t_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(jac));
auto t_diff_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(diffJac));
auto t_cof = getFTensor0FromVec(*(cofVals));
auto t_grad_u =
getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(commPtr->matGradPtr));
auto t_cauchy_stress =
getFTensor2SymmetricFromMat<SPACE_DIM>(*(commPtr->getMatCauchyStress()));
auto t_D =
getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*(commPtr->matDPtr));
for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
const double alpha = t_w * vol;
auto t_det = determinantTensor(t_jac);
CHKERR invertTensor(t_jac, t_det, t_inv_jac);
-(t_inv_jac(
i,
I) * t_diff_jac(
I,
J)) * t_inv_jac(
J,
j);
t_diff_grad(
i,
j) = t_grad_u(
i,
k) * t_diff_inv_jac(
k,
j);
t_diff_strain(
i,
j) = t_diff_symm(
i,
j,
k,
l) * t_diff_grad(
k,
l);
t_diff_stress(
i,
j) = t_D(
i,
j,
k,
l) * t_diff_strain(
k,
l);
auto t_nf = OP::template getNf<SPACE_DIM>();
int rr = 0;
t_diff_row_grad(
k) = t_row_grad(
j) * t_diff_inv_jac(
j,
k);
t_nf(
j) += alpha * t_diff_row_grad(
i) * t_cauchy_stress(
i,
j);
t_nf(
j) += (alpha * t_cof) * t_row_grad(
i) * t_cauchy_stress(
i,
j);
t_nf(
j) += alpha * t_row_grad(
i) * t_diff_stress(
i,
j);
++t_row_grad;
++t_nf;
}
for (; rr < OP::nbRowBaseFunctions; ++rr) {
++t_row_grad;
}
++t_grad_u;
++t_cauchy_stress;
++t_jac;
++t_diff_jac;
++t_cof;
++t_w;
}
}
using OP = ForcesAndSourcesCore::UserDataOperator;
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)
auto nb_gauss_pts = getGaussPts().size2();
auto objective_ptr = boost::make_shared<VectorDouble>(nb_gauss_pts);
auto objective_dstress =
boost::make_shared<MatrixDouble>(symm_size, nb_gauss_pts);
auto objective_dstrain =
boost::make_shared<MatrixDouble>(symm_size, nb_gauss_pts);
auto obj_grad = boost::make_shared<MatrixDouble>(
SPACE_DIM, nb_gauss_pts);
auto evaluate_python = [&]() {
auto &coords = OP::getCoordsAtGaussPts();
objective_ptr);
objective_dstress);
objective_dstrain);
auto t_grad_u =
auto t_D =
auto t_obj_dstress =
auto t_obj_dstrain =
auto vol = OP::getMeasure();
auto t_w = getFTensor0IntegrationWeight();
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
-(t_inv_jac(
i,
I) * t_diff_jac(
I,
J)) * t_inv_jac(
J,
j);
t_diff_grad(
i,
j) = t_grad_u(
i,
k) * t_diff_inv_jac(
k,
j);
t_d_strain(
i,
j) = t_diff_symm(
i,
j,
k,
l) * (
+
);
auto alpha = t_w * vol;
(*globObjectivePtr) += alpha * t_obj;
(*globObjectiveGradPtr) +=
alpha *
(
t_obj_dstress(
i,
j) * (t_D(
i,
j,
k,
l) * t_d_strain(
k,
l))
+
t_obj_dstrain(
i,
j) * t_d_strain(
i,
j)
+
t_obj * t_cof
);
++t_w;
++t_jac;
++t_diff_jac;
++t_cof;
++t_obj;
++t_obj_dstress;
++t_obj_dstrain;
++t_grad_u;
++t_d_grad;
}
};
}
private:
boost::shared_ptr<ObjectiveFunctionData>
pythonPtr;
boost::shared_ptr<HookeOps::CommonData>
commPtr;
boost::shared_ptr<MatrixDouble>
jacPtr;
boost::shared_ptr<VectorDouble>
cofVals;
boost::shared_ptr<MatrixDouble>
dGradPtr;
boost::shared_ptr<MatrixDouble>
uPtr;
};
Vec objective_function_gradient,
Vec adjoint_vector) {
auto get_essential_fe = [this]() {
auto post_proc_rhs = boost::make_shared<FEMethod>();
auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
post_proc_rhs, 0)();
};
post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
return post_proc_rhs;
};
auto get_fd_direvative_fe = [&]() {
auto fe = boost::make_shared<DomainEle>(
mField);
fe->getRuleHook = [](int, int, int p_data) {
return 2 * p_data + p_data - 1;
};
auto &pip = fe->getOpPtrVector();
constexpr bool debug =
false;
mField, pip,
"U",
"MAT_ELASTIC", Sev::noisy);
pip.push_back(
} else {
mField, pip,
"U",
"MAT_ELASTIC", Sev::noisy);
}
return fe;
};
auto calulate_fd_residual = [&](
auto eps,
auto diff_vec,
auto fd_vec) {
constexpr bool debug =
false;
double nrm2 = 0.0;
auto norm2_field = [&](const double val) {
nrm2 += val * val;
return val;
};
MPI_Allreduce(MPI_IN_PLACE, &nrm2, 1, MPI_DOUBLE, MPI_SUM,
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Geometry norm: " << sqrt(nrm2);
};
initial_current_geometry, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecAssemblyBegin(initial_current_geometry);
CHKERR VecAssemblyEnd(initial_current_geometry);
auto perturb_geometry = [&](
auto eps,
auto diff_vec) {
CHKERR VecCopy(initial_current_geometry, current_geometry);
CHKERR VecAXPY(current_geometry,
eps, diff_vec);
current_geometry, INSERT_VALUES, SCATTER_REVERSE);
};
auto fe = get_fd_direvative_fe();
auto calc_impl = [&](
auto f,
auto eps) {
simple->getDomainFEName(), fe);
CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
auto post_proc_rhs = get_essential_fe();
post_proc_rhs.get());
};
CHKERR VecWAXPY(fd_vec, -1.0, fm, fp);
initial_current_geometry, INSERT_VALUES, SCATTER_REVERSE);
};
auto get_direvative_fe = [&](auto diff_vec) {
auto fe_adjoint = boost::make_shared<DomainEle>(
mField);
fe_adjoint->getRuleHook = [](int, int, int p_data) {
return 2 * p_data + p_data - 1;
};
auto &pip = fe_adjoint->getOpPtrVector();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
auto cof_ptr = boost::make_shared<VectorDouble>();
using OpCoFactor =
"GEOMETRY", jac_ptr));
"U", diff_jac_ptr,
diff_vec));
pip.push_back(new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
mField, pip,
"U",
"MAT_ELASTIC", Sev::noisy);
"U", common_ptr, jac_ptr, diff_jac_ptr, cof_ptr));
return fe_adjoint;
};
auto get_objective_fe = [&](auto diff_vec, auto grad_vec,
auto glob_objective_ptr,
auto glob_objective_grad_ptr) {
auto fe_adjoint = boost::make_shared<DomainEle>(
mField);
fe_adjoint->getRuleHook = [](int, int, int p_data) {
return 2 * p_data + p_data - 1;
};
auto &pip = fe_adjoint->getOpPtrVector();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
auto cof_ptr = boost::make_shared<VectorDouble>();
auto d_grad_ptr = boost::make_shared<MatrixDouble>();
auto u_ptr = boost::make_shared<MatrixDouble>();
using OpCoFactor =
"GEOMETRY", jac_ptr));
"U", diff_jac_ptr,
diff_vec));
pip.push_back(new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
"U", d_grad_ptr, grad_vec));
mField, pip,
"U",
"MAT_ELASTIC", Sev::noisy);
pythonPtr, common_ptr, jac_ptr, diff_jac_ptr, cof_ptr, d_grad_ptr,
u_ptr, glob_objective_ptr, glob_objective_grad_ptr));
return fe_adjoint;
};
auto adjoint_fe = get_direvative_fe(dm_diff_vec);
auto glob_objective_ptr = boost::make_shared<double>(0.0);
auto glob_objective_grad_ptr = boost::make_shared<double>(0.0);
auto objective_fe = get_objective_fe(dm_diff_vec, d, glob_objective_ptr,
glob_objective_grad_ptr);
auto set_variance_of_geometry = [&](auto mode, auto mod_vec) {
SCATTER_REVERSE);
dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
};
auto calculate_variance_internal_forces = [&](auto mode, auto mod_vec) {
CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
auto post_proc_rhs = get_essential_fe();
post_proc_rhs.get());
#ifndef NDEBUG
constexpr bool debug =
false;
double norm0;
CHKERR VecNorm(f, NORM_2, &norm0);
CHKERR calulate_fd_residual(
eps, dm_diff_vec, fd_check);
double nrm;
CHKERR VecAXPY(fd_check, -1.0, f);
CHKERR VecNorm(fd_check, NORM_2, &nrm);
<< " FD check for internal forces [ " << mode << " ]: " << nrm
<< " / " << norm0 << " ( " << (nrm / norm0) << " )";
}
#endif
};
auto calulate_variance_of_displacemnte = [&](auto mode, auto mod_vec) {
CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
};
auto calculate_variance_of_objective_function = [&](auto mode, auto mod_vec) {
*glob_objective_ptr = 0.0;
*glob_objective_grad_ptr = 0.0;
objective_fe);
CHKERR VecSetValue(objective_function_gradient, mode,
*glob_objective_grad_ptr, ADD_VALUES);
std::array<double, 2> array = {*glob_objective_ptr,
*glob_objective_grad_ptr};
MPI_Allreduce(MPI_IN_PLACE, array.data(), 2, MPI_DOUBLE, MPI_SUM,
*glob_objective_ptr = array[0];
*glob_objective_grad_ptr = array[1];
(*objective_function_value) += *glob_objective_ptr;
MOFEM_LOG(
"WORLD", Sev::verbose) <<
" Objective gradient [ " << mode
<< " ]: " << *glob_objective_grad_ptr;
};
CHKERR VecZeroEntries(objective_function_gradient);
CHKERR VecZeroEntries(adjoint_vector);
*objective_function_value = 0.0;
int mode = 0;
CHKERR set_variance_of_geometry(mode, mod_vec);
CHKERR calculate_variance_internal_forces(mode, mod_vec);
CHKERR calulate_variance_of_displacemnte(mode, mod_vec);
CHKERR calculate_variance_of_objective_function(mode, mod_vec);
CHKERR VecAXPY(adjoint_vector, *glob_objective_grad_ptr, dm_diff_vec);
++mode;
}
(*objective_function_value) /=
modeVecs.size();
<< "Objective function: " << *glob_objective_ptr;
CHKERR VecAssemblyBegin(objective_function_gradient);
CHKERR VecAssemblyEnd(objective_function_gradient);
CHKERR VecAssemblyBegin(adjoint_vector);
CHKERR VecAssemblyEnd(adjoint_vector);
CHKERR VecGhostUpdateBegin(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
}
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
Py_Initialize();
np::initialize();
const char param_file[] = "param_file.petsc";
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
core_log->add_sink(
LogManager::createSink(LogManager::getStrmSync(), "FieldEvaluator"));
LogManager::setLog("FieldEvaluator");
try {
DMType dm_name = "DMMOFEM";
DMType dm_name_mg = "DMMOFEM_MG";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
if (Py_FinalizeEx() < 0) {
exit(120);
}
}
boost::shared_ptr<MatrixDouble> u_ptr,
boost::shared_ptr<MatrixDouble> stress_ptr,
boost::shared_ptr<MatrixDouble> strain_ptr,
boost::shared_ptr<VectorDouble> o_ptr);
boost::shared_ptr<MatrixDouble> u_ptr,
boost::shared_ptr<MatrixDouble> stress_ptr,
boost::shared_ptr<MatrixDouble> strain_ptr,
boost::shared_ptr<MatrixDouble> o_ptr);
boost::shared_ptr<MatrixDouble> u_ptr,
boost::shared_ptr<MatrixDouble> stress_ptr,
boost::shared_ptr<MatrixDouble> strain_ptr,
boost::shared_ptr<MatrixDouble> o_ptr);
boost::shared_ptr<MatrixDouble> u_ptr,
boost::shared_ptr<MatrixDouble> stress_ptr,
boost::shared_ptr<MatrixDouble> strain_ptr,
boost::shared_ptr<MatrixDouble> o_ptr);
std::array<double, 3> ¢roid,
std::array<double, 6> &bbodx, MatrixDouble &o_ptr);
private:
np::ndarray coords, np::ndarray u,
np::ndarray stress, np::ndarray strain, np::ndarray &o
);
np::ndarray coords, np::ndarray u,
np::ndarray stress, np::ndarray strain, np::ndarray &o
);
np::ndarray coords, np::ndarray u,
np::ndarray stress, np::ndarray strain, np::ndarray &o
);
np::ndarray coords, np::ndarray u,
np::ndarray stress, np::ndarray strain, np::ndarray &o
);
int block_id, np::ndarray coords, np::ndarray centroid, np::ndarray bbodx,
np::ndarray &o_ptr
);
int nb_gauss_pts);
};
boost::shared_ptr<ObjectiveFunctionData>
auto ptr = boost::make_shared<ObjectiveFunctionDataImpl>();
return ptr;
}
try {
auto main_module = bp::import("__main__");
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
MatrixDouble
f(9, s.size2());
for (int ii = 0; ii != s.size2(); ++ii) {
++t_f;
++t_s;
}
};
ptr + 0 * s.size2(), ptr + 1 * s.size2(),
ptr + 2 * s.size2(), ptr + 3 * s.size2(),
ptr + 4 * s.size2(), ptr + 5 * s.size2(),
ptr + 6 * s.size2(), ptr + 7 * s.size2(),
ptr + 8 * s.size2()
};
for (int ii = 0; ii != s.size2(); ++ii) {
t_s(
i,
j) = (t_f(
i,
j) || t_f(
j,
i)) / 2.0;
++t_f;
++t_s;
}
}
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) {
try {
auto np_coords =
auto np_u =
convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
auto np_stress =
convertToNumPy(full_stress.data(), full_stress.size1(),
full_stress.size2());
auto np_strain =
convertToNumPy(full_strain.data(), full_strain.size1(),
full_strain.size2());
np::ndarray np_output = np::empty(bp::make_tuple(strain_ptr->size2()),
np::dtype::get_builtin<double>());
np_output);
o_ptr->resize(stress_ptr->size2(), false);
double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
std::copy(val_ptr, val_ptr + strain_ptr->size2(), o_ptr->data().begin());
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
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) {
try {
auto np_coords =
auto np_u =
convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
auto np_stress =
convertToNumPy(full_stress.data(), full_stress.size1(),
full_stress.size2());
auto np_strain =
convertToNumPy(full_strain.data(), full_strain.size1(),
full_strain.size2());
np::ndarray np_output =
np::empty(bp::make_tuple(full_strain.size1(), full_strain.size2()),
np::dtype::get_builtin<double>());
np_output);
o_ptr->resize(stress_ptr->size1(), stress_ptr->size2(), false);
double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
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) {
try {
auto np_coords =
auto np_u =
convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
auto np_stress =
convertToNumPy(full_stress.data(), full_stress.size1(),
full_stress.size2());
auto np_strain =
convertToNumPy(full_strain.data(), full_strain.size1(),
full_strain.size2());
np::ndarray np_output =
np::empty(bp::make_tuple(full_strain.size1(), full_strain.size2()),
np::dtype::get_builtin<double>());
np_output);
o_ptr->resize(strain_ptr->size1(), strain_ptr->size2(), false);
double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
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) {
try {
auto np_coords =
auto np_u =
convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
auto np_stress =
convertToNumPy(full_stress.data(), full_stress.size1(),
full_stress.size2());
auto np_strain =
convertToNumPy(full_strain.data(), full_strain.size1(),
full_strain.size2());
np::ndarray np_output =
np::empty(bp::make_tuple(u_ptr->size1(), u_ptr->size2()),
np::dtype::get_builtin<double>());
np_output);
o_ptr->resize(u_ptr->size1(), u_ptr->size2(), false);
double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
std::copy(val_ptr, val_ptr + u_ptr->size1() * u_ptr->size2(),
o_ptr->data().begin());
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
int block_id, MatrixDouble &coords, std::array<double, 3> ¢roid,
std::array<double, 6> &bbodx, MatrixDouble &o_ptr) {
try {
auto np_coords =
np::ndarray np_output =
np::empty(bp::make_tuple(coords.size1(), nb_modes, 3),
np::dtype::get_builtin<double>());
np_output);
o_ptr.resize(nb_modes, coords.size1() * coords.size2(), false);
double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
std::copy(val_ptr, val_ptr + coords.size1() * coords.size2() * nb_modes,
o_ptr.data().begin());
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
np::ndarray coords, np::ndarray u,
np::ndarray stress, np::ndarray strain, np::ndarray &o
) {
try {
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
np::ndarray coords, np::ndarray u,
np::ndarray stress, np::ndarray strain, np::ndarray &o
) {
try {
o = bp::extract<np::ndarray>(
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
np::ndarray coords, np::ndarray u,
np::ndarray stress, np::ndarray strain, np::ndarray &o
) {
try {
o = bp::extract<np::ndarray>(
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
np::ndarray coords, np::ndarray u,
np::ndarray stress, np::ndarray strain, np::ndarray &o
) {
try {
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
int &modes) {
try {
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
np::ndarray coords,
np::ndarray centroid,
np::ndarray bbodx,
np::ndarray &o) {
try {
o = bp::extract<np::ndarray>(
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
inline np::ndarray
int nb_gauss_pts) {
auto dtype = np::dtype::get_builtin<double>();
auto size = bp::make_tuple(rows, nb_gauss_pts);
auto stride = bp::make_tuple(nb_gauss_pts * sizeof(double), sizeof(double));
return (np::from_data(data.data(), dtype, size, stride, bp::object()));
}
int s) {
auto dtype = np::dtype::get_builtin<double>();
auto size = bp::make_tuple(s);
auto stride = bp::make_tuple(sizeof(double));
return (np::from_data(ptr, dtype, size, stride, bp::object()));
}
const std::string block_name, int dim) {
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto bc : bcs) {
faces, true),
"get meshset ents");
}
for (
auto dd = dim - 1;
dd >= 0; --
dd) {
if (dd >= 0) {
moab::Interface::UNION),
"get adjs");
} else {
"get verts");
}
}
};
CHKERR moab.add_entities(*out_meshset, r);
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 >)
PetscBool is_plane_strain
constexpr int SPACE_DIM
[Define dimension]
constexpr IntegrationType I
FormsIntegrators< DomainEleOp >::Assembly< A >::OpBase DomainBaseOp
[Postprocess results]
boost::shared_ptr< ObjectiveFunctionData > create_python_objective_function(std::string py_file)
Range get_range_from_block(MoFEM::Interface &m_field, const std::string block_name, int dim)
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#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.
#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 Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HookeOps::CommonData > common_ptr, Sev sev)
Assemble domain LHS K factory (LHS first overload) Initializes and pushes operators to assemble the L...
auto commonDataFactory(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HookeOps::CommonData > common_ptr, Sev sev, bool is_non_linear=false)
Factory function to create and push internal force operators for the domain RHS. (RHS first overload)...
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.
static FTensor::Ddg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim23 > getFTensor4DdgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
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)
FTensor::Index< 'M', 3 > M
constexpr IntegrationType I
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double young_modulus
Young modulus.
NaturalBC< BoundaryEleOp >::Assembly< AT >::BiLinearForm< IT > BoundaryLhsBCs
NaturalBC< BoundaryEleOp >::Assembly< AT >::LinearForm< IT > BoundaryRhsBCs
#define EXECUTABLE_DIMENSION
PetscBool do_eval_field
Evaluate field.
double poisson_ratio
Poisson ratio.
DomainRhsBCs::OpFlux< PlasticOps::DomainBCs, 1, SPACE_DIM > OpDomainRhsBCs
NaturalBC< DomainEleOp >::Assembly< AT >::LinearForm< IT > DomainRhsBCs
BoundaryRhsBCs::OpFlux< PlasticOps::BoundaryBCs, 1, SPACE_DIM > OpBoundaryRhsBCs
BoundaryLhsBCs::OpFlux< PlasticOps::BoundaryBCs, 1, SPACE_DIM > OpBoundaryLhsBCs
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
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
MoFEMErrorCode readMesh()
[Run problem]
SmartPetscObj< KSP > kspElastic
std::vector< SmartPetscObj< Vec > > modeVecs
FieldApproximationBase base
std::vector< std::array< double, 3 > > modeCentroids
MoFEMErrorCode topologyModes()
[Boundary condition]
SmartPetscObj< Vec > initialGeometry
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
[Run problem]
MoFEMErrorCode calculateGradient(PetscReal *objective_function_value, Vec objective_function_gradient, Vec adjoint_vector)
MoFEMErrorCode setupAdJoint()
[Set up problem]
MoFEM::Interface & mField
std::vector< std::array< double, 6 > > modeBBoxes
boost::shared_ptr< ObjectiveFunctionData > pythonPtr
SmartPetscObj< DM > adjointDM
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode postprocessElastic(int iter, SmartPetscObj< Vec > adjoint_vector=nullptr)
[Solve]
MoFEMErrorCode solveElastic()
[Push operators to pipeline]
boost::shared_ptr< MatrixDouble > vectorFieldPtr
Simple interface for fast problem set-up.
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 dofs 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.
Interface for managing meshsets containing materials and boundary conditions.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Get values at integration pts for tensor field rank 1, i.e. vector field.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
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.
MoFEMErrorCode objectiveGradientStrainImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
MoFEMErrorCode numberOfModes(int block_id, int &modes)
MoFEMErrorCode blockModesImpl(int block_id, np::ndarray coords, np::ndarray centroid, np::ndarray bbodx, np::ndarray &o_ptr)
bp::object objectNumberOfModes
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)
bp::object objectiveGradientStress
ObjectiveFunctionDataImpl()=default
MoFEMErrorCode initPython(const std::string py_file)
MoFEMErrorCode blockModes(int block_id, MatrixDouble &coords, std::array< double, 3 > ¢roid, std::array< double, 6 > &bbodx, MatrixDouble &o_ptr)
bp::object objectiveGradientStrain
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)
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)
bp::object objectiveFunction
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)
bp::object objectiveGradientU
MatrixDouble copyToFull(MatrixDouble &s)
np::ndarray convertToNumPy(std::vector< double > &data, int rows, int nb_gauss_pts)
Converts a std::vector<double> to a NumPy ndarray.
MoFEMErrorCode objectiveGradientUImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
MoFEMErrorCode objectiveGradientStressImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
void copyToSymmetric(double *ptr, MatrixDouble &s)
bp::object objectBlockModes
MoFEMErrorCode objectiveFunctionImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
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
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
virtual MoFEMErrorCode numberOfModes(int block_id, int &modes)=0
virtual MoFEMErrorCode blockModes(int block_id, MatrixDouble &coords, std::array< double, 3 > ¢roid, std::array< double, 6 > &bbox, MatrixDouble &o_ptr)=0
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
boost::shared_ptr< MatrixDouble > dGradPtr
boost::shared_ptr< double > globObjectiveGradPtr
boost::shared_ptr< ObjectiveFunctionData > pythonPtr
boost::shared_ptr< double > globObjectivePtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
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
