9#ifndef EXECUTABLE_DIMENSION
10#define EXECUTABLE_DIMENSION 3
66 GAUSS>::OpMixDivTimesScalar<SPACE_DIM>;
80 GAUSS>::OpBaseTimesVector<3, SPACE_DIM, 1>;
87 GAUSS>::OpMixDivTimesU<3, 1, 2>;
95 GAUSS>::OpBaseTimesScalar<1>;
205 auto get_command_line_parameters = [&]() {
228 <<
"Reference_temperature " <<
ref_temp;
233 CHKERR get_command_line_parameters();
247 simple->getProblemName(),
"U");
249 simple->getProblemName(),
"FLUX");
261 auto boundary_marker =
262 bc_mng->getMergedBlocksMarker(vector<string>{
"HEATFLUX"});
273 auto time_scale = boost::make_shared<TimeScale>();
275 auto add_domain_ops = [&](
auto &pipeline) {
278 auto det_ptr = boost::make_shared<VectorDouble>();
279 auto jac_ptr = boost::make_shared<MatrixDouble>();
280 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
303 auto add_domain_rhs_ops = [&](
auto &pipeline) {
306 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
307 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
308 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
310 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
311 auto vec_temp_dot_ptr = boost::make_shared<VectorDouble>();
312 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
313 auto vec_temp_div_ptr = boost::make_shared<VectorDouble>();
319 "FLUX", vec_temp_div_ptr));
327 pipeline.push_back(
new OpStressThermal(
"U", mat_strain_ptr, vec_temp_ptr,
328 mDPtr, mat_stress_ptr));
330 pipeline.push_back(
new OpSetBc(
"FLUX",
true, boundary_marker));
334 [](
double,
double,
double)
constexpr {
return 1; }));
343 pipeline.push_back(
new OpHdivFlux(
"FLUX", mat_flux_ptr, resistance));
344 pipeline.push_back(
new OpHDivTemp(
"FLUX", vec_temp_ptr, unity));
345 pipeline.push_back(
new OpBaseDivFlux(
"T", vec_temp_div_ptr, unity));
346 pipeline.push_back(
new OpBaseDotT(
"T", vec_temp_dot_ptr, capacity));
349 pipeline,
mField,
"T", {time_scale},
"HEAT_SOURCE", Sev::inform);
351 pipeline_mng->getOpDomainRhsPipeline(),
mField,
"U", {time_scale},
352 "BODY_FORCE", Sev::inform);
354 pipeline.push_back(
new OpUnSetBc(
"FLUX"));
358 auto add_domain_lhs_ops = [&](
auto &pipeline) {
360 pipeline.push_back(
new OpSetBc(
"FLUX",
true, boundary_marker));
362 pipeline.push_back(
new OpKCauchy(
"U",
"U", mDPtr));
372 pipeline.push_back(
new OpHdivHdiv(
"FLUX",
"FLUX", resistance));
373 pipeline.push_back(
new OpHdivT(
374 "FLUX",
"T", []() {
return -1; },
true));
376 auto op_capacity =
new OpCapacity(
"T",
"T", capacity);
377 op_capacity->feScalingFun = [](
const FEMethod *fe_ptr) {
380 pipeline.push_back(op_capacity);
382 pipeline.push_back(
new OpUnSetBc(
"FLUX"));
386 auto add_boundary_rhs_ops = [&](
auto &pipeline) {
395 pipeline.push_back(
new OpSetBc(
"FLUX",
true, boundary_marker));
398 pipeline_mng->getOpBoundaryRhsPipeline(),
mField,
"U", {time_scale},
399 "FORCE", Sev::inform);
402 pipeline_mng->getOpBoundaryRhsPipeline(),
mField,
"FLUX", {time_scale},
403 "TEMPERATURE", Sev::inform);
405 pipeline.push_back(
new OpUnSetBc(
"FLUX"));
407 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
412 mField, pipeline,
simple->getProblemName(),
"FLUX", mat_flux_ptr,
418 auto add_boundary_lhs_ops = [&](
auto &pipeline) {
434 auto get_bc_hook_rhs = [&]() {
436 mField, pipeline_mng->getDomainRhsFE(), {time_scale});
439 auto get_bc_hook_lhs = [&]() {
441 mField, pipeline_mng->getDomainLhsFE(), {time_scale});
445 pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
446 pipeline_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
448 CHKERR add_domain_ops(pipeline_mng->getOpDomainRhsPipeline());
449 CHKERR add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
451 CHKERR add_domain_ops(pipeline_mng->getOpDomainLhsPipeline());
452 CHKERR add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
454 CHKERR add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
455 CHKERR add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
469 auto dm =
simple->getDM();
470 auto solver = pipeline_mng->createTSIM();
473 auto set_section_monitor = [&](
auto solver) {
476 CHKERR TSGetSNES(solver, &snes);
477 CHKERR SNESMonitorSet(snes,
480 (
void *)(snes_ctx_ptr.get()),
nullptr);
484 auto create_post_process_element = [&]() {
485 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
487 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
488 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
489 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
491 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
492 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
494 auto det_ptr = boost::make_shared<VectorDouble>();
495 auto jac_ptr = boost::make_shared<MatrixDouble>();
496 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
498 post_proc_fe->getOpPtrVector().push_back(
500 post_proc_fe->getOpPtrVector().push_back(
502 post_proc_fe->getOpPtrVector().push_back(
507 post_proc_fe->getOpPtrVector().push_back(
509 post_proc_fe->getOpPtrVector().push_back(
511 post_proc_fe->getOpPtrVector().push_back(
514 post_proc_fe->getOpPtrVector().push_back(
516 post_proc_fe->getOpPtrVector().push_back(
520 post_proc_fe->getOpPtrVector().push_back(
522 post_proc_fe->getOpPtrVector().push_back(
525 auto u_ptr = boost::make_shared<MatrixDouble>();
526 post_proc_fe->getOpPtrVector().push_back(
528 post_proc_fe->getOpPtrVector().push_back(
531 post_proc_fe->getOpPtrVector().push_back(
534 "U", mat_strain_ptr, vec_temp_ptr,
getMatDPtr(), mat_stress_ptr));
538 post_proc_fe->getOpPtrVector().push_back(
542 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
544 {{
"T", vec_temp_ptr}},
546 {{
"U", u_ptr}, {
"FLUX", mat_flux_ptr}},
550 {{
"STRAIN", mat_strain_ptr}, {
"STRESS", mat_stress_ptr}}
559 auto monitor_ptr = boost::make_shared<FEMethod>();
560 auto post_proc_fe = create_post_process_element();
562 auto set_time_monitor = [&](
auto dm,
auto solver) {
564 monitor_ptr->preProcessHook = [&]() {
567 CHKERR post_proc_fe->writeFile(
568 "out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
572 auto null = boost::shared_ptr<FEMethod>();
578 auto set_fieldsplit_preconditioner = [&](
auto solver) {
582 CHKERR TSGetSNES(solver, &snes);
584 CHKERR SNESGetKSP(snes, &ksp);
586 CHKERR KSPGetPC(ksp, &pc);
587 PetscBool is_pcfs = PETSC_FALSE;
588 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
591 if (is_pcfs == PETSC_TRUE) {
592 auto bc_mng = mField.getInterface<
BcManager>();
593 auto name_prb =
simple->getProblemName();
594 auto is_all_bc = bc_mng->getBlockIS(name_prb,
"HEATFLUX",
"FLUX", 0, 0);
596 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
598 <<
"Field split block size " << is_all_bc_size;
599 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
607 CHKERR TSSetSolution(solver,
D);
608 CHKERR TSSetFromOptions(solver);
609 CHKERR set_section_monitor(solver);
610 CHKERR set_time_monitor(dm, solver);
611 CHKERR set_fieldsplit_preconditioner(solver);
613 CHKERR TSSolve(solver, NULL);
620 auto set_matrial_stiffness = [&](
auto mDPtr) {
634 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mDPtr);
642 return set_matrial_stiffness(
648int main(
int argc,
char *argv[]) {
655 auto core_log = logging::core::get();
664 DMType dm_name =
"DMMOFEM";
669 moab::Core mb_instance;
670 moab::Interface &moab = mb_instance;
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#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 ...
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
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.
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.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Only used with Hooke equation (linear material model)]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Simple interface for fast problem set-up.
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)
Essential boundary conditions.
Class (Function) to enforce essential constrains.
structure for User Loop Methods on finite elements
default operator for TRI element
Definition of the heat flux bc data structure.
Section manager is used to create indexes and sections.
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.
Get vector field for H-div approximation.
Calculate divergence of vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Enforce essential constrains on lhs.
Enforce essential constrains on rhs.
Make Hdiv space from Hcurl space in 2d.
Post post-proc data at points from hash maps.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
transform local reference derivatives of shape function to global derivatives if higher order geometr...
PipelineManager interface.
MoFEM::VolumeElementForcesAndSourcesCore VolEle
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
MoFEMErrorCode runProblem()
[Run problem]
MoFEM::Interface & mField
MoFEMErrorCode createCommonData()
[Set up problem]
ThermoElasticProblem(MoFEM::Interface &m_field)
boost::shared_ptr< MatrixDouble > getMatDPtr()
[Solve]
MoFEMErrorCode setupProblem()
add fields
MoFEMErrorCode bC()
[Create common data]
MoFEMErrorCode OPs()
[Boundary condition]
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
double young_modulus
[Essential boundary conditions (Least square approach)]
#define EXECUTABLE_DIMENSION
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBaseDotT
Integrate Rhs base of temperature time heat capacity times heat rate (T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< SPACE_DIM > OpHdivT
Integrate Lhs div of base of flux time base of temperature (FLUX x T) and transpose of it,...
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, SPACE_DIM, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpCapacity
Integrate Lhs base of temperature times (heat capacity) times base of temperature (T x T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, SPACE_DIM > OpHdivHdiv
[Linear elastic problem]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, 2 > OpHDivTemp
Integrate Rhs div flux base times temperature (T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Linear elastic problem]
OpBaseDotT OpBaseDivFlux
Integrate Rhs base of temperature times divergence of flux (T)