#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
0>;
GAUSS>::OpMixDivTimesScalar<SPACE_DIM>;
GAUSS>::OpBaseTimesVector<3, SPACE_DIM, 1>;
GAUSS>::OpMixDivTimesU<3, 1, SPACE_DIM>;
GAUSS>::OpBaseTimesScalar<1>;
GAUSS>::OpEssentialRhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
GAUSS>::OpEssentialLhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
1;
CHKERR moab.add_entities(*out_meshset,
r);
CHKERR moab.write_file(name.c_str(),
"VTK",
"", out_meshset->get_ptr(), 1);
};
private:
std::array<double, SPACE_DIM> fieldEvalCoords;
boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> fieldEvalData;
boost::shared_ptr<VectorDouble> scalarFieldPtr;
boost::shared_ptr<MatrixDouble> vectorFieldPtr;
boost::shared_ptr<MatrixDouble> tensorFieldPtr;
struct BlockedParameters
: public boost::enable_shared_from_this<BlockedParameters> {
double coeffExpansion;
double heatConductivity;
double heatCapacity;
inline auto getDPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &mD);
}
inline auto getCoeffExpansionPtr() {
return boost::shared_ptr<double>(shared_from_this(), &coeffExpansion);
}
inline auto getHeatConductivityPtr() {
return boost::shared_ptr<double>(shared_from_this(), &heatConductivity);
}
inline auto getHeatCapacityPtr() {
return boost::shared_ptr<double>(shared_from_this(), &heatCapacity);
}
};
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::string block_elastic_name, std::string block_thermal_name,
boost::shared_ptr<BlockedParameters> blockedParamsPtr,
Sev sev);
};
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::string block_elastic_name, std::string block_thermal_name,
boost::shared_ptr<BlockedParameters> blockedParamsPtr,
Sev sev) {
OpMatElasticBlocks(boost::shared_ptr<MatrixDouble>
m,
double bulk_modulus_K,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
"Can not get data from block");
}
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
}
}
CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
}
private:
boost::shared_ptr<MatrixDouble> matDPtr;
double bulkModulusK;
double shearModulusG;
};
double bulkModulusKDefault;
double shearModulusGDefault;
std::vector<BlockData> blockData;
std::vector<const CubitMeshSets *> meshset_vec_ptr,
for (
auto m : meshset_vec_ptr) {
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() < 2) {
"Expected that block has two attributes");
}
auto get_block_ents = [&]() {
m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
return ents;
};
blockData.push_back(
}
}
auto set_material_stiffness = [&]() {
: 1;
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
};
set_material_stiffness();
}
};
double default_bulk_modulus_K =
double default_shear_modulus_G =
pipeline.push_back(new OpMatElasticBlocks(
blockedParamsPtr->getDPtr(), default_bulk_modulus_K,
default_bulk_modulus_K, mField, sev,
(boost::format("%s(.*)") % block_elastic_name).str()
))
));
OpMatThermalBlocks(boost::shared_ptr<double> expansion_ptr,
boost::shared_ptr<double> conductivity_ptr,
boost::shared_ptr<double> capacity_ptr,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
expansionPtr(expansion_ptr), conductivityPtr(conductivity_ptr),
capacityPtr(capacity_ptr) {
"Can not get data from block");
}
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
*expansionPtr = b.expansion;
*conductivityPtr = b.conductivity;
*capacityPtr = b.capacity;
}
}
}
private:
double expansion;
double conductivity;
double capacity;
};
std::vector<BlockData> blockData;
std::vector<const CubitMeshSets *> meshset_vec_ptr,
for (
auto m : meshset_vec_ptr) {
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() < 3) {
"Expected that block has two attributes");
}
auto get_block_ents = [&]() {
m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
return ents;
};
<<
m->getName() <<
": expansion = " << block_data[0]
<< " conductivity = " << block_data[1] << " capacity "
<< block_data[2];
blockData.push_back(
{block_data[0], block_data[1], block_data[2], get_block_ents()});
}
}
boost::shared_ptr<double> expansionPtr;
boost::shared_ptr<double> conductivityPtr;
boost::shared_ptr<double> capacityPtr;
};
pipeline.push_back(new OpMatThermalBlocks(
blockedParamsPtr->getCoeffExpansionPtr(),
blockedParamsPtr->getHeatConductivityPtr(),
blockedParamsPtr->getHeatCapacityPtr(), mField, sev,
(boost::format("%s(.*)") % block_thermal_name).str()
))
));
}
}
fieldEvalCoords.data(), &coords_dim,
scalarFieldPtr = boost::make_shared<VectorDouble>();
vectorFieldPtr = boost::make_shared<MatrixDouble>();
tensorFieldPtr = boost::make_shared<MatrixDouble>();
fieldEvalData =
fieldEvalData,
simple->getDomainFEName());
} else {
fieldEvalData,
simple->getDomainFEName());
}
fieldEvalData->setEvalPoints(fieldEvalCoords.data(), 1);
auto no_rule = [](
int,
int,
int) {
return -1; };
auto field_eval_fe_ptr = fieldEvalData->feMethodPtr.lock();
field_eval_fe_ptr->getRuleHook = no_rule;
field_eval_fe_ptr->getOpPtrVector().push_back(
field_eval_fe_ptr->getOpPtrVector().push_back(
field_eval_fe_ptr->getOpPtrVector().push_back(
"U", tensorFieldPtr));
}
}
auto get_command_line_parameters = [&]() {
PETSC_NULL);
PETSC_NULL);
<<
"Reference temperature " <<
ref_temp;
};
CHKERR get_command_line_parameters();
}
auto get_skin = [&]() {
CHKERR mField.get_moab().get_entities_by_dimension(0,
SPACE_DIM, body_ents);
Skinner skin(&mField.get_moab());
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
return skin_ents;
};
auto filter_flux_blocks = [&](auto skin) {
auto remove_cubit_blocks = [&](
auto c) {
) {
CHKERR mField.get_moab().get_entities_by_dimension(
skin = subtract(skin, ents);
}
};
auto remove_named_blocks = [&](
auto n) {
std::regex(
(boost::format(
"%s(.*)") %
n).str()
))
) {
CHKERR mField.get_moab().get_entities_by_dimension(
skin = subtract(skin, ents);
}
};
"remove_cubit_blocks");
"remove_cubit_blocks");
"remove_named_blocks");
return skin;
};
auto filter_true_skin = [&](auto skin) {
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
auto remove_flux_ents = filter_true_skin(filter_flux_blocks(get_skin()));
remove_flux_ents);
MOFEM_LOG(
"SYNC", Sev::noisy) << remove_flux_ents << endl;
#ifndef NDEBUG
mField.get_moab(),
(boost::format("flux_remove_%d.vtk") % mField.get_comm_rank()).str(),
remove_flux_ents);
#endif
simple->getProblemName(),
"FLUX", remove_flux_ents);
auto set_init_temp = [](boost::shared_ptr<FieldEntity> field_entity_ptr) {
field_entity_ptr->getEntFieldData()[0] =
init_temp;
return 0;
};
CHKERR mField.getInterface<
FieldBlas>()->fieldLambdaOnEntities(set_init_temp,
"T");
simple->getProblemName(),
"U");
simple->getProblemName(),
"FLUX",
false);
}
auto boundary_marker =
bc_mng->getMergedBlocksMarker(vector<string>{"HEATFLUX"});
};
auto block_params = boost::make_shared<BlockedParameters>();
auto mDPtr = block_params->getDPtr();
auto coeff_expansion_ptr = block_params->getCoeffExpansionPtr();
auto heat_conductivity_ptr = block_params->getHeatConductivityPtr();
auto heat_capacity_ptr = block_params->getHeatCapacityPtr();
auto time_scale = boost::make_shared<TimeScale>();
auto time_bodyforce_scaling =
boost::make_shared<TimeScale>("bodyforce_scale.txt");
auto time_heatsource_scaling =
boost::make_shared<TimeScale>("heatsource_scale.txt");
auto time_temperature_scaling =
boost::make_shared<TimeScale>("temperature_bc_scale.txt");
auto time_displacement_scaling =
boost::make_shared<TimeScale>("displacement_bc_scale.txt");
auto time_flux_scaling = boost::make_shared<TimeScale>("flux_bc_scale.txt");
auto time_force_scaling =
boost::make_shared<TimeScale>("force_bc_scale.txt");
auto add_domain_rhs_ops = [&](auto &pipeline) {
Sev::inform);
auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
auto vec_temp_ptr = boost::make_shared<VectorDouble>();
auto vec_temp_dot_ptr = boost::make_shared<VectorDouble>();
auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
auto vec_temp_div_ptr = boost::make_shared<VectorDouble>();
pipeline.push_back(
"FLUX", vec_temp_div_ptr));
pipeline.push_back(
"U", mat_grad_ptr));
pipeline.push_back(
pipeline.push_back(
new OpStressThermal(mat_strain_ptr, vec_temp_ptr, mDPtr,
coeff_expansion_ptr,
mat_stress_ptr));
pipeline.push_back(
new OpSetBc(
"FLUX",
true, boundary_marker));
"U", mat_stress_ptr,
[](double, double, double) constexpr { return 1; }));
auto resistance = [heat_conductivity_ptr](
const double,
const double,
return (1. / (*heat_conductivity_ptr));
};
auto capacity = [heat_capacity_ptr](
const double,
const double,
return -(*heat_capacity_ptr);
};
return -1.;
};
pipeline.push_back(
new OpHdivFlux(
"FLUX", mat_flux_ptr, resistance));
pipeline.push_back(
new OpHDivTemp(
"FLUX", vec_temp_ptr, unity));
pipeline.push_back(
new OpBaseDivFlux(
"T", vec_temp_div_ptr, unity));
pipeline.push_back(
new OpBaseDotT(
"T", vec_temp_dot_ptr, capacity));
pipeline, mField, "T", {time_scale, time_heatsource_scaling},
"HEAT_SOURCE", Sev::inform);
pipeline, mField, "U", {time_scale, time_bodyforce_scaling},
"BODY_FORCE", Sev::inform);
pipeline, mField, "T", vec_temp_ptr, "SETTEMP", Sev::inform);
};
auto add_domain_lhs_ops = [&](auto &pipeline) {
Sev::verbose);
pipeline.push_back(
new OpSetBc(
"FLUX",
true, boundary_marker));
pipeline.push_back(
new OpKCauchy(
"U",
"U", mDPtr));
"U", "T", mDPtr, coeff_expansion_ptr));
auto resistance = [heat_conductivity_ptr](
const double,
const double,
return (1. / (*heat_conductivity_ptr));
};
auto capacity = [heat_capacity_ptr](
const double,
const double,
return -(*heat_capacity_ptr);
};
pipeline.push_back(
new OpHdivHdiv(
"FLUX",
"FLUX", resistance));
"FLUX", "T", []() constexpr { return -1; }, true));
auto op_capacity =
new OpCapacity(
"T",
"T", capacity);
op_capacity->feScalingFun = [](
const FEMethod *fe_ptr) {
return fe_ptr->ts_a;
};
pipeline.push_back(op_capacity);
auto vec_temp_ptr = boost::make_shared<VectorDouble>();
pipeline, mField, "T", vec_temp_ptr, "SETTEMP", Sev::verbose);
};
auto add_boundary_rhs_ops = [&](auto &pipeline) {
pipeline.push_back(
new OpSetBc(
"FLUX",
true, boundary_marker));
pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U",
{time_scale, time_force_scaling}, "FORCE", Sev::inform);
pipeline_mng->getOpBoundaryRhsPipeline(), mField, "FLUX",
{time_scale, time_temperature_scaling}, "TEMPERATURE", Sev::inform);
auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
mField, pipeline,
simple->getProblemName(),
"FLUX", mat_flux_ptr,
{time_scale, time_flux_scaling});
};
auto add_boundary_lhs_ops = [&](auto &pipeline) {
mField, pipeline,
simple->getProblemName(),
"FLUX");
};
auto get_bc_hook_rhs = [&]() {
mField, pipeline_mng->getDomainRhsFE(),
{time_scale, time_displacement_scaling});
return hook;
};
auto get_bc_hook_lhs = [&]() {
mField, pipeline_mng->getDomainLhsFE(),
{time_scale, time_displacement_scaling});
return hook;
};
pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
pipeline_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
CHKERR add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
CHKERR add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
CHKERR add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
CHKERR add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
}
auto solver = pipeline_mng->createTSIM();
auto set_section_monitor = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
(void *)(snes_ctx_ptr.get()), nullptr);
};
auto create_post_process_element = [&]() {
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
auto block_params = boost::make_shared<BlockedParameters>();
auto mDPtr = block_params->getDPtr();
auto coeff_expansion_ptr = block_params->getCoeffExpansionPtr();
"MAT_THERMAL", block_params, Sev::verbose);
post_proc_fe->getOpPtrVector(), {H1, HDIV});
auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
auto vec_temp_ptr = boost::make_shared<VectorDouble>();
auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
auto u_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
mat_grad_ptr));
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
coeff_expansion_ptr, mat_stress_ptr));
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"T", vec_temp_ptr}},
{{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
{},
{{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
)
);
return post_proc_fe;
};
auto monitor_ptr = boost::make_shared<FEMethod>();
auto post_proc_fe = create_post_process_element();
auto set_time_monitor = [&](auto dm, auto solver) {
monitor_ptr->preProcessHook = [&]() {
post_proc_fe,
monitor_ptr->getCacheWeakPtr());
CHKERR post_proc_fe->writeFile(
"out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
".h5m");
->evalFEAtThePoint3D(
fieldEvalCoords.data(), 1e-12,
simple->getProblemName(),
simple->getDomainFEName(), fieldEvalData,
mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
} else {
->evalFEAtThePoint2D(
fieldEvalCoords.data(), 1e-12,
simple->getProblemName(),
simple->getDomainFEName(), fieldEvalData,
mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
}
if (scalarFieldPtr->size()) {
<< "Eval point T: " << t_temp;
}
if (vectorFieldPtr->size1()) {
auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
<<
"Eval point U magnitude: " << sqrt(t_disp(
i) * t_disp(
i));
}
if (tensorFieldPtr->size1()) {
auto t_disp_grad =
getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*tensorFieldPtr);
<<
"Eval point U_GRAD trace: " << t_disp_grad(
i,
i);
}
}
};
auto null = boost::shared_ptr<FEMethod>();
monitor_ptr, null);
};
auto set_fieldsplit_preconditioner = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
auto name_prb =
simple->getProblemName();
CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"U", 0,
CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"FLUX", 0, 0,
is_flux);
CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"T", 0, 0,
is_T);
IS is_tmp;
CHKERR ISExpand(is_T, is_flux, &is_tmp);
auto is_all_bc = bc_mng->getBlockIS(name_prb, "HEATFLUX", "FLUX", 0, 0);
int is_all_bc_size;
CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
<< "Field split block size " << is_all_bc_size;
if (is_all_bc_size) {
IS is_tmp2;
CHKERR ISDifference(is_TFlux, is_all_bc, &is_tmp2);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
is_all_bc);
}
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_TFlux);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
}
};
CHKERR TSSetFromOptions(solver);
CHKERR set_section_monitor(solver);
CHKERR set_fieldsplit_preconditioner(solver);
CHKERR set_time_monitor(dm, solver);
}
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
auto core_log = logging::core::get();
core_log->add_sink(
core_log->add_sink(
try {
DMType dm_name = "DMMOFEM";
}
}