22static char help[] =
"...\n\n";
56constexpr double c = 1;
57constexpr double k = 1;
81 "out_level_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
88 boost::shared_ptr<PostProcEle>
postProc;
170 1,
BLOCKSET, 2, inner_surface,
true);
171 if (!inner_surface.empty()) {
172 Range inner_surface_verts;
174 inner_surface_verts,
false);
176 init_u, MBVERTEX, inner_surface_verts,
"U");
188 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"ESSENTIAL",
191 auto &bc_map = bc_mng->getBcMapByBlockName();
192 boundaryMarker = boost::make_shared<std::vector<char unsigned>>();
193 for (
auto b : bc_map) {
194 if (std::regex_match(b.first, std::regex(
"(.*)ESSENTIAL(.*)"))) {
196 for (
int i = 0;
i != b.second->bcMarkers.size(); ++
i) {
197 (*boundaryMarker)[
i] |= b.second->bcMarkers[
i];
208 auto add_domain_lhs_ops = [&](
auto &pipeline) {
212 "U",
"U", [](
double,
double,
double) ->
double {
return k; }));
216 return c * fe_domain_lhs->ts_a;
222 auto add_domain_rhs_ops = [&](
auto &pipeline) {
225 auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
226 auto dot_u_at_gauss_pts = boost::make_shared<VectorDouble>();
228 "U", grad_u_at_gauss_pts));
232 "U", grad_u_at_gauss_pts,
233 [](
double,
double,
double) ->
double {
return k; }));
235 "U", dot_u_at_gauss_pts,
236 [](
const double,
const double,
const double) {
return c; }));
237 auto source_term = [&](
const double x,
const double y,
const double z) {
240 const auto t = fe_domain_lhs->ts_t;
241 return 1e1 * pow(M_E, -M_PI * M_PI *
t) * sin(1. * M_PI * x) *
248 auto add_boundary_lhs_ops = [&](
auto &pipeline) {
251 "U",
"U", [](
const double,
const double,
const double) {
return c; }));
254 auto add_boundary_rhs_ops = [&](
auto &pipeline) {
256 auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
257 auto boundary_function = [&](
const double x,
const double y,
261 const auto t = fe_rhs->ts_t;
269 [](
const double,
const double,
const double) {
return c; }));
275 add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
276 add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
277 add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
278 add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
284 static PetscErrorCode
set(TS ts, PetscReal
t, Vec u, Vec u_t, PetscReal a,
285 Mat A, Mat
B,
void *ctx) {
306 auto create_post_process_element = [&]() {
307 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
308 auto det_ptr = boost::make_shared<VectorDouble>();
309 auto jac_ptr = boost::make_shared<MatrixDouble>();
310 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
312 post_proc_fe->getOpPtrVector().push_back(
314 post_proc_fe->getOpPtrVector().push_back(
317 auto u_ptr = boost::make_shared<VectorDouble>();
318 post_proc_fe->getOpPtrVector().push_back(
323 post_proc_fe->getOpPtrVector().push_back(
325 new OpPPMap(post_proc_fe->getPostProcMesh(),
326 post_proc_fe->getMapGaussPts(),
339 auto set_time_monitor = [&](
auto dm,
auto solver) {
341 boost::shared_ptr<Monitor> monitor_ptr(
342 new Monitor(dm, create_post_process_element()));
343 boost::shared_ptr<ForcesAndSourcesCore> null;
345 monitor_ptr, null, null);
349 auto set_fieldsplit_preconditioner = [&](
auto solver) {
352 CHKERR TSGetSNES(solver, &snes);
354 CHKERR SNESGetKSP(snes, &ksp);
356 CHKERR KSPGetPC(ksp, &pc);
357 PetscBool is_pcfs = PETSC_FALSE;
358 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
360 if (is_pcfs == PETSC_TRUE) {
361 auto bc_mng = mField.getInterface<
BcManager>();
362 auto name_prb =
simple->getProblemName();
363 auto is_all_bc = bc_mng->getBlockIS(name_prb,
"ESSENTIAL",
"U", 0, 0);
365 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
367 <<
"Field split block size " << is_all_bc_size;
368 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
391 auto set_user_ts_jacobian = [&](
auto dm) {
397 auto dm =
simple->getDM();
401 auto solver = pipeline_mng->createTSIM(
404 CHKERR set_user_ts_jacobian(dm);
405 CHKERR set_time_monitor(dm, solver);
406 CHKERR TSSetSolution(solver,
D);
407 CHKERR TSSetFromOptions(solver);
408 CHKERR set_fieldsplit_preconditioner(solver);
439int main(
int argc,
char *argv[]) {
446 auto core_log = logging::core::get();
455 DMType dm_name =
"DMMOFEM";
459 moab::Core mb_instance;
460 moab::Interface &moab = mb_instance;
void simple(double P1[], double P2[], double P3[], double c[], const int N)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#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 ...
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.
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.
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
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
constexpr int SPACE_DIM
[Define dimension]
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
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 TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobina in TS solver.
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr double t
plate stiffness
static constexpr int approx_order
static PetscErrorCode set(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
MoFEM::Interface & mField
MoFEMErrorCode solveSystem()
MoFEMErrorCode outputResults()
HeatEquation(MoFEM::Interface &m_field)
MoFEMErrorCode setupProblem()
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode boundaryCondition()
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode runProgram()
MoFEMErrorCode initialCondition()
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
virtual moab::Interface & get_moab()=0
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)
structure for User Loop Methods on finite elements
default operator for TRI element
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.
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 value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
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]
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc)
static constexpr int saveEveryNthStep
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< PostProcEle > postProc
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode preProcess()
function is run at the beginning of loop