28#include <boost/math/quadrature/gauss_kronrod.hpp>
29using namespace boost::math::quadrature;
66 GAUSS>::OpSource<1, 3>;
68 GAUSS>::OpSource<1, 1>;
77constexpr double omega = 7.292 * 1e-5;
78constexpr double g = 9.80616;
79constexpr double mu = 1e4;
80constexpr double h0 = 1e4;
84constexpr double phi_0 = M_PI / 7;
86constexpr double phi_2 = M_PI / 4;
100 boost::shared_ptr<MatrixDouble> u_dot_ptr,
101 boost::shared_ptr<MatrixDouble> grad_u_ptr,
102 boost::shared_ptr<MatrixDouble> grad_h_ptr)
110 const double vol = getMeasure();
111 auto t_w = getFTensor0IntegrationWeight();
114 auto t_dot_u = getFTensor1FromMat<3>(*
uDotPtr);
115 auto t_u = getFTensor1FromMat<3>(*
uPtr);
116 auto t_grad_u = getFTensor2FromMat<3, 3>(*
uGradPtr);
117 auto t_grad_h = getFTensor1FromMat<3>(*
hGradPtr);
118 auto t_coords = getFTensor1CoordsAtGaussPts();
119 auto t_normal = getFTensor1NormalsAtGaussPts();
123 const double alpha = t_w * vol;
124 auto t_nf = getFTensor1FromArray<3, 3>(
locF);
129 const auto a = std::sqrt(t_coords(
i) * t_coords(
i));
130 const auto sin_fi = t_coords(2) /
a;
131 const auto f = 2 *
omega * sin_fi;
134 t_r(
i) = t_normal(
i);
139 t_P(
i,
j) = t_r(
i) * t_r(
j);
145 t_rhs(
m) = t_Q(
m,
i) * (t_dot_u(
i) + t_grad_u(
i,
j) * t_u(
j) +
146 f * t_A(
i,
j) * t_u(
j) +
g * t_grad_h(
i));
147 t_rhs_grad(
m,
j) = t_Q(
m,
i) * (
mu * t_grad_u(
i,
j));
149 t_rhs(
m) += t_P(
m,
j) * t_u(
j);
152 for (; rr !=
nbRows / 3; ++rr) {
153 t_nf(
i) +=
alpha * t_row_base * t_rhs(
i);
154 t_nf(
i) +=
alpha * t_row_diff_base(
j) * t_rhs_grad(
i,
j);
177 boost::shared_ptr<MatrixDouble>
uPtr;
185 OpULhs_dU(
const std::string field_name_row,
const std::string field_name_col,
186 boost::shared_ptr<MatrixDouble> u_ptr,
187 boost::shared_ptr<MatrixDouble> grad_u_ptr)
198 const double vol = getMeasure();
199 auto t_w = getFTensor0IntegrationWeight();
202 auto t_coords = getFTensor1CoordsAtGaussPts();
203 auto t_normal = getFTensor1NormalsAtGaussPts();
205 auto t_u = getFTensor1FromMat<3>(*
uPtr);
206 auto t_grad_u = getFTensor2FromMat<3, 3>(*
uGradPtr);
208 auto get_t_mat = [&](
const int rr) {
217 const auto ts_a = getFEMethod()->ts_a;
221 const auto a = std::sqrt(t_coords(
i) * t_coords(
i));
222 const auto sin_fi = t_coords(2) /
a;
223 const auto f = 2 *
omega * sin_fi;
226 t_r(
i) = t_normal(
i);
231 t_P(
i,
j) = t_r(
i) * t_r(
j);
239 t_Q(
m,
i) * (ts_a *
t_kd(
i,
j) + t_grad_u(
i,
j) +
f * t_A(
i,
j)) +
242 const double alpha = t_w * vol;
245 for (; rr !=
nbRows / 3; rr++) {
249 auto t_mat = get_t_mat(3 * rr);
251 for (
int cc = 0; cc !=
nbCols / 3; cc++) {
252 t_mat(
i,
j) += (
alpha * t_row_base * t_col_base) * t_rhs_du(
i,
j);
254 (t_row_diff_base(
m) * t_col_diff_base(
m));
278 boost::shared_ptr<MatrixDouble>
uPtr;
284 OpULhs_dH(
const std::string field_name_row,
const std::string field_name_col)
295 const double vol = getMeasure();
297 auto t_w = getFTensor0IntegrationWeight();
301 auto t_normal = getFTensor1NormalsAtGaussPts();
303 auto get_t_vec = [&](
const int rr) {
316 t_r(
i) = t_normal(
i);
321 t_P(
i,
j) = t_r(
i) * t_r(
j);
324 const double alpha = t_w * vol;
327 for (; rr !=
nbRows / 3; rr++) {
328 auto t_vec = get_t_vec(3 * rr);
330 const double a =
alpha *
g * t_row_base;
332 for (
int cc = 0; cc !=
nbCols; cc++) {
333 t_vec(
i) +=
a * (t_Q(
i,
m) * t_col_diff_base(
m));
431 PetscBool is_restart = PETSC_FALSE;
435 auto restart_vector = [&]() {
438 auto dm =
simple->getDM();
440 <<
"reading vector in binary from vector.dat ...";
442 PetscViewerBinaryOpen(PETSC_COMM_WORLD,
"vector.dat", FILE_MODE_READ,
456 const double e_n = exp(-4 / pow(
phi_1 -
phi_0, 2));
463 auto get_phi = [&](
const double x,
const double y,
const double z) {
465 const double r = std::sqrt(t_r(
i) * t_r(
i));
469 auto init_u_phi = [&](
const double phi) {
477 auto init_u = [&](
const double x,
const double y,
const double z) {
479 const double u_phi = init_u_phi(get_phi(x, y, z));
483 t_u(
i) = ((t_A(
i,
j) * t_r(
j)) * u_phi);
488 auto init_h = [&](
const double x,
const double y,
const double z) {
489 const double a = std::sqrt(x * x + y * y + z * z);
491 auto integral = [&](
const double fi) {
492 const double u_phi = init_u_phi(fi);
493 const auto f = 2 *
omega * sin(fi);
494 return a * u_phi * (
f + (tan(fi) /
a) * u_phi);
497 auto montain = [&](
const double lambda,
const double fi) {
505 const double fi = get_phi(x, y, z);
506 const double lambda = atan2(x, y);
510 h1 = gauss_kronrod<double, 32>::integrate(
511 integral,
phi_0, fi, 0, std::numeric_limits<float>::epsilon());
513 return h0 + (montain(
lambda, fi) - (h1 /
g));
516 auto set_domain_general = [&](
auto &pipeline) {
517 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
523 auto set_domain_rhs = [&](
auto &pipeline) {
528 auto set_domain_lhs = [&](
auto &pipeline) {
530 new OpMassUU(
"U",
"U", [](
double,
double,
double) {
return 1; }));
532 new OpMassHH(
"H",
"H", [](
double,
double,
double) {
return 1; }));
535 auto post_proc = [&]() {
538 auto dm =
simple->getDM();
540 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
541 post_proc_fe->generateReferenceElementMesh();
543 auto det_ptr = boost::make_shared<VectorDouble>();
544 auto jac_ptr = boost::make_shared<MatrixDouble>();
545 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
547 post_proc_fe->getOpPtrVector().push_back(
549 post_proc_fe->getOpPtrVector().push_back(
551 post_proc_fe->getOpPtrVector().push_back(
553 post_proc_fe->getOpPtrVector().push_back(
555 post_proc_fe->addFieldValuesPostProc(
"U");
556 post_proc_fe->addFieldValuesPostProc(
"H");
557 post_proc_fe->addFieldValuesPostProc(
"HO_POSITIONS");
560 CHKERR post_proc_fe->writeFile(
"out_init.h5m");
565 auto solve_init = [&]() {
570 auto solver = pipeline_mng->createKSP();
571 CHKERR KSPSetFromOptions(solver);
573 CHKERR KSPGetPC(solver, &pc);
574 PetscBool is_pcfs = PETSC_FALSE;
575 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
576 if (is_pcfs == PETSC_TRUE) {
578 auto name_prb =
simple->getProblemName();
581 name_prb,
ROW,
"U", 0, 3, is_u);
582 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
583 CHKERR PCFieldSplitSetType(pc, PC_COMPOSITE_ADDITIVE);
588 auto dm =
simple->getDM();
593 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
594 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
608 set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
609 set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
610 set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
611 set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
617 pipeline_mng->getOpDomainRhsPipeline().clear();
618 pipeline_mng->getOpDomainLhsPipeline().clear();
631 auto set_domain_general = [&](
auto &pipeline) {
632 auto det_ptr = boost::make_shared<VectorDouble>();
633 auto jac_ptr = boost::make_shared<MatrixDouble>();
634 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
643 auto set_domain_rhs = [&](
auto &pipeline) {
644 auto dot_u_ptr = boost::make_shared<MatrixDouble>();
645 auto u_ptr = boost::make_shared<MatrixDouble>();
646 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
647 auto div_u_ptr = boost::make_shared<VectorDouble>();
648 auto dot_h_ptr = boost::make_shared<VectorDouble>();
649 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
662 "H", dot_h_ptr, [](
double,
double,
double) {
return 1.; }));
664 "H", div_u_ptr, [](
double,
double,
double) {
return h0; }));
665 pipeline.push_back(
new OpConvectiveH(
"H", u_ptr, grad_h_ptr));
667 new OpURhs(
"U", u_ptr, dot_u_ptr, grad_u_ptr, grad_h_ptr));
670 auto set_domain_lhs = [&](
auto &pipeline) {
671 auto u_ptr = boost::make_shared<MatrixDouble>();
672 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
673 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
680 pipeline.push_back(
new OpMassHH(
"H",
"H", [&](
double,
double,
double) {
684 "H",
"U", [](
const double,
const double,
const double) {
return h0; },
690 pipeline.push_back(
new OpULhs_dU(
"U",
"U", u_ptr, grad_u_ptr));
691 pipeline.push_back(
new OpULhs_dH(
"U",
"H"));
702 set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
703 set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
705 set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
706 set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
730 "out_step_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
732 <<
"writing vector in binary to vector.dat ...";
734 PetscViewerBinaryOpen(PETSC_COMM_WORLD,
"vector.dat", FILE_MODE_WRITE,
736 VecView(
ts_u, viewer);
737 PetscViewerDestroy(&viewer);
744 boost::shared_ptr<PostProcEle>
postProc;
753 auto dm =
simple->getDM();
755 auto set_initial_step = [&](
auto ts) {
760 CHKERR TSSetStepNumber(ts, step);
764 auto set_fieldsplit_preconditioner_ksp = [&](
auto ksp) {
767 CHKERR KSPGetPC(ksp, &pc);
768 PetscBool is_pcfs = PETSC_FALSE;
769 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
770 if (is_pcfs == PETSC_TRUE) {
772 auto name_prb =
simple->getProblemName();
775 name_prb,
ROW,
"U", 0, 3, is_u);
776 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
782 auto get_fe_post_proc = [&]() {
783 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
784 post_proc_fe->generateReferenceElementMesh();
786 auto det_ptr = boost::make_shared<VectorDouble>();
787 auto jac_ptr = boost::make_shared<MatrixDouble>();
788 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
790 post_proc_fe->getOpPtrVector().push_back(
792 post_proc_fe->getOpPtrVector().push_back(
794 post_proc_fe->getOpPtrVector().push_back(
796 post_proc_fe->getOpPtrVector().push_back(
798 post_proc_fe->getOpPtrVector().push_back(
800 post_proc_fe->addFieldValuesPostProc(
"U");
801 post_proc_fe->addFieldValuesPostProc(
"H");
802 post_proc_fe->addFieldValuesGradientPostProc(
"U");
803 post_proc_fe->addFieldValuesGradientPostProc(
"H");
804 post_proc_fe->addFieldValuesPostProc(
"HO_POSITIONS");
808 auto set_fieldsplit_preconditioner_ts = [&](
auto solver) {
811 CHKERR TSGetSNES(solver, &snes);
813 CHKERR SNESGetKSP(snes, &ksp);
814 CHKERR set_fieldsplit_preconditioner_ksp(ksp);
819 ts = pipeline_mng->createTSIM();
821 boost::shared_ptr<FEMethod> null_fe;
822 auto monitor_ptr = boost::make_shared<Monitor>(dm, get_fe_post_proc());
824 null_fe, monitor_ptr);
829 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
835 CHKERR TSSetFromOptions(ts);
836 CHKERR set_fieldsplit_preconditioner_ts(ts);
838 CHKERR set_initial_step(ts);
840 CHKERR TSGetTime(ts, &ftime);
862int main(
int argc,
char *argv[]) {
869 auto core_log = logging::core::get();
870 core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(),
"SW"));
871 LogManager::setLog(
"SW");
877 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
Tensor1< T, Tensor_Dim > normalize()
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#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 ...
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM, coord_type > OpMixScalarTimesDiv
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
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.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
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.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
DeprecatedCoreInterface Interface
const double D
diffusivity
const double r
rate factor
const double u0
inital vale on blocksets
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
static char help[]
[Check]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpMassHH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 3 > OpSourceU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBaseTimesDotH
constexpr double beta_montain
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpMassUU
FTensor::Index< 'l', 3 > l
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpConvectiveTermLhsDu< 1, 1, 3 > OpConvectiveH_dU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpConvectiveTermRhs< 1, 1, 3 > OpConvectiveH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< 3 > OpBaseDivU
OpBaseTimesDotH OpBaseTimesDivU
FTensor::Index< 'j', 3 > j
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpConvectiveTermLhsDy< 1, 1, 3 > OpConvectiveH_dGradH
FTensor::Index< 'i', 3 > i
FTensor::Index< 'm', 3 > m
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpSourceH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixVectorTimesGrad< 1, 3, 3 > OpBaseGradH
constexpr double alpha_montain
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
boost::shared_ptr< FEMethod > domianRhsFEPtr
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
boost::shared_ptr< FEMethod > domianLhsFEPtr
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.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
structure for User Loop Methods on finite elements
default operator for TRI element
Section manager is used to create indexes and sections.
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Calculate HO coordinates at gauss points.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector 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.
Calculate normals at Gauss points of triangle element.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
Simple interface for fast problem set-up.
PetscInt ts_step
time step
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc)
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< PostProcEle > postProc
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
OpULhs_dH(const std::string field_name_row, const std::string field_name_col)
boost::shared_ptr< MatrixDouble > uPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< MatrixDouble > uGradPtr
OpULhs_dU(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > grad_u_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
Class dedicated to integrate operator.
boost::shared_ptr< MatrixDouble > uDotPtr
boost::shared_ptr< MatrixDouble > uGradPtr
boost::shared_ptr< MatrixDouble > hGradPtr
boost::shared_ptr< MatrixDouble > uPtr
OpURhs(const std::string field_name, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > u_dot_ptr, boost::shared_ptr< MatrixDouble > grad_u_ptr, boost::shared_ptr< MatrixDouble > grad_h_ptr)