You can find many resources on the internet about shallow wave equations. We point to two of them, which has been used as a starting point for this tutorial.
#include <boost/math/quadrature/gauss_kronrod.hpp>
using namespace boost::math::quadrature;
};
using OpBaseGradH = FormsIntegrators<DomainEleOp>::Assembly<
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;
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)
uPtr(u_ptr), uDotPtr(u_dot_ptr), uGradPtr(grad_u_ptr),
hGradPtr(grad_h_ptr) {}
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
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();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
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;
}
for (; rr < nbRowBaseFunctions; ++rr) {
++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)
uPtr(u_ptr), uGradPtr(grad_u_ptr) {
this->sYmm = false;
}
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
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) {
&locMat(rr + 0, 0), &locMat(rr + 0, 1), &locMat(rr + 0, 2),
&locMat(rr + 1, 0), &locMat(rr + 1, 1), &locMat(rr + 1, 2),
&locMat(rr + 2, 0), &locMat(rr + 2, 1), &locMat(rr + 2, 2)};
};
const auto ts_a = getFEMethod()->ts_a;
for (int gg = 0; gg != nbIntegrationPts; gg++) {
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_col_base = col_data.getFTensor0N(gg, 0);
auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
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;
}
for (; rr < nbRowBaseFunctions; ++rr) {
++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_row_base = row_data.getFTensor0N();
auto t_normal = getFTensor1NormalsAtGaussPts();
auto get_t_vec = [&](const int rr) {
&locMat(rr + 0, 0),
&locMat(rr + 1, 0),
&locMat(rr + 2, 0)};
};
for (int gg = 0; gg != nbIntegrationPts; gg++) {
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);
auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
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;
}
for (; rr < nbRowBaseFunctions; ++rr)
++t_row_base;
++t_w;
++t_normal;
}
}
};
private:
boost::shared_ptr<FEMethod> domianLhsFEPtr;
boost::shared_ptr<FEMethod> domianRhsFEPtr;
};
}
}
auto simple = mField.getInterface<Simple>();
}
Simple *
simple = mField.getInterface<Simple>();
}
PetscBool is_restart = PETSC_FALSE;
PETSC_NULL);
auto restart_vector = [&]() {
auto simple = mField.getInterface<Simple>();
<< "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));
};
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>();
pipeline.push_back(new OpGetHONormalsOnFace("HO_POSITIONS"));
pipeline.push_back(new OpCalculateHOCoords("HO_POSITIONS"));
pipeline.push_back(new OpSetHOWeightsOnFace());
};
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 simple = mField.getInterface<Simple>();
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
post_proc_fe->generateReferenceElementMesh();
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(
new OpGetHONormalsOnFace("HO_POSITIONS"));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateHOCoords("HO_POSITIONS"));
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
post_proc_fe->addFieldValuesPostProc("U");
post_proc_fe->addFieldValuesPostProc("H");
post_proc_fe->addFieldValuesPostProc("HO_POSITIONS");
CHKERR post_proc_fe->writeFile(
"out_init.h5m");
};
auto solve_init = [&]() {
auto simple = mField.getInterface<Simple>();
auto pipeline_mng = mField.getInterface<PipelineManager>();
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 bc_mng = mField.getInterface<BcManager>();
auto name_prb =
simple->getProblemName();
SmartPetscObj<IS> is_u;
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
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);
};
auto pipeline_mng = mField.getInterface<PipelineManager>();
};
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>();
pipeline.push_back(new OpGetHONormalsOnFace("HO_POSITIONS"));
pipeline.push_back(new OpCalculateHOCoords("HO_POSITIONS"));
pipeline.push_back(new OpSetHOWeightsOnFace());
pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
};
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(new OpCalculateVectorFieldValuesDot<3>("U", dot_u_ptr));
pipeline.push_back(new OpCalculateVectorFieldValues<3>("U", u_ptr));
pipeline.push_back(
pipeline.push_back(
new OpCalculateDivergenceVectorFieldValues<3>("U", div_u_ptr));
pipeline.push_back(new OpCalculateScalarFieldGradient<3>("H", grad_h_ptr));
"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(new OpCalculateVectorFieldValues<3>("U", u_ptr));
pipeline.push_back(
pipeline.push_back(new OpCalculateScalarFieldGradient<3>("H", grad_h_ptr));
pipeline.push_back(
new OpMassHH(
"H",
"H", [&](
double,
double,
double) {
return domianLhsFEPtr->ts_a;
}));
"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));
};
auto *pipeline_mng = mField.getInterface<PipelineManager>();
};
set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
domianLhsFEPtr = pipeline_mng->getDomainLhsFE();
domianRhsFEPtr = pipeline_mng->getDomainRhsFE();
}
Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc)
: dM(dm), postProc(post_proc){};
"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);
VecView(ts_u, viewer);
PetscViewerDestroy(&viewer);
}
}
private:
SmartPetscObj<DM> dM;
boost::shared_ptr<PostProcEle> postProc;
};
auto *
simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
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 bc_mng = mField.getInterface<BcManager>();
auto name_prb =
simple->getProblemName();
SmartPetscObj<IS> is_u;
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
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);
post_proc_fe->generateReferenceElementMesh();
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(
new OpGetHONormalsOnFace("HO_POSITIONS"));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateHOCoords("HO_POSITIONS"));
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
post_proc_fe->addFieldValuesPostProc("U");
post_proc_fe->addFieldValuesPostProc("H");
post_proc_fe->addFieldValuesGradientPostProc("U");
post_proc_fe->addFieldValuesGradientPostProc("H");
post_proc_fe->addFieldValuesPostProc("HO_POSITIONS");
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();
try {
DMType dm_name = "DMMOFEM";
}
}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Tensor1< T, Tensor_Dim > normalize()
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
#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 >::OpMixScalarTimesDiv< SPACE_DIM, coord_type > OpMixScalarTimesDiv
auto smartCreateDMVector
Get smart vector from DM.
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.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face 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.
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)
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
const double r
rate factor
const double u0
inital vale on blocksets
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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBaseTimesDotH
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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpConvectiveTermLhsDy< 1, 1, 3 > OpConvectiveH_dGradH
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()
[Boundary condition]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
[Solve]
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode runProblem()
[Create common data]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode outputResults()
[Solve]
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)
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.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
[Push operators to pipeline]
static char help[] =
"...\n\n";
#include <BasicFiniteElements.hpp>
};
constexpr
double a2 =
a *
a;
constexpr
double a4 =
a2 *
a2;
constexpr
double A = 6371220;
auto res_J = [](
const double x,
const double y,
const double z) {
const double res = (x * x + y * y + z * z -
a2);
return res;
};
auto res_J_dx = [](
const double x,
const double y,
const double z) {
const double res =
res_J(x, y, z);
res * (2 * z)};
};
auto lhs_J_dx2 = [](
const double x,
const double y,
const double z) {
const double res =
res_J(x, y, z);
(res * 2 + (4 * x * x)),
(4 * y * x),
(4 * z * x),
(4 * x * y),
(2 * res + (4 * y * y)),
(4 * z * y),
(4 * x * z),
(4 * y * z),
(2 * res + (4 * z * z))};
};
OpRhs(
const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr,
boost::shared_ptr<MatrixDouble> dot_x_ptr)
xPtr(x_ptr), xDotPtr(dot_x_ptr) {}
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
auto t_x0 = getFTensor1CoordsAtGaussPts();
auto t_x = getFTensor1FromMat<3>(*xPtr);
auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
t_P(
i,
j) = t_n(
i) * t_n(
j);
auto t_J_res =
res_J_dx(t_x(0), t_x(1), t_x(2));
const double alpha = t_w;
auto t_nf = getFTensor1FromArray<3, 3>(locF);
double l = std::sqrt(t_normal(
i) * t_normal(
i));
alpha *
l * ((t_P(
i,
k) * t_J_res(
k) + t_Q(
i,
k) * t_dot_x(
k)));
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
t_nf(
j) += t_row_base * t_res(
j);
++t_row_base;
++t_nf;
}
for (; rr < nbRowBaseFunctions; ++rr) {
++t_row_base;
}
++t_w;
++t_x;
++t_dot_x;
++t_x0;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble> xPtr;
boost::shared_ptr<MatrixDouble> xDotPtr;
};
OpLhs(
const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr,
boost::shared_ptr<MatrixDouble> dot_x_ptr)
xPtr(x_ptr), xDotPtr(dot_x_ptr) {
this->sYmm = false;
}
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
auto t_x0 = getFTensor1CoordsAtGaussPts();
auto t_x = getFTensor1FromMat<3>(*xPtr);
auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
auto get_t_mat = [&](const int rr) {
&locMat(rr + 0, 0), &locMat(rr + 0, 1), &locMat(rr + 0, 2),
&locMat(rr + 1, 0), &locMat(rr + 1, 1), &locMat(rr + 1, 2),
&locMat(rr + 2, 0), &locMat(rr + 2, 1), &locMat(rr + 2, 2)};
};
const double ts_a = getTSa();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
t_P(
i,
j) = t_n(
i) * t_n(
j);
auto t_J_lhs =
lhs_J_dx2(t_x(0), t_x(1), t_x(2));
double l = std::sqrt(t_normal(
i) * t_normal(
i));
const double alpha = t_w;
(
alpha *
l) * (t_P(
i,
k) * t_J_lhs(
k,
j) + t_Q(
i,
j) * ts_a);
int rr = 0;
for (; rr != nbRows / 3; rr++) {
auto t_col_base = col_data.getFTensor0N(gg, 0);
auto t_mat = get_t_mat(3 * rr);
for (int cc = 0; cc != nbCols / 3; cc++) {
const double rc = t_row_base * t_col_base;
t_mat(
i,
j) += rc * t_lhs(
i,
j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
for (; rr < nbRowBaseFunctions; ++rr)
++t_row_base;
++t_w;
++t_x;
++t_x0;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble> xPtr;
boost::shared_ptr<MatrixDouble> xDotPtr;
};
OpError(
const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr)
xPtr(x_ptr) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
auto t_w = getFTensor0IntegrationWeight();
auto t_x = getFTensor1FromMat<3>(*xPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
auto nb_integration_pts = getGaussPts().size2();
double error = 0;
for (int gg = 0; gg != nb_integration_pts; gg++) {
double l = std::sqrt(t_normal(
i) * t_normal(
i));
error += t_w *
l * std::abs((t_x(
i) * t_x(
i) -
A *
A));
++t_w;
++t_x;
++t_normal;
}
CHKERR VecSetValue(errorVec, 0, error, ADD_VALUES);
}
static SmartPetscObj<Vec> errorVec;
private:
boost::shared_ptr<MatrixDouble> xPtr;
};
private:
};
}
}
auto simple = mField.getInterface<Simple>();
}
auto simple = mField.getInterface<Simple>();
}
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
};
auto x_ptr = boost::make_shared<MatrixDouble>();
auto dot_x_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto def_ops = [&](auto &pipeline) {
pipeline.push_back(
new OpCalculateVectorFieldValues<3>("HO_POSITIONS", x_ptr));
pipeline.push_back(
new OpCalculateVectorFieldValuesDot<3>("HO_POSITIONS", dot_x_ptr));
};
def_ops(pipeline_mng->getOpDomainRhsPipeline());
def_ops(pipeline_mng->getOpDomainLhsPipeline());
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpRhs(
"HO_POSITIONS", x_ptr, dot_x_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpLhs(
"HO_POSITIONS", x_ptr, dot_x_ptr));
}
Projection10NodeCoordsOnField ent_method_material(mField, "HO_POSITIONS");
CHKERR mField.loop_dofs(
"HO_POSITIONS", ent_method_material);
auto *
simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
ts = pipeline_mng->createTSIM();
double ftime = 1;
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
SCATTER_FORWARD);
CHKERR mField.getInterface<FieldBlas>()->fieldScale(
A,
"HO_POSITIONS");
}
auto simple = mField.getInterface<Simple>();
auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
post_proc_fe->generateReferenceElementMesh();
post_proc_fe->addFieldValuesPostProc("HO_POSITIONS");
CHKERR post_proc_fe->writeFile(
"out_approx.h5m");
auto error_fe = boost::make_shared<DomainEle>(mField);
auto x_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
error_fe->getOpPtrVector().push_back(
new OpGetHONormalsOnFace("HO_POSITIONS"));
error_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>("HO_POSITIONS", x_ptr));
error_fe->getOpPtrVector().push_back(
new OpError(
"HO_POSITIONS", x_ptr));
error_fe->preProcessHook = [&]() {
MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Create vec ";
mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
};
error_fe->postProcessHook = [&]() {
double error2;
<<
"Error " << std::sqrt(error2 / (4 * M_PI *
A *
A));
};
if (mField.get_comm_size() > 1)
CHKERR mField.get_moab().write_file(
"out_ho_mesh.h5m",
"MOAB",
"PARALLEL=WRITE_PART");
else
CHKERR mField.get_moab().write_file(
"out_ho_mesh.h5m");
}
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(
try {
DMType dm_name = "DMMOFEM";
}
}
FTensor::Index< 'k', 3 > k
auto createSmartVectorMPI
Create MPI Vector.
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode runProblem()
[Run programme]
MoFEMErrorCode setOPs()
[Set up problem]
MoFEMErrorCode getOptions()
[Run programme]
MoFEMErrorCode outputResults()
[Solve]
static SmartPetscObj< Vec > errorVec