v0.14.0
Loading...
Searching...
No Matches
reaction_diffusion.cpp
/**
* \file reaction_diffusion.cpp
* \example reaction_diffusion.cpp
*
**/
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();
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:
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();
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>();
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>();
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;
};
/**
* @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; }
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";
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";
CHKERR DMRegister_MoFEM(dm_name);
// 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(
// 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(
// 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;
}
static Index< 'p', 3 > p
std::string param_file
static char help[]
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
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:786
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:509
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1183
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:572
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:839
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1153
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto createKSP(MPI_Comm comm)
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:1042
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:868
auto createTS(MPI_Comm comm)
FaceElementForcesAndSourcesCore::UserDataOperator OpEle
const double D
diffusivity
const double r
rate factor
const double u0
inital vale on blocksets
const double k
caring capacity
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FaceElementForcesAndSourcesCore Ele
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition: Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
bool sYmm
If true assume that matrix is symmetric structure.
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.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Basic algebra on fields.
Definition: FieldBlas.hpp:21
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Interface for managing meshsets containing materials and boundary conditions.
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.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
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
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:362
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:341
intrusive_ptr for managing petsc objects
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
[Push operators to pipeline]
SmartPetscObj< KSP > ksp
Linear solver.
VectorDouble dot_val
Rate of values of field "u" at integration points.
VectorDouble val
Values of field "u" at integration points.
SmartPetscObj< Mat > M
Mass matrix.
MatrixDouble grad
Gradients of field "u" at integration points.
boost::shared_ptr< PostProcEle > postProc
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
boost::shared_ptr< CommonData > commonData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.