v0.14.0
reaction_diffusion.cpp
/**
* \file reaction_diffusion.cpp
* \example reaction_diffusion.cpp
*
**/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
const double D = 2e-3; ///< diffusivity
const double r = 1; ///< rate factor
const double k = 1; ///< caring capacity
const double u0 = 0.1; ///< inital vale on blocksets
const int save_every_nth_step = 1;
/**
* @brief Common data
*
* Common data are used to keep and pass data between elements
*
*/
struct CommonData {
MatrixDouble grad; ///< Gradients of field "u" at integration points
VectorDouble val; ///< Values of field "u" at integration points
VectorDouble dot_val; ///< Rate of values of field "u" at integration points
SmartPetscObj<Mat> M; ///< Mass matrix
SmartPetscObj<KSP> ksp; ///< Linear solver
};
/**
* @brief Assemble mass matrix
*/
struct OpAssembleMass : OpEle {
OpAssembleMass(boost::shared_ptr<CommonData> &data)
: OpEle("u", "u", OpEle::OPROWCOL), commonData(data) {
sYmm = true;
}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type, EntData &row_data,
EntData &col_data) {
const int nb_row_dofs = row_data.getIndices().size();
const int nb_col_dofs = col_data.getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
const int nb_integration_pts = getGaussPts().size2();
mat.resize(nb_row_dofs, nb_col_dofs, false);
mat.clear();
auto t_row_base = row_data.getFTensor0N();
auto t_w = getFTensor0IntegrationWeight();
const double vol = getMeasure();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double a = t_w * vol;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
auto t_col_base = col_data.getFTensor0N(gg, 0);
for (int cc = 0; cc != nb_col_dofs; ++cc) {
mat(rr, cc) += a * t_row_base * t_col_base;
++t_col_base;
}
++t_row_base;
}
++t_w;
}
CHKERR MatSetValues(commonData->M, row_data, col_data, &mat(0, 0),
ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
transMat.resize(nb_col_dofs, nb_row_dofs, false);
noalias(transMat) = trans(mat);
CHKERR MatSetValues(commonData->M, col_data, row_data, &transMat(0, 0),
ADD_VALUES);
}
}
}
private:
MatrixDouble mat, transMat;
boost::shared_ptr<CommonData> commonData;
};
/**
* @brief Assemble slow part
*
* Solve problem \f$ F(t,u,\dot{u}) = G(t,u) \f$ where here the right hand side
* \f$ G(t,u) \f$ is implemented.
*
*/
struct OpAssembleSlowRhs : OpEle {
OpAssembleSlowRhs(boost::shared_ptr<CommonData> &data)
: OpEle("u", OpEle::OPROW), commonData(data) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
const int nb_dofs = data.getIndices().size();
if (nb_dofs) {
vecF.resize(nb_dofs, false);
vecF.clear();
const int nb_integration_pts = getGaussPts().size2();
auto t_val = getFTensor0FromVec(commonData->val);
auto t_base = data.getFTensor0N();
auto t_w = getFTensor0IntegrationWeight();
const double vol = getMeasure();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double a = vol * t_w;
const double f = a * r * t_val * (1 - t_val / k);
for (int rr = 0; rr != nb_dofs; ++rr) {
const double b = f * t_base;
vecF[rr] += b;
++t_base;
}
++t_val;
++t_w;
}
CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
PETSC_TRUE);
CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
ADD_VALUES);
}
}
private:
boost::shared_ptr<CommonData> commonData;
};
/**
* @brief Assemble stiff part
*
* Solve problem \f$ F(t,u,\dot{u}) = G(t,u) \f$ where here the right hand side
* \f$ F(t,u,\dot{u}) \f$ is implemented.
*
*/
template <int DIM> struct OpAssembleStiffRhs : OpEle {
OpAssembleStiffRhs(boost::shared_ptr<CommonData> &data)
: OpEle("u", OpEle::OPROW), commonData(data) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
const int nb_dofs = data.getIndices().size();
if (nb_dofs) {
vecF.resize(nb_dofs, false);
vecF.clear();
const int nb_integration_pts = getGaussPts().size2();
auto t_dot_val = getFTensor0FromVec(commonData->dot_val);
auto t_grad = getFTensor1FromMat<DIM>(commonData->grad);
auto t_base = data.getFTensor0N();
auto t_diff_base = data.getFTensor1DiffN<DIM>();
auto t_w = getFTensor0IntegrationWeight();
const double vol = getMeasure();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double a = vol * t_w;
for (int rr = 0; rr != nb_dofs; ++rr) {
vecF[rr] += a * (t_base * t_dot_val + D * t_diff_base(i) * t_grad(i));
++t_diff_base;
++t_base;
}
++t_dot_val;
++t_grad;
++t_w;
}
CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
PETSC_TRUE);
CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
ADD_VALUES);
}
}
private:
boost::shared_ptr<CommonData> commonData;
};
/**
* @brief Assemble stiff part tangent
*
* Solve problem \f$ F(t,u,\dot{u}) = G(t,u) \f$ where here the right hand side
* \f$ \frac{\textrm{d} F}{\textrm{d} u^n} = a F_{\dot{u}}(t,u,\textrm{u}) +
* F_{u}(t,u,\textrm{u}) \f$ is implemented.
*
*/
template <int DIM> struct OpAssembleStiffLhs : OpEle {
OpAssembleStiffLhs(boost::shared_ptr<CommonData> &data)
: OpEle("u", "u", OpEle::OPROWCOL), commonData(data) {
sYmm = true;
}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type, EntData &row_data,
EntData &col_data) {
const int nb_row_dofs = row_data.getIndices().size();
const int nb_col_dofs = col_data.getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
mat.resize(nb_row_dofs, nb_col_dofs, false);
mat.clear();
const int nb_integration_pts = getGaussPts().size2();
auto t_row_base = row_data.getFTensor0N();
auto t_row_diff_base = row_data.getFTensor1DiffN<DIM>();
auto t_w = getFTensor0IntegrationWeight();
const double ts_a = getFEMethod()->ts_a;
const double vol = getMeasure();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double a = vol * t_w;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
auto t_col_base = col_data.getFTensor0N(gg, 0);
auto t_col_diff_base = col_data.getFTensor1DiffN<DIM>(gg, 0);
for (int cc = 0; cc != nb_col_dofs; ++cc) {
mat(rr, cc) += a * (t_row_base * t_col_base * ts_a +
D * t_row_diff_base(i) * t_col_diff_base(i));
++t_col_base;
++t_col_diff_base;
}
++t_row_base;
++t_row_diff_base;
}
++t_w;
}
CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
transMat.resize(nb_col_dofs, nb_row_dofs, false);
noalias(transMat) = trans(mat);
CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
&transMat(0, 0), ADD_VALUES);
}
}
}
private:
boost::shared_ptr<CommonData> commonData;
MatrixDouble mat, transMat;
};
/**
* @brief Monitor solution
*
* This functions is called by TS solver at the end of each step. It is used
* to output results to the hard drive.
*/
struct Monitor : public FEMethod {
Monitor(SmartPetscObj<DM> &dm, boost::shared_ptr<PostProcEle> &post_proc)
: dM(dm), postProc(post_proc){};
MoFEMErrorCode preProcess() { return 0; }
MoFEMErrorCode operator()() { return 0; }
MoFEMErrorCode postProcess() {
if (ts_step % save_every_nth_step == 0) {
CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc);
CHKERR postProc->writeFile(
"out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
}
}
private:
boost::shared_ptr<PostProcEle> postProc;
};
}; // namespace ReactionDiffusionEquation
using namespace ReactionDiffusionEquation;
int main(int argc, char *argv[]) {
// initialize petsc
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
try {
// Create moab and mofem instances
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
// Register DM Manager
DMType dm_name = "DMMOFEM";
// Simple interface
Simple *simple_interface;
CHKERR m_field.getInterface(simple_interface);
CHKERR simple_interface->getOptions();
CHKERR simple_interface->loadFile();
int order = 4; ///< approximation order
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
// add fields
1);
// set fields order
CHKERR simple_interface->setFieldOrder("u", order);
// setup problem
CHKERR simple_interface->setUp();
// Create common data structure
boost::shared_ptr<CommonData> data(new CommonData());
/// Alias pointers to data in common data structure
auto val_ptr = boost::shared_ptr<VectorDouble>(data, &data->val);
auto dot_val_ptr = boost::shared_ptr<VectorDouble>(data, &data->dot_val);
auto grad_ptr = boost::shared_ptr<MatrixDouble>(data, &data->grad);
// Create finite element instances to integrate the right-hand side of slow
// and stiff vector, and the tangent left-hand side for stiff part.
boost::shared_ptr<Ele> vol_ele_slow_rhs(new Ele(m_field));
boost::shared_ptr<Ele> vol_ele_stiff_rhs(new Ele(m_field));
boost::shared_ptr<Ele> vol_ele_stiff_lhs(new Ele(m_field));
// Push operators to integrate the slow right-hand side vector
vol_ele_slow_rhs->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("u", val_ptr));
vol_ele_slow_rhs->getOpPtrVector().push_back(new OpAssembleSlowRhs(data));
// PETSc IMAX and Explicit solver demans that g = M^-1 G is provided. So
// when the slow right-hand side vector (G) is assembled is solved for g
// vector.
auto solve_for_g = [&]() {
if (*(vol_ele_slow_rhs->vecAssembleSwitch)) {
CHKERR VecGhostUpdateBegin(vol_ele_slow_rhs->ts_F, ADD_VALUES,
SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(vol_ele_slow_rhs->ts_F, ADD_VALUES,
SCATTER_REVERSE);
CHKERR VecAssemblyBegin(vol_ele_slow_rhs->ts_F);
CHKERR VecAssemblyEnd(vol_ele_slow_rhs->ts_F);
*vol_ele_slow_rhs->vecAssembleSwitch = false;
}
CHKERR KSPSolve(data->ksp, vol_ele_slow_rhs->ts_F,
vol_ele_slow_rhs->ts_F);
};
// Add hook to the element to calculate g.
vol_ele_slow_rhs->postProcessHook = solve_for_g;
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
// Add operators to calculate the stiff right-hand side
vol_ele_stiff_rhs->getOpPtrVector().push_back(
new OpCalculateHOJac<2>(jac_ptr));
vol_ele_stiff_rhs->getOpPtrVector().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
vol_ele_stiff_rhs->getOpPtrVector().push_back(
new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
vol_ele_stiff_rhs->getOpPtrVector().push_back(
new OpCalculateScalarFieldValuesDot("u", dot_val_ptr));
vol_ele_stiff_rhs->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<2>("u", grad_ptr));
vol_ele_stiff_rhs->getOpPtrVector().push_back(
new OpAssembleStiffRhs<2>(data));
// Add operators to calculate the stiff left-hand side
vol_ele_stiff_lhs->getOpPtrVector().push_back(
new OpCalculateHOJac<2>(jac_ptr));
vol_ele_stiff_lhs->getOpPtrVector().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
vol_ele_stiff_lhs->getOpPtrVector().push_back(
new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
vol_ele_stiff_lhs->getOpPtrVector().push_back(
new OpAssembleStiffLhs<2>(data));
// Set integration rules
auto vol_rule = [](int, int, int p) -> int { return 2 * p; };
vol_ele_slow_rhs->getRuleHook = vol_rule;
vol_ele_stiff_rhs->getRuleHook = vol_rule;
vol_ele_stiff_lhs->getRuleHook = vol_rule;
// Crate element for post-processing
auto post_proc = boost::make_shared<PostProcEle>(m_field);
boost::shared_ptr<ForcesAndSourcesCore> null;
auto u_ptr = boost::make_shared<VectorDouble>();
post_proc->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("u", u_ptr));
post_proc->getOpPtrVector().push_back(
new OpPPMap(
post_proc->getPostProcMesh(), post_proc->getMapGaussPts(),
{{"u", u_ptr}},
{},
{},
{}
)
);
// Get PETSc discrete manager
auto dm = simple_interface->getDM();
// Get surface entities form blockset, set initial values in those
// blocksets. To keep it simple is assumed that inital values are on
// blockset 1
Range inner_surface;
CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
1, BLOCKSET, 2, inner_surface, true);
if (!inner_surface.empty()) {
Range inner_surface_verts;
CHKERR moab.get_connectivity(inner_surface, inner_surface_verts, false);
CHKERR m_field.getInterface<FieldBlas>()->setField(
u0, MBVERTEX, inner_surface_verts, "u");
}
}
// Get skin on the body, i.e. body boundary, and apply homogenous Dirichlet
// conditions on that boundary.
CHKERR moab.get_entities_by_dimension(0, 2, surface, false);
Skinner skin(&m_field.get_moab());
Range edges;
CHKERR skin.find_skin(0, surface, false, edges);
Range edges_part;
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
CHKERR pcomm->filter_pstatus(edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &edges_part);
Range edges_verts;
CHKERR moab.get_connectivity(edges_part, edges_verts, false);
// Since Dirichlet b.c. are essential boundary conditions, remove DOFs from
// the problem.
CHKERR m_field.getInterface<ProblemsManager>()->removeDofsOnEntities(
simple_interface->getProblemName(), "u",
unite(edges_verts, edges_part));
// Create mass matrix, calculate and assemble
CHKERR MatZeroEntries(data->M);
boost::shared_ptr<Ele> vol_mass_ele(new Ele(m_field));
vol_mass_ele->getOpPtrVector().push_back(new OpAssembleMass(data));
vol_mass_ele);
CHKERR MatAssemblyBegin(data->M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(data->M, MAT_FINAL_ASSEMBLY);
// Create and septup KSP (linear solver), we need this to calculate g(t,u) =
// M^-1G(t,u)
data->ksp = createKSP(m_field.get_comm());
CHKERR KSPSetOperators(data->ksp, data->M, data->M);
CHKERR KSPSetFromOptions(data->ksp);
CHKERR KSPSetUp(data->ksp);
// Create and setup TS solver
auto ts = createTS(m_field.get_comm());
// Use IMEX solver, i.e. implicit/explicit solver
CHKERR TSSetType(ts, TSARKIMEX);
CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
// Add element to calculate lhs of stiff part
CHKERR DMMoFEMTSSetIJacobian(dm, simple_interface->getDomainFEName(),
vol_ele_stiff_lhs, null, null);
// Add element to calculate rhs of stiff part
CHKERR DMMoFEMTSSetIFunction(dm, simple_interface->getDomainFEName(),
vol_ele_stiff_rhs, null, null);
// Add element to calculate rhs of slow (nonlinear) part
vol_ele_slow_rhs, null, null);
// Add monitor to time solver
boost::shared_ptr<Monitor> monitor_ptr(new Monitor(dm, post_proc));
CHKERR DMMoFEMTSSetMonitor(dm, ts, simple_interface->getDomainFEName(),
monitor_ptr, null, null);
// Create solution vector
CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_FORWARD);
// Solve problem
double ftime = 1;
CHKERR TSSetDM(ts, dm);
CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
CHKERR TSSetSolution(ts, X);
CHKERR TSSetFromOptions(ts);
CHKERR TSSolve(ts, X);
}
// finish work cleaning memory, getting statistics, etc.
return 0;
}
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
main
int main(int argc, char *argv[])
Definition: reaction_diffusion.cpp:296
surface
Definition: surface.py:1
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MoFEM::createTS
auto createTS(MPI_Comm comm)
Definition: PetscSmartObj.hpp:249
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::OpSetHOInvJacToScalarBases< 2 >
Definition: HODataOperators.hpp:78
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEM.hpp
ReactionDiffusionEquation
Definition: reaction_diffusion.cpp:12
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:261
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROWCOL
@ OPROWCOL
operator doWork is executed on FE rows &columns
Definition: ForcesAndSourcesCore.hpp:569
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::OpCalculateHOJac< 2 >
Definition: HODataOperators.hpp:273
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1167
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:747
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
ReactionDiffusionEquation::EntData
EntitiesFieldData::EntData EntData
Definition: reaction_diffusion.cpp:16
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1197
save_every_nth_step
int save_every_nth_step
Definition: photon_diffusion.cpp:67
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: FaceElementForcesAndSourcesCore.hpp:86
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::DMMoFEMTSSetIJacobian
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, 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 TS Jacobian evaluation function
Definition: DMMoFEM.cpp:853
ReactionDiffusionEquation::u0
const double u0
inital vale on blocksets
Definition: reaction_diffusion.cpp:24
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
help
static char help[]
Definition: reaction_diffusion.cpp:10
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
convert.type
type
Definition: convert.py:64
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
OpEle
FaceElementForcesAndSourcesCore::UserDataOperator OpEle
Definition: quad_polynomial_approximation.cpp:14
MoFEM::Simple::addDomainField
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:264
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
surface.surface
def surface(x, y, z, eta)
Definition: surface.py:3
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:368
ReactionDiffusionEquation::Ele
FaceElementForcesAndSourcesCore Ele
Definition: reaction_diffusion.cpp:14
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
FTensor::Index< 'i', DIM >
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:545
Range
MoFEM::DMMoFEMTSSetRHSFunction
PetscErrorCode DMMoFEMTSSetRHSFunction(DM dm, 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 TS the right hand side function
Definition: DMMoFEM.cpp:882
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
Ele
MoFEM::PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomainEle Ele
Definition: hdiv_check_approx_in_3d.cpp:18
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
Monitor
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:783
CommonData
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:22
MoFEM::DMMoFEMTSSetMonitor
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.
Definition: DMMoFEM.cpp:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::DMMoFEMTSSetIFunction
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMoFEM.cpp:800
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3249
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:21
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
ReactionDiffusionEquation::k
const double k
caring capacity
Definition: reaction_diffusion.cpp:22
MoFEM::SmartPetscObj< Mat >
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:683
MoFEM::MeshsetsManager::checkMeshset
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
Definition: MeshsetsManager.cpp:360
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:389
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698