14static char help[] =
"...\n\n";
79constexpr double a0 = 0.98;
81constexpr double mu_m = 0.0101;
83constexpr double mu_p = 0.000182;
85constexpr double W = 0.25;
87template <
int T>
constexpr int powof2() {
95constexpr double h = 0.025;
100constexpr double md = 1e-2;
101constexpr double eps = 1e-12;
102constexpr double tol = std::numeric_limits<float>::epsilon();
121auto my_max = [](
const double x) {
return (x - 1 + std::abs(x + 1)) / 2; };
122auto my_min = [](
const double x) {
return (x + 1 - std::abs(x - 1)) / 2; };
125 if (
h >= -1 &&
h < 1)
139auto get_f = [](
const double h) {
return 4 *
W *
h * (
h *
h - 1); };
147 return md * (1 -
h *
h);
157 const double h2 =
h *
h;
158 const double h3 = h2 *
h;
160 return md * (2 * h3 - 3 * h2 + 1);
162 return md * (-2 * h3 - 3 * h2 + 1);
184 constexpr double R = 0.0125;
185 constexpr double A = R * 0.2;
186 const double theta = atan2(r, y);
187 const double w = R +
A * cos(
n * theta);
188 const double d = std::sqrt(r * r + y * y);
189 return tanh((w - d) / (
eta * std::sqrt(2)));
193 constexpr double y0 = 0.5;
194 constexpr double R = 0.4;
196 const double d = std::sqrt(r * r + y * y);
197 return tanh((R - d) / (
eta * std::sqrt(2)));
200auto init_h = [](
double r,
double y,
double theta) {
207auto save_range = [](moab::Interface &moab,
const std::string name,
211 CHKERR moab.create_meshset(MESHSET_SET, out_meshset);
212 CHKERR moab.add_entities(out_meshset, r);
213 CHKERR moab.write_file(name.c_str(),
"VTK",
"", &out_meshset, 1);
214 CHKERR moab.delete_entities(&out_meshset, 1);
220 std::vector<EntityHandle> ents_vec;
221 ents_vec.reserve(prb_ptr->numeredRowDofsPtr->size());
222 for (
auto dof : *prb_ptr->numeredRowDofsPtr) {
223 ents_vec.push_back(dof->getEnt());
225 std::sort(ents_vec.begin(), ents_vec.end());
226 auto it = std::unique(ents_vec.begin(), ents_vec.end());
228 r.insert_list(ents_vec.begin(), it);
325 auto dm =
simple->getDM();
327 auto h_ptr = boost::make_shared<VectorDouble>();
328 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
329 auto g_ptr = boost::make_shared<VectorDouble>();
330 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
332 auto set_generic = [&](
auto &pipeline,
auto &fe) {
333 auto det_ptr = boost::make_shared<VectorDouble>();
334 auto jac_ptr = boost::make_shared<MatrixDouble>();
335 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
351 auto post_proc = [&]() {
353 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
355 set_generic(post_proc_fe->getOpPtrVector(), post_proc_fe);
359 post_proc_fe->getOpPtrVector().push_back(
362 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
364 {{
"H", h_ptr}, {
"G", g_ptr}},
366 {{
"GRAD_H", grad_h_ptr}, {
"GRAD_G", grad_g_ptr}},
377 CHKERR post_proc_fe->writeFile(
"out_init.h5m");
382 auto solve_init = [&]() {
385 auto set_domain_rhs = [&](
auto &pipeline,
auto &fe) {
386 set_generic(pipeline, fe);
387 pipeline.push_back(
new OpRhsH<true>(
"H",
nullptr,
nullptr, h_ptr,
388 grad_h_ptr, grad_g_ptr));
389 pipeline.push_back(
new OpRhsG<true>(
"G", h_ptr, grad_h_ptr, g_ptr));
392 auto set_domain_lhs = [&](
auto &pipeline,
auto &fe) {
393 set_generic(pipeline, fe);
394 pipeline.push_back(
new OpLhsH_dH<true>(
"H",
nullptr, h_ptr, grad_g_ptr));
400 auto create_subdm = [&]() {
402 CHKERR DMCreate(mField.get_comm(), &subdm);
403 CHKERR DMSetType(subdm,
"DMMOFEM");
416 auto subdm = create_subdm();
420 pipeline_mng->getDomainRhsFE().reset();
421 pipeline_mng->getDomainLhsFE().reset();
425 set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline(),
426 pipeline_mng->getDomainRhsFE());
427 set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline(),
428 pipeline_mng->getDomainLhsFE());
431 auto snes = pipeline_mng->createSNES(subdm);
434 auto set_section_monitor = [&](
auto solver) {
436 CHKERR SNESMonitorSet(snes,
439 (
void *)(snes_ctx_ptr.get()),
nullptr);
440 auto section = mField.getInterface<
ISManager>()->sectionCreate(
"SUB");
442 CHKERR PetscSectionGetNumFields(section, &num_fields);
443 for (
int f = 0;
f < num_fields; ++
f) {
447 <<
"Field " <<
f <<
" " << std::string(
field_name);
452 CHKERR set_section_monitor(snes);
454 CHKERR SNESSetFromOptions(snes);
455 CHKERR SNESSolve(snes, PETSC_NULL,
D);
457 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
458 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
467 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYMETRY",
469 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYMETRY",
471 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"U",
473 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"L",
475 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"ZERO",
479 pipeline_mng->getOpDomainRhsPipeline().clear();
480 pipeline_mng->getOpDomainLhsPipeline().clear();
490 auto dot_u_ptr = boost::make_shared<MatrixDouble>();
491 auto u_ptr = boost::make_shared<MatrixDouble>();
492 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
493 auto dot_h_ptr = boost::make_shared<VectorDouble>();
494 auto h_ptr = boost::make_shared<VectorDouble>();
495 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
496 auto g_ptr = boost::make_shared<VectorDouble>();
497 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
498 auto lambda_ptr = boost::make_shared<VectorDouble>();
499 auto p_ptr = boost::make_shared<VectorDouble>();
500 auto div_u_ptr = boost::make_shared<VectorDouble>();
504 auto set_domain_general = [&](
auto &pipeline,
auto &fe) {
505 auto det_ptr = boost::make_shared<VectorDouble>();
506 auto jac_ptr = boost::make_shared<MatrixDouble>();
507 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
538 auto set_domain_rhs = [&](
auto &pipeline,
auto &fe) {
539 set_domain_general(pipeline, fe);
540 pipeline.push_back(
new OpRhsU(
"U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
541 grad_h_ptr, g_ptr, p_ptr));
542 pipeline.push_back(
new OpRhsH<false>(
"H", u_ptr, dot_h_ptr, h_ptr,
543 grad_h_ptr, grad_g_ptr));
544 pipeline.push_back(
new OpRhsG<false>(
"G", h_ptr, grad_h_ptr, g_ptr));
546 "P", div_u_ptr, [](
const double r,
const double,
const double) {
550 "P", p_ptr, [](
const double r,
const double,
const double) {
555 auto set_domain_lhs = [&](
auto &pipeline,
auto &fe) {
556 set_domain_general(pipeline, fe);
557 pipeline.push_back(
new OpLhsU_dU(
"U", u_ptr, grad_u_ptr, h_ptr));
559 new OpLhsU_dH(
"U",
"H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
560 pipeline.push_back(
new OpLhsU_dG(
"U",
"G", grad_h_ptr));
562 pipeline.push_back(
new OpLhsH_dU(
"H",
"U", grad_h_ptr));
571 [](
const double r,
const double,
const double) {
581 auto set_boundary_rhs = [&](
auto &pipeline,
auto &fe) {
589 auto set_boundary_lhs = [&](
auto &pipeline,
auto &fe) {
601 set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline(),
602 pipeline_mng->getDomainRhsFE());
603 set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline(),
604 pipeline_mng->getDomainLhsFE());
606 set_boundary_rhs(pipeline_mng->getOpBoundaryRhsPipeline(),
607 pipeline_mng->getBoundaryRhsFE());
608 set_boundary_lhs(pipeline_mng->getOpBoundaryLhsPipeline(),
609 pipeline_mng->getBoundaryLhsFE());
626 boost::shared_ptr<PostProcEdgeEle> post_proc_edge,
627 std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
638 "out_step_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
643 "out_step_bdy_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
653 <<
" lift vec x: " << (*liftVec)[0] <<
" y: " << (*liftVec)[1];
660 boost::shared_ptr<PostProcEle>
postProc;
672 auto dm =
simple->getDM();
674 auto get_fe_post_proc = [&]() {
675 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
677 auto det_ptr = boost::make_shared<VectorDouble>();
678 auto jac_ptr = boost::make_shared<MatrixDouble>();
679 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
681 auto u_ptr = boost::make_shared<MatrixDouble>();
682 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
683 auto h_ptr = boost::make_shared<VectorDouble>();
684 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
685 auto p_ptr = boost::make_shared<VectorDouble>();
686 auto g_ptr = boost::make_shared<VectorDouble>();
687 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
689 post_proc_fe->getOpPtrVector().push_back(
691 post_proc_fe->getOpPtrVector().push_back(
693 post_proc_fe->getOpPtrVector().push_back(
696 post_proc_fe->getOpPtrVector().push_back(
698 post_proc_fe->getOpPtrVector().push_back(
702 post_proc_fe->getOpPtrVector().push_back(
704 post_proc_fe->getOpPtrVector().push_back(
707 post_proc_fe->getOpPtrVector().push_back(
709 post_proc_fe->getOpPtrVector().push_back(
711 post_proc_fe->getOpPtrVector().push_back(
716 post_proc_fe->getOpPtrVector().push_back(
719 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
721 {{
"H", h_ptr}, {
"P", p_ptr}, {
"G", g_ptr}},
723 {{
"U", u_ptr}, {
"H_GRAD", grad_h_ptr}, {
"G_GRAD", grad_g_ptr}},
725 {{
"GRAD_U", grad_u_ptr}},
736 auto get_bdy_post_proc_fe = [&]() {
737 auto post_proc_fe = boost::make_shared<PostProcEdgeEle>(mField);
739 auto u_ptr = boost::make_shared<MatrixDouble>();
740 auto p_ptr = boost::make_shared<VectorDouble>();
741 auto lambda_ptr = boost::make_shared<VectorDouble>();
743 post_proc_fe->getOpPtrVector().push_back(
745 post_proc_fe->getOpPtrVector().push_back(
747 post_proc_fe->getOpPtrVector().push_back(
752 post_proc_fe->getOpPtrVector().push_back(
754 new OpPPMap(post_proc_fe->getPostProcMesh(),
755 post_proc_fe->getMapGaussPts(),
772 auto get_lift_fe = [&]() {
773 auto fe = boost::make_shared<BoundaryEle>(mField);
774 auto lift_ptr = boost::make_shared<VectorDouble>();
775 auto p_ptr = boost::make_shared<VectorDouble>();
776 auto ents_ptr = boost::make_shared<Range>();
778 std::vector<const CubitMeshSets *> vec_ptr;
780 std::regex(
"LIFT"), vec_ptr);
781 for (
auto m_ptr : vec_ptr) {
784 CHKERR mField.get_moab().get_entities_by_dimension(meshset,
SPACE_DIM - 1,
786 ents_ptr->merge(ents);
789 MOFEM_LOG(
"FS", Sev::noisy) <<
"Lift ents " << (*ents_ptr);
791 fe->getOpPtrVector().push_back(
793 fe->getOpPtrVector().push_back(
796 return std::make_pair(fe, lift_ptr);
799 auto set_ts = [&](
auto solver) {
802 CHKERR TSGetSNES(solver, &snes);
804 CHKERR SNESGetKSP(snes, &ksp);
808 auto ts = pipeline_mng->createTSIM();
809 CHKERR TSSetType(ts, TSALPHA);
811 auto set_post_proc_monitor = [&](
auto dm) {
813 boost::shared_ptr<FEMethod> null_fe;
814 auto monitor_ptr = boost::make_shared<Monitor>(
815 dm, get_fe_post_proc(), get_bdy_post_proc_fe(), get_lift_fe());
817 null_fe, monitor_ptr);
820 CHKERR set_post_proc_monitor(dm);
825 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
830 CHKERR TSSetSolution(ts, T);
831 CHKERR TSSetFromOptions(ts);
835 auto print_fields_in_section = [&]() {
838 auto section = mField.getInterface<
ISManager>()->sectionCreate(
839 simple->getProblemName());
841 CHKERR PetscSectionGetNumFields(section, &num_fields);
842 for (
int f = 0;
f < num_fields; ++
f) {
846 <<
"Field " <<
f <<
" " << std::string(
field_name);
851 CHKERR print_fields_in_section();
854 CHKERR TSGetTime(ts, &ftime);
859int main(
int argc,
char *argv[]) {
866 auto core_log = logging::core::get();
874 DMType dm_name =
"DMMOFEM";
879 moab::Core mb_instance;
880 moab::Interface &moab = mb_instance;
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define BITREFLEVEL_SIZE
max number of refinements
#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.
FTensor::Index< 'n', SPACE_DIM > n
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 >::OpBaseTimesScalar< 1, 1 > OpBaseTimesScalar
constexpr double rho_diff
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, 1 > OpDomainSourceH
constexpr int order
approximation order
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM, coord_type > OpMixScalarTimesDiv
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 DMMoFEMAddSubFieldRow(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
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.
auto smartCreateDMVector(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.
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomianParentEle DomianParentEle
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
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.
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto getProblemPtr(DM dm)
get problem pointer from DM
auto smartGetDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
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 makeRefProblem()
MoFEMErrorCode solveSystem()
[Solve]
boost::weak_ptr< CacheTuple > getCacheWeakPtr() const
Get the cache weak ptr object.
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.
Base face element used to integrate on skeleton.
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
Base face element used to integrate on skeleton.
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 values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Set inverse jacobian to base functions.
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.
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
boost::shared_ptr< PostProcEdgeEle > postProcEdge
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, boost::shared_ptr< PostProcEdgeEle > post_proc_edge, std::pair< boost::shared_ptr< BoundaryEle >, boost::shared_ptr< VectorDouble > > p)
boost::shared_ptr< PostProcEle > postProc