26 using namespace MoFEM;
28 static char help[] =
"...\n\n";
30 #include <BasicFiniteElements.hpp>
94 constexpr
double a0 = 0.98;
96 constexpr
double mu_m = 0.0101;
97 constexpr
double rho_p = 0.0012;
98 constexpr
double mu_p = 0.000182;
100 constexpr
double W = 0.25;
103 constexpr
double h = 0.02;
108 constexpr
double md = 1e-2;
109 constexpr
double eps = 1e-12;
110 constexpr
double tol = std::numeric_limits<float>::epsilon();
131 auto my_max = [](
const double x) {
return (x - 1 + std::abs(x + 1)) / 2; };
132 auto my_min = [](
const double x) {
return (x + 1 - std::abs(x - 1)) / 2; };
137 if (
h >= -1 &&
h < 1)
151 auto get_f = [](
const double h) {
return 4 *
W *
h * (
h *
h - 1); };
152 auto get_f_dh = [](
const double h) {
return 4 *
W * (3 *
h *
h - 1); };
159 return md * (1 -
h *
h);
169 const double h2 =
h *
h;
170 const double h3 = h2 *
h;
172 return md * (2 * h3 - 3 * h2 + 1);
174 return md * (-2 * h3 - 3 * h2 + 1);
196 constexpr
double R = 0.0125;
197 constexpr
double A = R * 0.2;
198 const double theta = atan2(
r, y);
199 const double w = R +
A * cos(
n * theta);
200 const double d = std::sqrt(
r *
r + y * y);
201 return tanh((
w -
d) / (
eta * std::sqrt(2)));
205 constexpr
double y0 = 0.5;
206 constexpr
double R = 0.4;
208 const double d = std::sqrt(
r *
r + y * y);
209 return tanh((R -
d) / (
eta * std::sqrt(2)));
212 auto init_h = [](
double r,
double y,
double theta) {
242 CHKERR boundaryCondition();
284 constexpr
int order = 3;
303 auto dm =
simple->getDM();
305 auto h_ptr = boost::make_shared<VectorDouble>();
306 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
307 auto g_ptr = boost::make_shared<VectorDouble>();
308 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
310 auto post_proc = [&]() {
312 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
313 post_proc_fe->generateReferenceElementMesh();
315 auto det_ptr = boost::make_shared<VectorDouble>();
316 auto jac_ptr = boost::make_shared<MatrixDouble>();
317 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
319 post_proc_fe->addFieldValuesPostProc(
"H");
320 post_proc_fe->addFieldValuesPostProc(
"G");
321 post_proc_fe->addFieldValuesGradientPostProc(
"G", 2);
322 post_proc_fe->addFieldValuesGradientPostProc(
"H", 2);
325 CHKERR post_proc_fe->writeFile(
"out_init.h5m");
330 auto solve_init = [&]() {
333 auto set_generic = [&](
auto &pipeline) {
334 auto det_ptr = boost::make_shared<VectorDouble>();
335 auto jac_ptr = boost::make_shared<MatrixDouble>();
336 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
352 auto set_domain_rhs = [&](
auto &pipeline) {
353 pipeline.push_back(
new OpRhsH<true>(
"H",
nullptr,
nullptr, h_ptr,
354 grad_h_ptr, grad_g_ptr));
355 pipeline.push_back(
new OpRhsG<true>(
"G", h_ptr, grad_h_ptr, g_ptr));
358 auto set_domain_lhs = [&](
auto &pipeline) {
359 pipeline.push_back(
new OpLhsH_dH<true>(
"H",
nullptr, h_ptr, grad_g_ptr));
365 auto create_subdm = [&]() {
367 CHKERR DMCreate(mField.get_comm(), &subdm);
368 CHKERR DMSetType(subdm,
"DMMOFEM");
380 auto subdm = create_subdm();
384 set_generic(pipeline_mng->getOpDomainRhsPipeline());
385 set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
386 set_generic(pipeline_mng->getOpDomainLhsPipeline());
387 set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
390 auto snes = pipeline_mng->createSNES(subdm);
392 auto set_section_monitor = [&](
auto solver) {
394 PetscViewerAndFormat *vf;
395 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
396 PETSC_VIEWER_DEFAULT, &vf);
400 void *))SNESMonitorFields,
402 auto section = mField.getInterface<
ISManager>()->sectionCreate(
"SUB");
404 CHKERR PetscSectionGetNumFields(section, &num_fields);
405 for (
int f = 0;
f < num_fields; ++
f) {
406 const char *field_name;
407 CHKERR PetscSectionGetFieldName(section,
f, &field_name);
409 <<
"Field " <<
f <<
" " << std::string(field_name);
414 CHKERR set_section_monitor(snes);
416 CHKERR SNESSetFromOptions(snes);
417 CHKERR SNESSolve(snes, PETSC_NULL,
D);
419 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
420 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
429 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
430 "SYMETRY",
"U", 0, 0);
431 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
432 "SYMETRY",
"L", 0, 0);
433 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
435 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
438 pipeline_mng->getOpDomainRhsPipeline().clear();
439 pipeline_mng->getOpDomainLhsPipeline().clear();
449 auto dot_u_ptr = boost::make_shared<MatrixDouble>();
450 auto u_ptr = boost::make_shared<MatrixDouble>();
451 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
452 auto dot_h_ptr = boost::make_shared<VectorDouble>();
453 auto h_ptr = boost::make_shared<VectorDouble>();
454 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
455 auto g_ptr = boost::make_shared<VectorDouble>();
456 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
457 auto lambda_ptr = boost::make_shared<VectorDouble>();
458 auto p_ptr = boost::make_shared<VectorDouble>();
459 auto div_u_ptr = boost::make_shared<VectorDouble>();
463 auto set_domain_general = [&](
auto &pipeline) {
464 auto det_ptr = boost::make_shared<VectorDouble>();
465 auto jac_ptr = boost::make_shared<MatrixDouble>();
466 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
494 auto set_domain_rhs = [&](
auto &pipeline) {
495 pipeline.push_back(
new OpRhsU(
"U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
496 grad_h_ptr, g_ptr, p_ptr));
497 pipeline.push_back(
new OpRhsH<false>(
"H", u_ptr, dot_h_ptr, h_ptr,
498 grad_h_ptr, grad_g_ptr));
499 pipeline.push_back(
new OpRhsG<false>(
"G", h_ptr, grad_h_ptr, g_ptr));
501 "P", div_u_ptr, [](
const double r,
const double,
const double) {
505 "P", p_ptr, [](
const double r,
const double,
const double) {
510 auto set_domain_lhs = [&](
auto &pipeline) {
511 pipeline.push_back(
new OpLhsU_dU(
"U", u_ptr, grad_u_ptr, h_ptr));
513 new OpLhsU_dH(
"U",
"H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
514 pipeline.push_back(
new OpLhsU_dG(
"U",
"G", grad_h_ptr));
516 pipeline.push_back(
new OpLhsH_dU(
"H",
"U", grad_h_ptr));
525 [](
const double r,
const double,
const double) {
535 auto set_boundary_rhs = [&](
auto &pipeline) {
543 auto set_boundary_lhs = [&](
auto &pipeline) {
554 set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
555 set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
556 set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
557 set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
558 set_boundary_rhs(pipeline_mng->getOpBoundaryRhsPipeline());
559 set_boundary_lhs(pipeline_mng->getOpBoundaryLhsPipeline());
561 domianLhsFEPtr = pipeline_mng->getDomainLhsFE();
576 std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
578 : dM(dm), postProc(post_proc), liftFE(
p.first), liftVec(
p.second) {}
584 this->getCacheWeakPtr());
585 CHKERR postProc->writeFile(
586 "out_step_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
601 MPI_Allreduce(MPI_IN_PLACE, &(*liftVec)[0],
SPACE_DIM, MPI_DOUBLE, MPI_SUM,
604 <<
"Step " << ts_step <<
" time " << ts_t
605 <<
" lift vec x: " << (*liftVec)[0] <<
" y: " << (*liftVec)[1];
612 boost::shared_ptr<PostProcEle> postProc;
623 auto dm =
simple->getDM();
625 auto get_fe_post_proc = [&]() {
626 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
627 post_proc_fe->generateReferenceElementMesh();
628 auto det_ptr = boost::make_shared<VectorDouble>();
629 auto jac_ptr = boost::make_shared<MatrixDouble>();
630 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
632 post_proc_fe->getOpPtrVector().push_back(
634 post_proc_fe->getOpPtrVector().push_back(
636 post_proc_fe->getOpPtrVector().push_back(
639 post_proc_fe->addFieldValuesPostProc(
"U");
640 post_proc_fe->addFieldValuesPostProc(
"H");
641 post_proc_fe->addFieldValuesPostProc(
"P");
642 post_proc_fe->addFieldValuesPostProc(
"G");
643 post_proc_fe->addFieldValuesGradientPostProc(
"U", 2);
644 post_proc_fe->addFieldValuesGradientPostProc(
"H", 2);
645 post_proc_fe->addFieldValuesGradientPostProc(
"G", 2);
649 auto get_lift_fe = [&]() {
650 auto fe = boost::make_shared<BoundaryEle>(mField);
651 auto lift_ptr = boost::make_shared<VectorDouble>();
652 auto p_ptr = boost::make_shared<VectorDouble>();
653 auto ents_ptr = boost::make_shared<Range>();
655 std::vector<const CubitMeshSets *> vec_ptr;
657 std::regex(
"LIFT"), vec_ptr);
658 for (
auto m_ptr : vec_ptr) {
661 CHKERR mField.get_moab().get_entities_by_dimension(meshset,
SPACE_DIM - 1,
663 ents_ptr->merge(ents);
666 MOFEM_LOG(
"FS", Sev::inform) <<
"Lift ents " << (*ents_ptr);
668 fe->getOpPtrVector().push_back(
670 fe->getOpPtrVector().push_back(
673 return std::make_pair(fe, lift_ptr);
676 auto set_ts = [&](
auto solver) {
679 CHKERR TSGetSNES(solver, &snes);
681 CHKERR SNESGetKSP(snes, &ksp);
685 auto ts = pipeline_mng->createTSIM();
686 CHKERR TSSetType(ts, TSALPHA);
688 auto set_post_proc_monitor = [&](
auto dm) {
690 boost::shared_ptr<FEMethod> null_fe;
692 boost::make_shared<Monitor>(dm, get_fe_post_proc(), get_lift_fe());
694 null_fe, monitor_ptr);
697 CHKERR set_post_proc_monitor(dm);
699 auto set_section_monitor = [&](
auto solver) {
702 CHKERR TSGetSNES(solver, &snes);
703 PetscViewerAndFormat *vf;
704 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
705 PETSC_VIEWER_DEFAULT, &vf);
708 (
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
711 auto section = mField.getInterface<
ISManager>()->sectionCreate(
712 simple->getProblemName());
714 CHKERR PetscSectionGetNumFields(section, &num_fields);
715 for (
int f = 0;
f < num_fields; ++
f) {
716 const char *field_name;
717 CHKERR PetscSectionGetFieldName(section,
f, &field_name);
719 <<
"Field " <<
f <<
" " << std::string(field_name);
728 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
734 CHKERR TSSetFromOptions(ts);
736 CHKERR set_section_monitor(ts);
739 CHKERR TSGetTime(ts, &ftime);
744 int main(
int argc,
char *argv[]) {
751 auto core_log = logging::core::get();
759 DMType dm_name =
"DMMOFEM";
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
CoordinateTypes
Coodinate system.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
constexpr int U_FIELD_DIM
FTensor::Index< 'j', SPACE_DIM > j
constexpr CoordinateTypes coord_type
select coordinate system <CARTESIAN, CYLINDRICAL>;
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, U_FIELD_DIM > OpDomainSourceU
FTensor::Index< 'k', SPACE_DIM > k
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
FTensor::Index< 'l', SPACE_DIM > l
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, 1 > OpDomainMassH
OpDomainMassH OpDomainMassP
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1, 1 > OpBaseTimesScalarField
constexpr double rho_diff
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, 1 > OpDomainSourceH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM, coord_type > OpMixScalarTimesDiv
auto smartCreateDMVector
Get smart vector from DM.
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
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.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
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.
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
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.
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
const double D
diffusivity
const double r
rate factor
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
PipelineManager::FaceEle DomainEle
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode runProblem()
[Run programme]
FreeSurface(MoFEM::Interface &m_field)
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode readMesh()
[Run programme]
MoFEMErrorCode assembleSystem()
[Boundary condition]
MoFEM::Interface & mField
boost::shared_ptr< FEMethod > domianLhsFEPtr
MoFEMErrorCode solveSystem()
[Solve]
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.
default operator for EDGE element
Data on single entity (This is passed as argument to DataOperator::doWork)
structure for User Loop Methods on finite elements
default operator for TRI element
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.
Interface for managing meshsets containing materials and boundary conditions.
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Get field gradients at integration pts for scalar filed rank 0, i.e. 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.
Approximate field valuse for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
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 with switches.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
[Push operators to pipeline]
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< VectorDouble > liftVec
boost::shared_ptr< BoundaryEle > liftFE
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc, std::pair< boost::shared_ptr< BoundaryEle >, boost::shared_ptr< VectorDouble >> p)