32#include <BasicFiniteElements.hpp>
36static char help[] =
"...\n\n";
70constexpr double c = 1;
71constexpr double k = 1;
95 "out_level_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
102 boost::shared_ptr<PostProcEle>
postProc;
184 1,
BLOCKSET, 2, inner_surface,
true);
185 if (!inner_surface.empty()) {
186 Range inner_surface_verts;
188 inner_surface_verts,
false);
190 init_u, MBVERTEX, inner_surface_verts,
"U");
202 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"ESSENTIAL",
205 auto &bc_map = bc_mng->getBcMapByBlockName();
206 boundaryMarker = boost::make_shared<std::vector<char unsigned>>();
207 for (
auto b : bc_map) {
208 if (std::regex_match(
b.first, std::regex(
"(.*)ESSENTIAL(.*)"))) {
210 for (
int i = 0;
i !=
b.second->bcMarkers.size(); ++
i) {
211 (*boundaryMarker)[
i] |=
b.second->bcMarkers[
i];
222 auto add_domain_base_ops = [&](
auto &pipeline) {
223 auto det_ptr = boost::make_shared<VectorDouble>();
224 auto jac_ptr = boost::make_shared<MatrixDouble>();
225 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
232 auto add_domain_lhs_ops = [&](
auto &pipeline) {
235 "U",
"U", [](
double,
double,
double) ->
double {
return k; }));
239 return c * fe_domain_lhs->ts_a;
245 auto add_domain_rhs_ops = [&](
auto &pipeline) {
247 auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
248 auto dot_u_at_gauss_pts = boost::make_shared<VectorDouble>();
250 "U", grad_u_at_gauss_pts));
254 "U", grad_u_at_gauss_pts,
255 [](
double,
double,
double) ->
double {
return k; }));
257 "U", dot_u_at_gauss_pts,
258 [](
const double,
const double,
const double) {
return c; }));
259 auto source_term = [&](
const double x,
const double y,
const double z) {
262 const auto t = fe_domain_lhs->ts_t;
263 return 1e1 * pow(M_E, -M_PI * M_PI *
t) * sin(1. * M_PI * x) *
270 auto add_boundary_base_ops = [&](
auto &pipeline) {};
272 auto add_lhs_base_ops = [&](
auto &pipeline) {
275 "U",
"U", [](
const double,
const double,
const double) {
return c; }));
278 auto add_rhs_base_ops = [&](
auto &pipeline) {
280 auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
281 auto boundary_function = [&](
const double x,
const double y,
285 const auto t = fe_rhs->ts_t;
293 [](
const double,
const double,
const double) {
return c; }));
299 add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
300 add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
301 add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
302 add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
304 add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
305 add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
306 add_lhs_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
307 add_rhs_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
318 auto create_post_process_element = [&]() {
319 auto post_froc_fe = boost::make_shared<PostProcEle>(
mField);
320 post_froc_fe->generateReferenceElementMesh();
321 auto det_ptr = boost::make_shared<VectorDouble>();
322 auto jac_ptr = boost::make_shared<MatrixDouble>();
323 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
325 post_froc_fe->getOpPtrVector().push_back(
327 post_froc_fe->getOpPtrVector().push_back(
329 post_froc_fe->addFieldValuesPostProc(
"U");
333 auto set_time_monitor = [&](
auto dm,
auto solver) {
335 boost::shared_ptr<Monitor> monitor_ptr(
336 new Monitor(dm, create_post_process_element()));
337 boost::shared_ptr<ForcesAndSourcesCore> null;
339 monitor_ptr, null, null);
343 auto set_fieldsplit_preconditioner = [&](
auto solver) {
346 CHKERR TSGetSNES(solver, &snes);
348 CHKERR SNESGetKSP(snes, &ksp);
350 CHKERR KSPGetPC(ksp, &pc);
351 PetscBool is_pcfs = PETSC_FALSE;
352 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
354 if (is_pcfs == PETSC_TRUE) {
356 auto name_prb =
simple->getProblemName();
357 auto is_all_bc = bc_mng->getBlockIS(name_prb,
"ESSENTIAL",
"U", 0, 0);
359 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
361 <<
"Field split block size " << is_all_bc_size;
362 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
368 auto dm =
simple->getDM();
372 auto solver = pipeline_mng->createTSIM();
373 CHKERR TSSetSolution(solver,
D);
374 CHKERR set_time_monitor(dm, solver);
375 CHKERR TSSetSolution(solver,
D);
376 CHKERR TSSetFromOptions(solver);
377 CHKERR set_fieldsplit_preconditioner(solver);
379 CHKERR TSSolve(solver, NULL);
381 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
382 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
411int main(
int argc,
char *argv[]) {
418 auto core_log = logging::core::get();
420 LogManager::createSink(LogManager::getStrmWorld(),
"EXAMPLE"));
421 LogManager::setLog(
"EXAMPLE");
427 DMType dm_name =
"DMMOFEM";
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
MoFEM::FaceElementForcesAndSourcesCore FaceEle
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
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
#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 smartCreateDMVector(DM dm)
Get smart vector from DM.
#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
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
[Define dimension]
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBoundaryTimeScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 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, 2 > OpDomainGradGrad
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 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 PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
const double D
diffusivity
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBoundaryTimeScalarField
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
constexpr double t
plate stiffness
static constexpr int approx_order
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()
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
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.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Simple interface for fast problem set-up.
PetscInt ts_step
time step
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
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