#include <boost/math/quadrature/gauss_kronrod.hpp>
using namespace boost::math::quadrature;
};
constexpr double omega = 7.292 * 1e-5;
constexpr double g = 9.80616;
constexpr double mu = 1e4;
constexpr double h0 = 1e4;
constexpr double h_hat = 120;
constexpr double u_max = 80;
constexpr double phi_0 = M_PI / 7;
constexpr double phi_2 = M_PI / 4;
boost::shared_ptr<MatrixDouble> u_dot_ptr,
boost::shared_ptr<MatrixDouble> grad_u_ptr,
boost::shared_ptr<MatrixDouble> grad_h_ptr)
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_dot_u = getFTensor1FromMat<3>(*
uDotPtr);
auto t_u = getFTensor1FromMat<3>(*
uPtr);
auto t_grad_u = getFTensor2FromMat<3, 3>(*
uGradPtr);
auto t_grad_h = getFTensor1FromMat<3>(*
hGradPtr);
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_normal = getFTensor1NormalsAtGaussPts();
const double alpha = t_w * vol;
auto t_nf = getFTensor1FromArray<3, 3>(
locF);
const auto a = std::sqrt(t_coords(
i) * t_coords(
i));
const auto sin_fi = t_coords(2) /
a;
const auto f = 2 *
omega * sin_fi;
t_P(
i,
j) = t_r(
i) * t_r(
j);
t_rhs(
m) = t_Q(
m,
i) * (t_dot_u(
i) + t_grad_u(
i,
j) * t_u(
j) +
f * t_A(
i,
j) * t_u(
j) +
g * t_grad_h(
i));
t_rhs_grad(
m,
j) = t_Q(
m,
i) * (
mu * t_grad_u(
i,
j));
t_rhs(
m) += t_P(
m,
j) * t_u(
j);
int rr = 0;
for (; rr !=
nbRows / 3; ++rr) {
t_nf(
i) +=
alpha * t_row_base * t_rhs(
i);
t_nf(
i) +=
alpha * t_row_diff_base(
j) * t_rhs_grad(
i,
j);
++t_row_base;
++t_row_diff_base;
++t_nf;
}
++t_row_base;
++t_row_diff_base;
}
++t_w;
++t_u;
++t_dot_u;
++t_grad_u;
++t_grad_h;
++t_coords;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble>
uPtr;
boost::shared_ptr<MatrixDouble>
uDotPtr;
boost::shared_ptr<MatrixDouble>
uGradPtr;
boost::shared_ptr<MatrixDouble>
hGradPtr;
};
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)
this->sYmm = false;
}
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_normal = getFTensor1NormalsAtGaussPts();
auto t_u = getFTensor1FromMat<3>(*
uPtr);
auto t_grad_u = getFTensor2FromMat<3, 3>(*
uGradPtr);
auto get_t_mat = [&](const int rr) {
};
const auto ts_a = getFEMethod()->ts_a;
const auto a = std::sqrt(t_coords(
i) * t_coords(
i));
const auto sin_fi = t_coords(2) /
a;
const auto f = 2 *
omega * sin_fi;
t_P(
i,
j) = t_r(
i) * t_r(
j);
t_Q(
m,
i) * (ts_a *
t_kd(
i,
j) + t_grad_u(
i,
j) +
f * t_A(
i,
j)) +
const double alpha = t_w * vol;
int rr = 0;
for (; rr !=
nbRows / 3; rr++) {
auto t_mat = get_t_mat(3 * rr);
for (
int cc = 0; cc !=
nbCols / 3; cc++) {
t_mat(
i,
j) += (
alpha * t_row_base * t_col_base) * t_rhs_du(
i,
j);
(t_row_diff_base(
m) * t_col_diff_base(
m));
++t_col_diff_base;
++t_col_base;
++t_mat;
}
++t_row_base;
++t_row_diff_base;
}
++t_row_base;
++t_row_diff_base;
}
++t_w;
++t_coords;
++t_normal;
++t_u;
++t_grad_u;
}
}
private:
boost::shared_ptr<MatrixDouble>
uPtr;
boost::shared_ptr<MatrixDouble>
uGradPtr;
};
OpULhs_dH(
const std::string field_name_row,
const std::string field_name_col)
this->sYmm = false;
}
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_normal = getFTensor1NormalsAtGaussPts();
auto get_t_vec = [&](const int rr) {
};
t_P(
i,
j) = t_r(
i) * t_r(
j);
const double alpha = t_w * vol;
int rr = 0;
for (; rr !=
nbRows / 3; rr++) {
auto t_vec = get_t_vec(3 * rr);
const double a =
alpha *
g * t_row_base;
for (
int cc = 0; cc !=
nbCols; cc++) {
t_vec(
i) +=
a * (t_Q(
i,
m) * t_col_diff_base(
m));
++t_vec;
++t_col_diff_base;
}
++t_row_base;
}
++t_row_base;
++t_w;
++t_normal;
}
}
};
private:
};
}
}
}
}
PetscBool is_restart = PETSC_FALSE;
PETSC_NULL);
auto restart_vector = [&]() {
<< "reading vector in binary from vector.dat ...";
PetscViewer viewer;
PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_READ,
&viewer);
};
if (is_restart) {
} else {
const double e_n = exp(-4 / pow(
phi_1 -
phi_0, 2));
auto get_phi = [&](const double x, const double y, const double z) {
const double r = std::sqrt(t_r(
i) * t_r(
i));
return atan2(z, r);
};
auto init_u_phi = [&](
const double phi) {
} else {
return 0.;
}
};
auto init_u = [&](
const double x,
const double y,
const double z) {
const double u_phi = init_u_phi(get_phi(x, y, z));
if (u_phi > 0) {
t_u(
i) = ((t_A(
i,
j) * t_r(
j)) * u_phi);
}
return t_u;
};
auto init_h = [&](
const double x,
const double y,
const double z) {
const double a = std::sqrt(x * x + y * y + z * z);
auto integral = [&](const double fi) {
const double u_phi = init_u_phi(fi);
const auto f = 2 *
omega * sin(fi);
return a * u_phi * (
f + (tan(fi) /
a) * u_phi);
};
auto montain = [&](
const double lambda,
const double fi) {
else
return 0.;
};
const double fi = get_phi(x, y, z);
const double lambda = atan2(x, y);
double h1 = 0;
h1 = gauss_kronrod<double, 32>::integrate(
integral,
phi_0, fi, 0, std::numeric_limits<float>::epsilon());
return h0 + (montain(
lambda, fi) - (h1 /
g));
};
auto set_domain_general = [&](auto &pipeline) {
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
};
auto set_domain_rhs = [&](auto &pipeline) {
};
auto set_domain_lhs = [&](auto &pipeline) {
pipeline.push_back(
new OpMassUU(
"U",
"U", [](
double,
double,
double) {
return 1; }));
pipeline.push_back(
new OpMassHH(
"H",
"H", [](
double,
double,
double) {
return 1; }));
};
auto 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>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
auto u_ptr = boost::make_shared<MatrixDouble>();
auto h_ptr = boost::make_shared<VectorDouble>();
auto pos_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(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{{"H", h_ptr}},
{{"U", u_ptr}, {"HO_POSITIONS", pos_ptr}},
{}, {}
)
);
CHKERR post_proc_fe->writeFile(
"out_init.h5m");
};
auto solve_init = [&]() {
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
auto name_prb =
simple->getProblemName();
name_prb,
ROW,
"U", 0, 3, is_u);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
CHKERR PCFieldSplitSetType(pc, PC_COMPOSITE_ADDITIVE);
}
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
};
};
set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
pipeline_mng->getOpDomainRhsPipeline().clear();
pipeline_mng->getOpDomainLhsPipeline().clear();
}
}
auto set_domain_general = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
};
auto set_domain_rhs = [&](auto &pipeline) {
auto dot_u_ptr = boost::make_shared<MatrixDouble>();
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto div_u_ptr = boost::make_shared<VectorDouble>();
auto dot_h_ptr = boost::make_shared<VectorDouble>();
auto grad_h_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
pipeline.push_back(
"H", dot_h_ptr, [](double, double, double) { return 1.; }));
"H", div_u_ptr, [](
double,
double,
double) {
return h0; }));
pipeline.push_back(
new OpURhs(
"U", u_ptr, dot_u_ptr, grad_u_ptr, grad_h_ptr));
};
auto set_domain_lhs = [&](auto &pipeline) {
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto grad_h_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
pipeline.push_back(
new OpMassHH(
"H",
"H", [&](
double,
double,
double) {
}));
"H",
"U", [](
const double,
const double,
const double) {
return h0; },
false, false));
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(
new OpULhs_dU(
"U",
"U", u_ptr, grad_u_ptr));
};
};
set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
}
"out_step_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
<< "writing vector in binary to vector.dat ...";
PetscViewer viewer;
PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_WRITE,
&viewer);
PetscViewerDestroy(&viewer);
}
}
private:
boost::shared_ptr<PostProcEle>
postProc;
};
auto set_initial_step = [&](auto ts) {
int step = 0;
PETSC_NULL);
CHKERR TSSetStepNumber(ts, step);
};
auto set_fieldsplit_preconditioner_ksp = [&](auto ksp) {
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
auto name_prb =
simple->getProblemName();
name_prb,
ROW,
"U", 0, 3, is_u);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
}
};
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>();
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(
auto u_ptr = boost::make_shared<MatrixDouble>();
auto h_ptr = boost::make_shared<VectorDouble>();
auto pos_ptr = boost::make_shared<MatrixDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto grad_h_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(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"H", h_ptr}},
{{"U", u_ptr}, {"HO_POSITIONS", pos_ptr}, {"GRAD_H", grad_h_ptr}},
{{"GRAD_U", grad_u_ptr}}, {}
)
);
return post_proc_fe;
};
auto set_fieldsplit_preconditioner_ts = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
CHKERR set_fieldsplit_preconditioner_ksp(ksp);
};
ts = pipeline_mng->createTSIM();
boost::shared_ptr<FEMethod> null_fe;
auto monitor_ptr = boost::make_shared<Monitor>(dm, get_fe_post_proc());
null_fe, monitor_ptr);
double ftime = 1;
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
SCATTER_FORWARD);
CHKERR set_fieldsplit_preconditioner_ts(ts);
}
}
}
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "SW"));
LogManager::setLog("SW");
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 >::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 ...
constexpr double omega
Save field DOFS on vertices/tags.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'm', SPACE_DIM > m
auto init_h
Initialisation function.
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::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 createDMVector(DM dm)
Get smart vector from DM.
#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
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 > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
const double u0
inital vale on blocksets
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpMassHH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 3 > OpSourceU
constexpr double beta_montain
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpMassUU
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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBaseTimesDotH
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()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
[Solve]
MoFEMErrorCode createCommonData()
[Set up problem]
MoFEMErrorCode runProblem()
[Run problem]
MoFEM::Interface & mField
boost::shared_ptr< FEMethod > domianRhsFEPtr
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode outputResults()
[Solve]
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.
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 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.
Calculate normals at Gauss points of triangle element.
Post post-proc data at points from hash maps.
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.
[Push operators to pipeline]
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.
boost::shared_ptr< MatrixDouble > uPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< MatrixDouble > uGradPtr
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