#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 2
#endif
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
GAUSS>::OpMixDivTimesScalar<SPACE_DIM>;
GAUSS>::OpBaseTimesVector<3, 3, 1>;
GAUSS>::OpMixDivTimesU<3, 1, 2>;
GAUSS>::OpBaseTimesScalarField<1>;
GAUSS>::OpEssentialRhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
GAUSS>::OpEssentialLhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
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;
struct BlockedParameters
: public boost::enable_shared_from_this<BlockedParameters> {
double fluidConductivity;
inline auto getDPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &mD);
}
inline auto getConductivityPtr() {
return boost::shared_ptr<double>(shared_from_this(), &fluidConductivity);
}
};
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()
))
));
OpMatFluidBlocks(boost::shared_ptr<double> conductivity_ptr,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
conductivityPtr(conductivity_ptr) {
"Can not get data from block");
}
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
*conductivityPtr = b.conductivity;
}
}
}
private:
double conductivity;
};
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() < 1) {
"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() <<
": conductivity = " << block_data[0];
blockData.push_back({block_data[0], get_block_ents()});
}
}
boost::shared_ptr<double> expansionPtr;
boost::shared_ptr<double> conductivityPtr;
boost::shared_ptr<double> capacityPtr;
};
pipeline.push_back(new OpMatFluidBlocks(
blockedParamsPtr->getConductivityPtr(), mField, sev,
(boost::format("%s(.*)") % block_thermal_name).str()
))
));
}
}
}
auto get_command_line_parameters = [&]() {
fieldEvalCoords.data(), &coords_dim,
};
CHKERR get_command_line_parameters();
}
simple->getProblemName(),
"U");
simple->getProblemName(),
"FLUX",
false);
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");
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 boundary_marker =
bc_mng->getMergedBlocksMarker(vector<string>{"FLUIDFLUX"});
auto u_grad_ptr = boost::make_shared<MatrixDouble>();
auto dot_u_grad_ptr = boost::make_shared<MatrixDouble>();
auto trace_dot_u_grad_ptr = boost::make_shared<VectorDouble>();
auto h_ptr = boost::make_shared<VectorDouble>();
auto dot_h_ptr = boost::make_shared<VectorDouble>();
auto flux_ptr = boost::make_shared<MatrixDouble>();
auto div_flux_ptr = boost::make_shared<VectorDouble>();
auto strain_ptr = boost::make_shared<MatrixDouble>();
auto stress_ptr = boost::make_shared<MatrixDouble>();
auto time_scale = boost::make_shared<TimeScale>();
auto block_params = boost::make_shared<BlockedParameters>();
auto mDPtr = block_params->getDPtr();
auto conductivity_ptr = block_params->getConductivityPtr();
};
auto add_domain_base_ops = [&](auto &pip) {
Sev::inform);
"U", u_grad_ptr));
"U", dot_u_grad_ptr));
trace_dot_u_grad_ptr));
"FLUX", div_flux_ptr));
};
auto add_domain_ops_lhs_mechanical = [&](auto &pip) {
pip.push_back(
new OpKCauchy(
"U",
"U", mDPtr));
"H", "U",
[](const double, const double, const double) { return -9.81; }, true,
true));
};
auto add_domain_ops_rhs_mechanical = [&](auto &pip) {
pip, mField, "U", {time_scale}, "BODY_FORCE", Sev::inform);
"U", strain_ptr, stress_ptr, mDPtr));
pip.push_back(
};
auto add_domain_ops_lhs_seepage = [&](auto &pip, auto &fe) {
auto resistance = [conductivity_ptr](
const double,
const double,
return (1. / *(conductivity_ptr));
};
auto unity = []() constexpr { return -1; };
pip.push_back(
new OpHdivHdiv(
"FLUX",
"FLUX", resistance));
pip.push_back(
new OpHdivQ(
"FLUX",
"H", unity,
true));
"H", "U", [](double, double, double) constexpr { return -1; }, false,
false);
op_base_div_u->feScalingFun = [](
const FEMethod *fe_ptr) {
return fe_ptr->ts_a;
};
pip.push_back(op_base_div_u);
};
auto add_domain_ops_rhs_seepage = [&](auto &pip) {
return (1. / (*conductivity_ptr));
};
pip.push_back(
new OpHdivFlux(
"FLUX", flux_ptr, resistance));
pip.push_back(
new OpHDivH(
"FLUX", h_ptr, minus_one));
pip.push_back(
new OpBaseDotH(
"H", trace_dot_u_grad_ptr, minus_one));
};
auto add_boundary_rhs_ops = [&](auto &pip) {
pip.push_back(
new OpSetBc(
"FLUX",
true, boundary_marker));
pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
"FORCE", Sev::inform);
pipeline_mng->getOpBoundaryRhsPipeline(), mField, "FLUX", {time_scale},
"PRESSURE", Sev::inform);
auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
mField, pip,
simple->getProblemName(),
"FLUX", mat_flux_ptr,
{time_scale});
};
auto add_boundary_lhs_ops = [&](auto &pip) {
mField, pip,
simple->getProblemName(),
"FLUX");
};
CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
CHKERR add_domain_ops_lhs_mechanical(pipeline_mng->getOpDomainLhsPipeline());
CHKERR add_domain_ops_lhs_seepage(pipeline_mng->getOpDomainLhsPipeline(),
pipeline_mng->getDomainLhsFE());
CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
CHKERR add_domain_ops_rhs_mechanical(pipeline_mng->getOpDomainRhsPipeline());
CHKERR add_domain_ops_rhs_seepage(pipeline_mng->getOpDomainRhsPipeline());
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();
"MAT_FLUID", 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 h_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(
"U", mat_strain_ptr, mat_stress_ptr, mDPtr));
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"H", h_ptr}},
{{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
{},
{{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
)
);
return post_proc_fe;
};
auto create_creaction_fe = [&]() {
auto fe_ptr = boost::make_shared<DomainEle>(mField);
fe_ptr->getRuleHook = [](
int,
int,
int o) {
return 2 * o; };
auto &pip = fe_ptr->getOpPtrVector();
auto block_params = boost::make_shared<BlockedParameters>();
auto mDPtr = block_params->getDPtr();
Sev::verbose);
auto u_grad_ptr = boost::make_shared<MatrixDouble>();
auto strain_ptr = boost::make_shared<MatrixDouble>();
auto stress_ptr = boost::make_shared<MatrixDouble>();
"U", u_grad_ptr));
"U", strain_ptr, stress_ptr, mDPtr));
fe_ptr->postProcessHook =
return fe_ptr;
};
auto monitor_ptr = boost::make_shared<FEMethod>();
auto post_proc_fe = create_post_process_element();
auto rections_fe = create_creaction_fe();
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");
rections_fe->f = res;
rections_fe,
monitor_ptr->getCacheWeakPtr());
double nrm;
CHKERR VecNorm(res, NORM_2, &nrm);
<< "Residual norm " << nrm << " at time step "
<< monitor_ptr->ts_step;
auto scalar_field_ptr = boost::make_shared<VectorDouble>();
auto vector_field_ptr = boost::make_shared<MatrixDouble>();
auto tensor_field_ptr = boost::make_shared<MatrixDouble>();
->getData<DomainEle>();
field_eval_data,
simple->getDomainFEName());
} else {
field_eval_data,
simple->getDomainFEName());
}
field_eval_data->setEvalPoints(fieldEvalCoords.data(), 1);
auto no_rule = [](
int,
int,
int) {
return -1; };
auto field_eval_ptr = field_eval_data->feMethodPtr.lock();
field_eval_ptr->getRuleHook = no_rule;
auto u_grad_ptr = boost::make_shared<MatrixDouble>();
auto strain_ptr = boost::make_shared<MatrixDouble>();
auto stress_ptr = boost::make_shared<MatrixDouble>();
auto h_ptr = boost::make_shared<VectorDouble>();
auto block_params = boost::make_shared<BlockedParameters>();
auto mDPtr = block_params->getDPtr();
"MAT_FLUID", block_params, Sev::noisy);
field_eval_ptr->getOpPtrVector(), {H1, HDIV});
field_eval_ptr->getOpPtrVector().push_back(
"U", u_grad_ptr));
field_eval_ptr->getOpPtrVector().push_back(
field_eval_ptr->getOpPtrVector().push_back(
field_eval_ptr->getOpPtrVector().push_back(
"U", strain_ptr, stress_ptr, mDPtr));
->evalFEAtThePoint3D(
fieldEvalCoords.data(), 1e-12,
simple->getProblemName(),
simple->getDomainFEName(), field_eval_data,
mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
} else {
->evalFEAtThePoint2D(
fieldEvalCoords.data(), 1e-12,
simple->getProblemName(),
simple->getDomainFEName(), field_eval_data,
mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
}
<< "Eval point hydrostatic hight: " << *h_ptr;
<< "Eval point skeleton stress pressure: " << *stress_ptr;
}
};
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,
"H", 0, 0,
is_H);
IS is_tmp;
CHKERR ISExpand(is_H, is_flux, &is_tmp);
auto is_all_bc = bc_mng->getBlockIS(name_prb, "FLUIDFLUX", "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_Flux, is_all_bc, &is_tmp2);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
is_all_bc);
}
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_Flux);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
}
};
auto pre_proc_ptr = boost::make_shared<FEMethod>();
auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
auto time_scale = boost::make_shared<TimeScale>();
auto get_bc_hook_rhs = [this, pre_proc_ptr, time_scale]() {
{time_scale}, false);
return hook;
};
auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
mField, post_proc_rhs_ptr, 1.)();
};
auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
post_proc_lhs_ptr, 1.);
};
pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
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(
LogManager::createSink(LogManager::getStrmWorld(), "Seepage"));
LogManager::setLog("Seepage");
core_log->add_sink(
LogManager::createSink(LogManager::getStrmSync(), "SeepageSync"));
LogManager::setLog("SeepageSync");
try {
DMType dm_name = "DMMOFEM";
}
}