Using PipelineManager interface calculate the divergence of base functions, and integral of flux on the boundary. Since the h-div space is used, volume integral and boundary integral should give the same result.
static char help[] =
"...\n\n";
};
constexpr double a0 = 0.98;
constexpr double rho_m = 0.998;
constexpr double mu_m = 0.0101;
constexpr double rho_p = 0.0012;
constexpr double mu_p = 0.000182;
constexpr double lambda = 7.4;
constexpr double W = 0.25;
template <
int T>
constexpr int powof2() {
return 1;
else
};
constexpr double h = 0.025;
constexpr double eta =
h;
constexpr double md = 1e-2;
constexpr double eps = 1e-12;
constexpr double tol = std::numeric_limits<float>::epsilon();
return 2 * M_PI * r;
else
return 1.;
};
auto my_max = [](
const double x) {
return (x - 1 + std::abs(x + 1)) / 2; };
auto my_min = [](
const double x) {
return (x + 1 - std::abs(x - 1)) / 2; };
return 1.;
else
return 0.;
};
auto phase_function = [](
const double h,
const double diff,
const double ave) {
};
};
auto get_f = [](
const double h) {
return 4 *
W *
h * (
h *
h - 1); };
auto get_f_dh = [](
const double h) {
return 4 *
W * (3 *
h *
h - 1); };
auto get_M2 = [](
auto h_tmp) {
};
};
auto get_M3 = [](
auto h_tmp) {
const double h3 = h2 *
h;
return md * (2 * h3 - 3 * h2 + 1);
else
return md * (-2 * h3 - 3 * h2 + 1);
};
else
};
auto get_D = [](
const double A) {
return t_D;
};
constexpr double R = 0.0125;
constexpr double A = R * 0.2;
const double theta = atan2(r, y);
const double w = R +
A * cos(
n * theta);
const double d = std::sqrt(r * r + y * y);
return tanh((w - d) / (
eta * std::sqrt(2)));
};
constexpr double y0 = 0.5;
constexpr double R = 0.4;
y -= y0;
const double d = std::sqrt(r * r + y * y);
return tanh((R - d) / (
eta * std::sqrt(2)));
};
auto init_h = [](
double r,
double y,
double theta) {
};
auto save_range = [](moab::Interface &moab,
const std::string name,
CHKERR moab.create_meshset(MESHSET_SET, out_meshset);
CHKERR moab.add_entities(out_meshset, r);
CHKERR moab.write_file(name.c_str(),
"VTK",
"", &out_meshset, 1);
CHKERR moab.delete_entities(&out_meshset, 1);
};
std::vector<EntityHandle> ents_vec;
ents_vec.reserve(prb_ptr->numeredRowDofsPtr->size());
for (auto dof : *prb_ptr->numeredRowDofsPtr) {
ents_vec.push_back(dof->getEnt());
}
std::sort(ents_vec.begin(), ents_vec.end());
auto it = std::unique(ents_vec.begin(), ents_vec.end());
r.insert_list(ents_vec.begin(), it);
};
private:
};
}
}
}
auto h_ptr = boost::make_shared<VectorDouble>();
auto grad_h_ptr = boost::make_shared<MatrixDouble>();
auto g_ptr = boost::make_shared<VectorDouble>();
auto grad_g_ptr = boost::make_shared<MatrixDouble>();
auto set_generic = [&](auto &pipeline, auto &fe) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(
};
auto post_proc = [&]() {
auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
set_generic(post_proc_fe->getOpPtrVector(), post_proc_fe);
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"H", h_ptr}, {"G", g_ptr}},
{{"GRAD_H", grad_h_ptr}, {"GRAD_G", grad_g_ptr}},
{},
{}
)
);
CHKERR post_proc_fe->writeFile(
"out_init.h5m");
};
auto solve_init = [&]() {
auto set_domain_rhs = [&](auto &pipeline, auto &fe) {
set_generic(pipeline, fe);
pipeline.push_back(
new OpRhsH<true>(
"H",
nullptr,
nullptr, h_ptr,
grad_h_ptr, grad_g_ptr));
pipeline.push_back(
new OpRhsG<true>(
"G", h_ptr, grad_h_ptr, g_ptr));
};
auto set_domain_lhs = [&](auto &pipeline, auto &fe) {
set_generic(pipeline, fe);
};
auto create_subdm = [&]() {
DM subdm;
CHKERR DMCreate(mField.get_comm(), &subdm);
CHKERR DMSetType(subdm,
"DMMOFEM");
};
auto subdm = create_subdm();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getDomainLhsFE().reset();
set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline(),
pipeline_mng->getDomainRhsFE());
set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline(),
pipeline_mng->getDomainLhsFE());
auto snes = pipeline_mng->createSNES(subdm);
auto set_section_monitor = [&](auto solver) {
void *))MoFEMSNESMonitorFields,
(void *)(snes_ctx_ptr.get()), nullptr);
auto section = mField.getInterface<
ISManager>()->sectionCreate(
"SUB");
PetscInt num_fields;
CHKERR PetscSectionGetNumFields(section, &num_fields);
for (
int f = 0;
f < num_fields; ++
f) {
}
};
CHKERR set_section_monitor(snes);
CHKERR SNESSetFromOptions(snes);
CHKERR SNESSolve(snes, PETSC_NULL,
D);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
};
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYMETRY",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYMETRY",
"L", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"U",
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"L",
0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"ZERO",
"L", 0, 0);
pipeline_mng->getOpDomainRhsPipeline().clear();
pipeline_mng->getOpDomainLhsPipeline().clear();
}
auto dot_u_ptr = boost::make_shared<MatrixDouble>();
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto dot_h_ptr = boost::make_shared<VectorDouble>();
auto h_ptr = boost::make_shared<VectorDouble>();
auto grad_h_ptr = boost::make_shared<MatrixDouble>();
auto g_ptr = boost::make_shared<VectorDouble>();
auto grad_g_ptr = boost::make_shared<MatrixDouble>();
auto lambda_ptr = boost::make_shared<VectorDouble>();
auto p_ptr = boost::make_shared<VectorDouble>();
auto div_u_ptr = boost::make_shared<VectorDouble>();
auto set_domain_general = [&](auto &pipeline, auto &fe) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(
grad_u_ptr));
pipeline.push_back(
"U", div_u_ptr));
pipeline.push_back(
pipeline.push_back(
};
auto set_domain_rhs = [&](auto &pipeline, auto &fe) {
set_domain_general(pipeline, fe);
pipeline.push_back(
new OpRhsU(
"U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
grad_h_ptr, g_ptr, p_ptr));
pipeline.push_back(
new OpRhsH<false>(
"H", u_ptr, dot_h_ptr, h_ptr,
grad_h_ptr, grad_g_ptr));
pipeline.push_back(
new OpRhsG<false>(
"G", h_ptr, grad_h_ptr, g_ptr));
"P", div_u_ptr, [](const double r, const double, const double) {
}));
"P", p_ptr, [](const double r, const double, const double) {
}));
};
auto set_domain_lhs = [&](auto &pipeline, auto &fe) {
set_domain_general(pipeline, fe);
pipeline.push_back(
new OpLhsU_dU(
"U", u_ptr, grad_u_ptr, h_ptr));
pipeline.push_back(
new OpLhsU_dH(
"U",
"H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
pipeline.push_back(
new OpLhsU_dG(
"U",
"G", grad_h_ptr));
pipeline.push_back(
new OpLhsH_dU(
"H",
"U", grad_h_ptr));
"P", "U",
[](const double r, const double, const double) {
},
true, false));
pipeline.push_back(
}));
};
auto set_boundary_rhs = [&](auto &pipeline, auto &fe) {
pipeline.push_back(
};
auto set_boundary_lhs = [&](auto &pipeline, auto &fe) {
};
set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline(),
pipeline_mng->getDomainRhsFE());
set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline(),
pipeline_mng->getDomainLhsFE());
set_boundary_rhs(pipeline_mng->getOpBoundaryRhsPipeline(),
pipeline_mng->getBoundaryRhsFE());
set_boundary_lhs(pipeline_mng->getOpBoundaryLhsPipeline(),
pipeline_mng->getBoundaryLhsFE());
}
boost::shared_ptr<PostProcEdgeEle> post_proc_edge,
std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
"out_step_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
"out_step_bdy_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
}
MPI_COMM_WORLD);
<< " lift vec x: " << (*liftVec)[0] << " y: " << (*liftVec)[1];
}
private:
boost::shared_ptr<PostProcEle>
postProc;
boost::shared_ptr<BoundaryEle>
liftFE;
boost::shared_ptr<VectorDouble>
liftVec;
};
auto get_fe_post_proc = [&]() {
auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto h_ptr = boost::make_shared<VectorDouble>();
auto grad_h_ptr = boost::make_shared<MatrixDouble>();
auto p_ptr = boost::make_shared<VectorDouble>();
auto g_ptr = boost::make_shared<VectorDouble>();
auto grad_g_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
grad_u_ptr));
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"H", h_ptr}, {"P", p_ptr}, {"G", g_ptr}},
{{"U", u_ptr}, {"H_GRAD", grad_h_ptr}, {"G_GRAD", grad_g_ptr}},
{{"GRAD_U", grad_u_ptr}},
{}
)
);
return post_proc_fe;
};
auto get_bdy_post_proc_fe = [&]() {
auto post_proc_fe = boost::make_shared<PostProcEdgeEle>(mField);
auto u_ptr = boost::make_shared<MatrixDouble>();
auto p_ptr = boost::make_shared<VectorDouble>();
auto lambda_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
OpPPMap::DataMapVec{{"L", lambda_ptr}, {"P", p_ptr}},
OpPPMap::DataMapMat{{"U", u_ptr}},
OpPPMap::DataMapMat(),
OpPPMap::DataMapMat()
)
);
return post_proc_fe;
};
auto get_lift_fe = [&]() {
auto fe = boost::make_shared<BoundaryEle>(mField);
auto lift_ptr = boost::make_shared<VectorDouble>();
auto p_ptr = boost::make_shared<VectorDouble>();
auto ents_ptr = boost::make_shared<Range>();
std::vector<const CubitMeshSets *> vec_ptr;
std::regex("LIFT"), vec_ptr);
for (auto m_ptr : vec_ptr) {
ents, true);
ents_ptr->merge(ents);
}
MOFEM_LOG(
"FS", Sev::noisy) <<
"Lift ents " << (*ents_ptr);
fe->getOpPtrVector().push_back(
fe->getOpPtrVector().push_back(
return std::make_pair(fe, lift_ptr);
};
auto set_ts = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
};
auto ts = pipeline_mng->createTSIM();
CHKERR TSSetType(ts, TSALPHA);
auto set_post_proc_monitor = [&](auto dm) {
boost::shared_ptr<FEMethod> null_fe;
auto monitor_ptr = boost::make_shared<Monitor>(
dm, get_fe_post_proc(), get_bdy_post_proc_fe(), get_lift_fe());
null_fe, monitor_ptr);
};
CHKERR set_post_proc_monitor(dm);
double ftime = 1;
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
SCATTER_FORWARD);
auto print_fields_in_section = [&]() {
auto section = mField.getInterface<
ISManager>()->sectionCreate(
PetscInt num_fields;
CHKERR PetscSectionGetNumFields(section, &num_fields);
for (
int f = 0;
f < num_fields; ++
f) {
}
};
CHKERR print_fields_in_section();
}
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "FS"));
LogManager::setLog("FS");
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
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
constexpr CoordinateTypes coord_type
select coordinate system <CARTESIAN, CYLINDRICAL>;
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, U_FIELD_DIM > OpDomainSourceU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
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 DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
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.
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
#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
auto marker
set bit to marker
FTensor::Index< 'i', SPACE_DIM > i
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomianParentEle DomianParentEle
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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.
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)
auto getProblemPtr(DM dm)
get problem pointer from DM
auto smartGetDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
const double r
rate factor
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr auto field_name
PostProcBrokenMeshInMoab< EdgeElementForcesAndSourcesCore > PostProcEdgeEle
OpBaseImpl< PETSC, EdgeEleOp > OpBase
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode runProblem()
[Run programme]
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)
default operator for TRI element
Base face element used to integrate on skeleton.
Section manager is used to create indexes and sections.
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 rate of scalar field at integration points.
Get value at integration points for scalar 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.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
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
boost::shared_ptr< PostProcEle > postProc