Implementation of mix-element for Large strains.
#define SINGULARITY
using namespace boost::multi_index;
using namespace boost::multiprecision;
using namespace boost::numeric;
#include <cmath>
#ifdef ENABLE_PYTHON_BINDING
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#endif
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
auto core_log = logging::core::get();
core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "EP"));
LogManager::setLog("EP");
core_log->add_sink(
LogManager::createSink(LogManager::getStrmSelf(), "EPSELF"));
LogManager::setLog("EPSELF");
#ifdef ENABLE_PYTHON_BINDING
Py_Initialize();
np::initialize();
MOFEM_LOG(
"EP", Sev::inform) <<
"Python initialised";
#else
MOFEM_LOG(
"EP", Sev::inform) <<
"Python NOT initialised";
#endif
core_log->add_sink(
LogManager::createSink(LogManager::getStrmSync(), "EPSYNC"));
LogManager::setLog("EPSYNC");
try {
PetscBool flg = PETSC_TRUE;
255, &flg);
255, &flg);
char restart_file[255];
PetscBool restart_flg = PETSC_TRUE;
CHKERR PetscOptionsGetString(PETSC_NULLPTR,
"",
"-restart", restart_file, 255,
&restart_flg);
double time = 0;
CHKERR PetscOptionsGetScalar(PETSC_NULLPTR,
"",
"-time", &time, PETSC_NULLPTR);
DMType dm_name = "DMMOFEM";
CHKERR DMRegister_MoFEM(dm_name);
DMType dm_name_mg = "DMMOFEM_MG";
CHKERR DMRegister_MGViaApproxOrders(dm_name_mg);
moab::Core moab_core;
moab::Interface &moab = moab_core;
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
if (pcomm == NULL)
pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
PetscBool fully_distributed = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULLPTR,
"",
"-fully_distributed",
&fully_distributed, PETSC_NULLPTR);
if (fully_distributed) {
const char *option;
if (pcomm->proc_config().proc_size() == 1)
option = "";
else
option = "PARALLEL=READ_PART;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION";
} else {
CHKERR CommInterface::loadFileRootProcAllRestDistributed(
}
MOFEM_LOG(
"EP", Sev::inform) <<
"Initialise MoFEM database";
MOFEM_LOG(
"EP", Sev::inform) <<
"Initialise MoFEM database <- done";
CHKERR DMRegister_MoFEM(
"DMMOFEM");
BitRefLevel bit_level0 = BitRefLevel().set(0);
0, 3, bit_level0);
auto meshset_ptr = get_temp_meshset_ptr(moab);
*meshset_ptr, CommInterface::getPartEntities(moab, pcomm->rank()));
auto get_adj = [&](
Range ents,
int dim) {
CHKERR moab.get_adjacencies(ents, dim,
false, adj,
moab::Interface::UNION);
return adj;
};
*meshset_ptr,
get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
2));
*meshset_ptr,
get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
1));
*meshset_ptr,
get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
0));
enum MaterialModel {
StVenantKirchhoff,
MooneyRivlin,
Hencky,
Neohookean,
LastMaterial
};
const char *list_materials[LastMaterial] = {
"stvenant_kirchhoff", "mooney_rivlin", "hencky", "neo_hookean"};
PetscInt choice_material = MooneyRivlin;
CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL,
"-material", list_materials,
LastMaterial, &choice_material, PETSC_NULLPTR);
switch (choice_material) {
case StVenantKirchhoff:
MOFEM_LOG(
"EP", Sev::inform) <<
"StVenantKirchhoff material model";
break;
case MooneyRivlin:
MOFEM_LOG(
"EP", Sev::inform) <<
"MooneyRivlin material model";
MOFEM_LOG(
"EP", Sev::inform) <<
"MU(1, 0)/2 = " <<
MU(1, 0) / 2;
break;
case Hencky:
MOFEM_LOG(
"EP", Sev::inform) <<
"Hencky material model";
break;
case Neohookean:
MOFEM_LOG(
"EP", Sev::inform) <<
"Neo-Hookean material model";
break;
default:
break;
}
f_elastic = vectorDuplicate(x_elastic);
CHKERR VecSetOption(f_elastic, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
PetscInt start_step = 0;
if (restart_flg) {
PetscViewer viewer;
CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, restart_file,
FILE_MODE_READ, &viewer);
CHKERR VecLoad(x_elastic, viewer);
CHKERR PetscViewerDestroy(&viewer);
SCATTER_REVERSE);
std::string restart_file_str(restart_file);
std::regex restart_pattern(R"(restart_([1-9]\d*)\.dat)");
std::smatch match;
if (std::regex_match(restart_file_str, match, restart_pattern)) {
start_step = std::stoi(match[1]);
} else {
"Restart file name must be in the format restart_##.dat");
}
}
auto ts_elastic = createTS(PETSC_COMM_WORLD);
CHKERR TSSetType(ts_elastic, TSBEULER);
TSAdapt adapt;
CHKERR TSGetAdapt(ts_elastic, &adapt);
CHKERR TSAdaptSetType(adapt, TSADAPTNONE);
CHKERR TSSetTime(ts_elastic, time);
CHKERR TSSetStepNumber(ts_elastic, start_step);
} else {
}
}
#ifdef ENABLE_PYTHON_BINDING
if (Py_FinalizeEx() < 0) {
exit(120);
}
#endif
}
Eshelbian plasticity interface.
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
@ MOFEM_DATA_INCONSISTENCY
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
implementation of Data Operators for Forces and Sources
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.
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.
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.