31using namespace boost::multi_index;
32using namespace boost::multiprecision;
33using namespace boost::numeric;
37#ifdef ENABLE_PYTHON_BINDING
38 #include <boost/python.hpp>
39 #include <boost/python/def.hpp>
40 #include <boost/python/numpy.hpp>
41namespace bp = boost::python;
42namespace np = boost::python::numpy;
48static char help[] =
"...\n\n";
50int main(
int argc,
char *argv[]) {
53 const char param_file[] =
"param_file.petsc";
57 auto core_log = logging::core::get();
66 #ifdef ENABLE_PYTHON_BINDING
69 MOFEM_LOG(
"EP", Sev::inform) <<
"Python initialised";
71 MOFEM_LOG(
"EP", Sev::inform) <<
"Python NOT initialised";
82 PetscBool flg = PETSC_TRUE;
88 char restart_file[255];
89 PetscBool restart_flg = PETSC_TRUE;
96 DMType dm_name =
"DMMOFEM";
98 DMType dm_name_mg =
"DMMOFEM_MG";
102 moab::Core moab_core;
103 moab::Interface &moab = moab_core;
105 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
107 pcomm =
new ParallelComm(&moab, PETSC_COMM_WORLD);
109 PetscBool fully_distributed = PETSC_FALSE;
111 &fully_distributed, PETSC_NULLPTR);
112 if (fully_distributed) {
114 if (pcomm->proc_config().proc_size() == 1)
117 option =
"PARALLEL=READ_PART;"
118 "PARALLEL_RESOLVE_SHARED_ENTS;"
119 "PARTITION=PARALLEL_PARTITION";
127 MOFEM_LOG(
"EP", Sev::inform) <<
"Initialise MoFEM database";
130 MOFEM_LOG(
"EP", Sev::inform) <<
"Initialise MoFEM database <- done";
147 auto get_adj = [&](
Range ents,
int dim) {
149 CHKERR moab.get_adjacencies(ents, dim,
false, adj,
150 moab::Interface::UNION);
194 const char *list_materials[LastMaterial] = {
195 "stvenant_kirchhoff",
"mooney_rivlin",
"hencky",
"neo_hookean"};
196 PetscInt choice_material = MooneyRivlin;
198 LastMaterial, &choice_material, PETSC_NULLPTR);
200 switch (choice_material) {
201 case StVenantKirchhoff:
202 MOFEM_LOG(
"EP", Sev::inform) <<
"StVenantKirchhoff material model";
207 MOFEM_LOG(
"EP", Sev::inform) <<
"MooneyRivlin material model";
208 MOFEM_LOG(
"EP", Sev::inform) <<
"MU(1, 0)/2 = " <<
MU(1, 0) / 2;
214 MOFEM_LOG(
"EP", Sev::inform) <<
"Hencky material model";
218 MOFEM_LOG(
"EP", Sev::inform) <<
"Neo-Hookean material model";
232 CHKERR VecSetOption(f_elastic, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
234 PetscInt start_step = 0;
237 CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, restart_file,
238 FILE_MODE_READ, &viewer);
239 CHKERR VecLoad(x_elastic, viewer);
240 CHKERR PetscViewerDestroy(&viewer);
244 std::string restart_file_str(restart_file);
245 std::regex restart_pattern(R
"(restart_([1-9]\d*)\.dat)");
247 if (std::regex_match(restart_file_str, match, restart_pattern)) {
248 start_step = std::stoi(match[1]);
251 "Restart file name must be in the format restart_##.dat");
255 auto ts_elastic =
createTS(PETSC_COMM_WORLD);
256 CHKERR TSSetType(ts_elastic, TSBEULER);
259 CHKERR TSGetAdapt(ts_elastic, &adapt);
260 CHKERR TSAdaptSetType(adapt, TSADAPTNONE);
261 CHKERR TSSetTime(ts_elastic, time);
262 CHKERR TSSetStepNumber(ts_elastic, start_step);
276#ifdef ENABLE_PYTHON_BINDING
277 if (Py_FinalizeEx() < 0) {
Eshelbian plasticity interface.
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
@ MOFEM_DATA_INCONSISTENCY
#define CHKERR
Inline error check.
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.
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
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.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createTS(MPI_Comm comm)
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
PetscErrorCode TSAdaptCreateMoFEM(TSAdapt adapt)
Craete MOFEM adapt.
MoFEMErrorCode setElasticElementOps(const int tag)
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0)
MoFEMErrorCode addFields(const EntityHandle meshset=0)
MoFEMErrorCode createCrackSurfaceMeshset()
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
MoFEMErrorCode getSpatialRotationBc()
MoFEMErrorCode addMaterial_Hencky(double E, double nu)
MoFEMErrorCode projectGeometry(const EntityHandle meshset=0)
MoFEMErrorCode setBlockTagsOnSkin()
MoFEMErrorCode addMaterial_HMHHStVenantKirchhoff(const int tape, const double lambda, const double mu, const double sigma_y)
MoFEMErrorCode addMaterial_HMHMooneyRivlin(const int tape, const double alpha, const double beta, const double lambda, const double sigma_y)
static PetscBool dynamicRelaxation
MoFEMErrorCode solveElastic(TS ts, Vec x)
MoFEMErrorCode setElasticElementToTs(DM dm)
MoFEMErrorCode solveDynamicRelaxation(TS ts, Vec x)
MoFEMErrorCode getSpatialTractionFreeBc(const EntityHandle meshset=0)
MoFEMErrorCode getExternalStrain()
MoFEMErrorCode getSpatialTractionBc()
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0), const EntityHandle meshset=0)
MoFEMErrorCode getSpatialDispBc()
MoFEMErrorCode createExchangeVectors(Sev sev)
MoFEMErrorCode addMaterial_HMHNeohookean(const int tape, const double c10, const double K)
SmartPetscObj< DM > dmElastic
Elastic problem.
static Range getPartEntities(moab::Interface &moab, int part)
static MoFEMErrorCode loadFileRootProcAllRestDistributed(moab::Interface &moab, const char *file_name, int dim, LoadFileFun proc_skin_fun=defaultProcSkinFun, const char *options="PARALLEL=BCAST;PARTITION=")
Root proc has whole mesh, other procs only part of it.
virtual MPI_Comm & get_comm() const =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.
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.
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.
Interface for managing meshsets containing materials and boundary conditions.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.