v0.15.5
Loading...
Searching...
No Matches
Macros | Functions | Variables
ep.cpp File Reference
#include <MoFEM.hpp>
#include <cmath>
#include <cholesky.hpp>
#include <EshelbianPlasticity.hpp>

Go to the source code of this file.

Macros

#define SINGULARITY
 

Functions

int main (int argc, char *argv[])
 

Variables

constexpr int adolc_tag = 1
 
static char help [] = "...\n\n"
 

Macro Definition Documentation

◆ SINGULARITY

#define SINGULARITY

Definition at line 27 of file ep.cpp.

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 50 of file ep.cpp.

50 {
51
52 // initialize petsc
53 const char param_file[] = "param_file.petsc";
54 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
55
56 // Add logging channel for example
57 auto core_log = logging::core::get();
58 core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "EP"));
59 LogManager::setLog("EP");
60 MOFEM_LOG_TAG("EP", "ep");
61 core_log->add_sink(
62 LogManager::createSink(LogManager::getStrmSelf(), "EPSELF"));
63 LogManager::setLog("EPSELF");
64 MOFEM_LOG_TAG("EPSELF", "ep");
65
66 #ifdef ENABLE_PYTHON_BINDING
67 Py_Initialize();
68 np::initialize();
69 MOFEM_LOG("EP", Sev::inform) << "Python initialised";
70#else
71 MOFEM_LOG("EP", Sev::inform) << "Python NOT initialised";
72#endif
73
74 core_log->add_sink(
75 LogManager::createSink(LogManager::getStrmSync(), "EPSYNC"));
76 LogManager::setLog("EPSYNC");
77 MOFEM_LOG_TAG("EPSYNC", "ep");
78
79 try {
80
81 // Get mesh file
82 PetscBool flg = PETSC_TRUE;
83 char mesh_file_name[255];
84 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-my_file", mesh_file_name,
85 255, &flg);
86 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-file_name", mesh_file_name,
87 255, &flg);
88 char restart_file[255] = "no_restart_file";
89 PetscBool restart_flg = PETSC_TRUE;
90 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-restart", restart_file, 255,
91 &restart_flg);
92 double time = 0;
93 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-time", &time, PETSC_NULLPTR);
94
95
96 MOFEM_LOG("EP", Sev::inform) << "Mesh file: " << mesh_file_name;
97 MOFEM_LOG("EP", Sev::inform) << "Restart file: " << restart_file;
98 MOFEM_LOG("EP", Sev::inform) << "Time: " << time;
99
100 // Register DM Manager
101 DMType dm_name = "DMMOFEM";
102 CHKERR DMRegister_MoFEM(dm_name);
103 DMType dm_name_mg = "DMMOFEM_MG";
104 CHKERR DMRegister_MGViaApproxOrders(dm_name_mg);
105
106 // Create MoAB database
107 moab::Core moab_core;
108 moab::Interface &moab = moab_core;
109
110 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
111 if (pcomm == NULL)
112 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
113 // Read mesh to MOAB
114 PetscBool fully_distributed = PETSC_FALSE;
115 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-fully_distributed",
116 &fully_distributed, PETSC_NULLPTR);
117 if (fully_distributed) {
118 const char *option;
119 if (pcomm->proc_config().proc_size() == 1)
120 option = "";
121 else
122 option = "PARALLEL=READ_PART;"
123 "PARALLEL_RESOLVE_SHARED_ENTS;"
124 "PARTITION=PARALLEL_PARTITION";
125 CHKERR moab.load_file(mesh_file_name, 0, option);
126 } else {
127 CHKERR CommInterface::loadFileRootProcAllRestDistributed(
128 moab, mesh_file_name, SPACE_DIM);
129 }
130
131 // Create MoFEM database and link it to MoAB
132 MOFEM_LOG("EP", Sev::inform) << "Initialise MoFEM database";
133 MoFEM::Core mofem_core(moab);
134 MoFEM::Interface &m_field = mofem_core;
135 MOFEM_LOG("EP", Sev::inform) << "Initialise MoFEM database <- done";
136
137 CHKERR m_field.getInterface<MeshsetsManager>()->setMeshsetFromFile();
138
139 // Register mofem DM
140 CHKERR DMRegister_MoFEM("DMMOFEM");
141
142 BitRefLevel bit_level0 = BitRefLevel().set(0);
143 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
144 0, 3, bit_level0);
145
146 // Data stuctures
147 EshelbianCore ep(m_field);
148
149 auto meshset_ptr = get_temp_meshset_ptr(moab);
150 CHKERR moab.add_entities(
151 *meshset_ptr, CommInterface::getPartEntities(moab, pcomm->rank()));
152 auto get_adj = [&](Range ents, int dim) {
153 Range adj;
154 CHKERR moab.get_adjacencies(ents, dim, false, adj,
155 moab::Interface::UNION);
156 return adj;
157 };
158 CHKERR moab.add_entities(
159 *meshset_ptr,
160 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
161 .subset_by_dimension(SPACE_DIM),
162 2));
163 CHKERR moab.add_entities(
164 *meshset_ptr,
165 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
166 .subset_by_dimension(SPACE_DIM),
167 1));
168 CHKERR moab.add_entities(
169 *meshset_ptr,
170 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
171 .subset_by_dimension(SPACE_DIM),
172 0));
173
174 // create growing crack surface meshset
175 CHKERR ep.createCrackSurfaceMeshset();
176
177 CHKERR ep.getSpatialDispBc();
178 CHKERR ep.getSpatialRotationBc();
179 CHKERR ep.getSpatialTractionBc();
180 CHKERR ep.getSpatialTractionFreeBc();
181 CHKERR ep.getExternalStrain();
182 CHKERR ep.setBlockTagsOnSkin();
183
184 CHKERR ep.createExchangeVectors(Sev::inform);
185
186 CHKERR ep.projectInternalStress();
187 CHKERR ep.addFields(*meshset_ptr);
188 CHKERR ep.projectGeometry(*meshset_ptr);
189 CHKERR ep.addVolumeFiniteElement(*meshset_ptr);
190 CHKERR ep.addBoundaryFiniteElement(*meshset_ptr);
191 CHKERR ep.addDMs();
192
193 const char *list_materials[EshelbianCore::MaterialModel::LastMaterial] = {
194 "stvenant_kirchhoff", "mooney_rivlin", "hencky", "neo_hookean"};
195 PetscInt choice_material = EshelbianCore::materialModel;
196 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-material",
197 list_materials,
198 EshelbianCore::MaterialModel::LastMaterial,
199 &choice_material, PETSC_NULLPTR);
200 EshelbianCore::materialModel =
201 static_cast<EshelbianCore::MaterialModel>(choice_material);
202
203 switch (choice_material) {
204 case EshelbianCore::MaterialModel::StVenantKirchhoff:
205 MOFEM_LOG("EP", Sev::inform) << "StVenantKirchhoff material model";
206 CHKERR ep.addMaterial_HMHHStVenantKirchhoff(adolc_tag, LAMBDA(1, 0.25),
207 MU(1, 0.25), 0);
208 break;
209 case EshelbianCore::MaterialModel::MooneyRivlin:
210 MOFEM_LOG("EP", Sev::inform) << "MooneyRivlin material model";
211 MOFEM_LOG("EP", Sev::inform) << "MU(1, 0)/2 = " << MU(1, 0) / 2;
212 MOFEM_LOG("EP", Sev::inform) << "LAMBDA(1, 0) = " << LAMBDA(1, 0);
213 CHKERR ep.addMaterial_HMHMooneyRivlin(adolc_tag, MU(1, 0) / 2., 0,
214 LAMBDA(1, 0), 0);
215 break;
216 case EshelbianCore::MaterialModel::Hencky:
217 MOFEM_LOG("EP", Sev::inform) << "Hencky material model";
218 CHKERR ep.addMaterial_Hencky(5., 0.25);
219 break;
220 case EshelbianCore::MaterialModel::Neohookean:
221 MOFEM_LOG("EP", Sev::inform) << "Neo-Hookean material model";
222 CHKERR ep.addMaterial_HMHNeohookean(adolc_tag, 1.0, 1.0);
223 break;
224 default:
225 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown material");
226 break;
227 }
228
229 CHKERR ep.setElasticElementOps(adolc_tag);
230 CHKERR ep.setElasticElementToTs(ep.dmElastic);
231
232 auto x_elastic = createDMVector(ep.dmElastic);
233 auto f_elastic = vectorDuplicate(x_elastic);
234 CHKERR VecSetOption(f_elastic, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
235
236 PetscInt start_step = 0;
237 if (restart_flg) {
238 PetscViewer viewer;
239 CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, restart_file,
240 FILE_MODE_READ, &viewer);
241 CHKERR VecLoad(x_elastic, viewer);
242 CHKERR PetscViewerDestroy(&viewer);
243 CHKERR VecGhostUpdateBegin(x_elastic, INSERT_VALUES, SCATTER_FORWARD);
244 CHKERR VecGhostUpdateEnd(x_elastic, INSERT_VALUES, SCATTER_FORWARD);
245
246 std::string restart_file_str(restart_file);
247 std::regex restart_pattern(R"(restart_([1-9]\d*)\.dat)");
248 std::smatch match;
249 if (std::regex_search(restart_file_str, match, restart_pattern)) {
250 start_step = std::stoi(match[1]);
251 } else {
252 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
253 "Restart file name must be in the format restart_##.dat");
254 }
255 }
256
257 auto ts_elastic = createTS(PETSC_COMM_WORLD);
258 CHKERR TSSetType(ts_elastic, TSBEULER);
259 CHKERR TSAdaptRegister(TSADAPTMOFEM, TSAdaptCreateMoFEM);
260 TSAdapt adapt;
261 CHKERR TSGetAdapt(ts_elastic, &adapt);
262 CHKERR TSAdaptSetType(adapt, TSADAPTNONE);
263
264 switch (EshelbianCore::solverType) {
265 case EshelbianCore::SolverType::TimeSolver:
266 MOFEM_LOG("EP", Sev::inform) << "Solver type: TS";
267 CHKERR TSSetTime(ts_elastic, time);
268 CHKERR TSSetStepNumber(ts_elastic, start_step);
269 CHKERR ep.solveElastic(ts_elastic, x_elastic);
270 break;
271 case EshelbianCore::SolverType::DynamicRelaxation:
272 MOFEM_LOG("EP", Sev::inform) << "Solver type: Dynamic Relaxation";
273 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
274 SCATTER_REVERSE);
275 CHKERR ep.solveDynamicRelaxation(ts_elastic, x_elastic, start_step, time);
276 break;
277 case EshelbianCore::SolverType::Cohesive:
278 MOFEM_LOG("EP", Sev::inform) << "Solver type: Cohesive";
279 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
280 SCATTER_REVERSE);
281 CHKERR ep.solveCohesiveCrackGrowth(ts_elastic, x_elastic, start_step,
282 time);
283 break;
284 default:
285 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
286 "Unknown solver type");
287 break;
288 }
289
290 }
291 CATCH_ERRORS;
292
293 // finish work cleaning memory, getting statistics, etc.
294 MoFEM::Core::Finalize();
295
296#ifdef ENABLE_PYTHON_BINDING
297 if (Py_FinalizeEx() < 0) {
298 exit(120);
299 }
300#endif
301}
TSADAPTMOFEM
#define TSADAPTMOFEM
Definition TsCtx.hpp:10
SPACE_DIM
constexpr int SPACE_DIM
Definition child_and_parent.cpp:17
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
help
static char help[]
Definition ep.cpp:48
adolc_tag
constexpr int adolc_tag
Definition ep.cpp:25
MU
#define MU(E, NU)
Definition fem_tools.h:23
LAMBDA
#define LAMBDA(E, NU)
Definition fem_tools.h:22
MoFEM::DMoFEMMeshToLocalVector
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:514
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEM::DMRegister_MGViaApproxOrders
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
Definition PCMGSetUpViaApproxOrders.cpp:302
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1234
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition LogManager.cpp:392
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition LogManager.hpp:347
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition DeprecatedPetsc.hpp:182
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition DeprecatedPetsc.hpp:162
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition PetscSmartObj.hpp:223
MoFEM::createTS
auto createTS(MPI_Comm comm)
Definition PetscSmartObj.hpp:257
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition DeprecatedPetsc.hpp:203
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition DeprecatedPetsc.hpp:172
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition Templates.hpp:1980
MoFEM::TSAdaptCreateMoFEM
PetscErrorCode TSAdaptCreateMoFEM(TSAdapt adapt)
Craete MOFEM adapt.
Definition TsCtx.cpp:829
EshelbianCore::solverType
static enum SolverType solverType
Definition EshelbianCore.hpp:24
EshelbianCore::materialModel
static enum MaterialModel materialModel
Definition EshelbianCore.hpp:25
MoFEM::BitRefManager
Managing BitRefLevels.
Definition BitRefManager.hpp:21
MoFEM::CommInterface::getPartEntities
static Range getPartEntities(moab::Interface &moab, int part)
Definition CommInterface.cpp:1419
MoFEM::CommInterface::loadFileRootProcAllRestDistributed
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.
Definition CommInterface.cpp:1206
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:82
MoFEM::CoreTmp< 0 >::Initialize
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
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition LogManager.cpp:293
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition LogManager.cpp:339
MoFEM::LogManager::getStrmSync
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition LogManager.cpp:343
MoFEM::LogManager::getStrmSelf
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.
Definition LogManager.cpp:335
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition MeshsetsManager.hpp:208
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
restart_file
char restart_file[255]
Definition thermoplastic.cpp:190
restart_flg
PetscBool restart_flg
Definition thermoplastic.cpp:191

Variable Documentation

◆ adolc_tag

constexpr int adolc_tag = 1
constexpr

Definition at line 25 of file ep.cpp.

◆ help

char help[] = "...\n\n"
static

Definition at line 48 of file ep.cpp.

Generated by Doxygen 1.9.8 and hosted at University of Glasgow