v0.16.3
Loading...
Searching...
No Matches
Macros | Functions | Variables
ep.cpp File Reference
#include <MoFEM.hpp>
#include <EshelbianRestart.hpp>
#include <cmath>
#include <cholesky.hpp>
#include <EshelbianPlasticity.hpp>
#include <MatElastic.hpp>
#include <MatOps.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 25 of file ep.cpp.

Function Documentation

◆ main()

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

Definition at line 51 of file ep.cpp.

51 {
52
53 // initialize petsc
54 const char param_file[] = "param_file.petsc";
55 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
56
57 // Add logging channel for example
58 auto core_log = logging::core::get();
59 core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "EP"));
60 LogManager::setLog("EP");
61 MOFEM_LOG_TAG("EP", "ep");
62 core_log->add_sink(
63 LogManager::createSink(LogManager::getStrmSelf(), "EPSELF"));
64 LogManager::setLog("EPSELF");
65 MOFEM_LOG_TAG("EPSELF", "ep");
66
67#ifdef ENABLE_PYTHON_BINDING
68 Py_Initialize();
69 np::initialize();
70 MOFEM_LOG("EP", Sev::inform) << "Python initialised";
71#else
72 MOFEM_LOG("EP", Sev::inform) << "Python NOT initialised";
73#endif
74
75 core_log->add_sink(
76 LogManager::createSink(LogManager::getStrmSync(), "EPSYNC"));
77 LogManager::setLog("EPSYNC");
78 MOFEM_LOG_TAG("EPSYNC", "ep");
79
80 try {
81
82 CHKERR JsonConfigManager::setMeshFileFromJson();
83
84 // Get mesh file
85 PetscBool flg = PETSC_FALSE;
86 char mesh_file_name[255] = "";
87 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-my_file", mesh_file_name,
88 255, &flg);
89 if (!flg) {
90 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-file_name",
91 mesh_file_name, 255, &flg);
92 }
93 if (!flg) {
94 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
95 "Missing required mesh option: use -file_name or -my_file");
96 }
97 char restart_file[255] = "no_restart_file";
98 PetscBool restart_flg = PETSC_TRUE;
99 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-restart", restart_file,
100 255, &restart_flg);
101 double time = 0;
102 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-time", &time,
103 PETSC_NULLPTR);
104
105 MOFEM_LOG("EP", Sev::inform) << "Mesh file: " << mesh_file_name;
106 MOFEM_LOG("EP", Sev::inform) << "Restart file: " << restart_file;
107 MOFEM_LOG("EP", Sev::inform) << "Time: " << time;
108
109 // Register DM Manager
110 DMType dm_name = "DMMOFEM";
111 CHKERR DMRegister_MoFEM(dm_name);
112 DMType dm_name_mg = "DMMOFEM_MG";
113 CHKERR DMRegister_MGViaApproxOrders(dm_name_mg);
114
115 // Create MoAB database
116 moab::Core moab_core;
117 moab::Interface &moab = moab_core;
118
119 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
120 if (pcomm == NULL)
121 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
122 // Read mesh to MOAB
123 PetscBool fully_distributed = PETSC_FALSE;
124 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-fully_distributed",
125 &fully_distributed, PETSC_NULLPTR);
126 if (fully_distributed) {
127 const char *option;
128 if (pcomm->proc_config().proc_size() == 1)
129 option = "";
130 else
131 option = "PARALLEL=READ_PART;"
132 "PARALLEL_RESOLVE_SHARED_ENTS;"
133 "PARTITION=PARALLEL_PARTITION";
134 CHKERR moab.load_file(mesh_file_name, 0, option);
135 } else {
136 CHKERR CommInterface::loadFileRootProcAllRestDistributed(
137 moab, mesh_file_name, SPACE_DIM);
138 }
139
140 // Create MoFEM database and link it to MoAB
141 MOFEM_LOG("EP", Sev::inform) << "Initialise MoFEM database";
142 MoFEM::Core mofem_core(moab);
143 MoFEM::Interface &m_field = mofem_core;
144 MOFEM_LOG("EP", Sev::inform) << "Initialise MoFEM database <- done";
145 CHKERR JsonConfigManager::addCanonicalAttributeNames(
146 {"kn", "kt"}); // TODO: remove this line and add kn kt to the core
147 // JsonConfigManagerCanonicalAttributes.hpp
148 CHKERR m_field.getInterface<MeshsetsManager>()->setMeshsetFromFile();
149
150 // Register mofem DM
151 CHKERR DMRegister_MoFEM("DMMOFEM");
152
153 BitRefLevel bit_level0 = BitRefLevel().set(0);
154 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
155 0, 3, bit_level0);
156
157 // Data stuctures
158 EshelbianCore ep(m_field);
159 CHKERR PhysicalEquations::create(ep);
160
161 auto meshset_ptr = get_temp_meshset_ptr(moab);
162 CHKERR moab.add_entities(
163 *meshset_ptr, CommInterface::getPartEntities(moab, pcomm->rank()));
164 auto get_adj = [&](Range ents, int dim) {
165 Range adj;
166 CHKERR moab.get_adjacencies(ents, dim, false, adj,
167 moab::Interface::UNION);
168 return adj;
169 };
170 CHKERR moab.add_entities(
171 *meshset_ptr,
172 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
173 .subset_by_dimension(SPACE_DIM),
174 2));
175 CHKERR moab.add_entities(
176 *meshset_ptr,
177 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
178 .subset_by_dimension(SPACE_DIM),
179 1));
180 CHKERR moab.add_entities(
181 *meshset_ptr,
182 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
183 .subset_by_dimension(SPACE_DIM),
184 0));
185
186 // create growing crack surface meshset
187 CHKERR ep.createCrackSurfaceMeshset();
188
189 CHKERR ep.getSpatialDispBc();
190 CHKERR ep.getSpatialRotationBc();
191 CHKERR ep.getSpatialTractionBc();
192 CHKERR ep.getSpatialTractionFreeBc();
193 CHKERR ep.getExternalStrain();
194 CHKERR ep.setBlockTagsOnSkin();
195
196 CHKERR ep.createExchangeVectors(Sev::inform);
197
198 CHKERR ep.projectMaterialTags();
199 CHKERR ep.resolveDissipationEntities(*meshset_ptr);
200 CHKERR ep.addFields(*meshset_ptr);
201 CHKERR ep.projectGeometry(*meshset_ptr, 0.0);
202 CHKERR ep.addVolumeFiniteElement(*meshset_ptr);
203 CHKERR ep.addBoundaryFiniteElement(*meshset_ptr);
204 CHKERR ep.addDMs();
205
206 CHKERR ep.setElasticElementOps(adolc_tag);
207 CHKERR ep.setElasticElementToTs(ep.dmElastic);
208
209 auto x_elastic = createDMVector(ep.dmElastic);
210
211 PetscInt start_step = 0;
212 if (restart_flg) {
213 if (EshelbianCore::solverType ==
214 EshelbianCore::SolverType::IncrementalOptimization ||
215 EshelbianCore::solverType ==
216 EshelbianCore::SolverType::TestIncrementalOptimizationLayout ||
217 EshelbianCore::solverType ==
218 EshelbianCore::SolverType::
219 TestIncrementalOptimizationTransaction ||
220 EshelbianCore::solverType ==
221 EshelbianCore::SolverType::
222 TestIncrementalOptimizationObjectiveDerivative ||
223 EshelbianCore::solverType ==
224 EshelbianCore::SolverType::
225 TestIncrementalOptimizationConstraintDerivative)
226 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
227 "Incremental optimization does not support -restart");
228 CHKERR validateRestartLayout(ep, x_elastic, restart_file);
229 PetscViewer viewer;
230 CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, restart_file,
231 FILE_MODE_READ, &viewer);
232 CHKERR VecLoad(x_elastic, viewer);
233 CHKERR PetscViewerDestroy(&viewer);
234 CHKERR VecGhostUpdateBegin(x_elastic, INSERT_VALUES, SCATTER_FORWARD);
235 CHKERR VecGhostUpdateEnd(x_elastic, INSERT_VALUES, SCATTER_FORWARD);
236 // Solver setup gathers its initial vector from mesh coefficients.
237 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
238 SCATTER_REVERSE);
239 CHKERR VecCopy(x_elastic, ep.solTSStep);
240 CHKERR VecGhostUpdateBegin(ep.solTSStep, INSERT_VALUES, SCATTER_FORWARD);
241 CHKERR VecGhostUpdateEnd(ep.solTSStep, INSERT_VALUES, SCATTER_FORWARD);
242
243 std::string restart_file_str(restart_file);
244 std::regex restart_pattern(R"(restart_([1-9]\d*)\.dat)");
245 std::smatch match;
246 if (std::regex_search(restart_file_str, match, restart_pattern)) {
247 start_step = std::stoi(match[1]);
248 } else {
249 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
250 "Restart file name must be in the format restart_##.dat");
251 }
252 }
253
254 auto ts_elastic = createTS(PETSC_COMM_WORLD);
255 CHKERR TSSetType(ts_elastic, TSBEULER);
256 CHKERR TSAdaptRegister(TSADAPTMOFEM, TSAdaptCreateMoFEM);
257 TSAdapt adapt;
258 CHKERR TSGetAdapt(ts_elastic, &adapt);
259 CHKERR TSAdaptSetType(adapt, TSADAPTNONE);
260
261 switch (EshelbianCore::solverType) {
262 case EshelbianCore::SolverType::TimeSolver:
263 MOFEM_LOG("EP", Sev::inform) << "Solver type: TS";
264 CHKERR TSSetTime(ts_elastic, time);
265 CHKERR TSSetStepNumber(ts_elastic, start_step);
266 CHKERR ep.solveElastic(ts_elastic, x_elastic);
267 break;
268 case EshelbianCore::SolverType::DynamicRelaxation:
269 MOFEM_LOG("EP", Sev::inform) << "Solver type: Dynamic Relaxation";
270 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
271 SCATTER_REVERSE);
272 CHKERR ep.solveDynamicRelaxation(ts_elastic, x_elastic, start_step, time);
273 break;
274 case EshelbianCore::SolverType::IncrementalOptimization:
275 MOFEM_LOG("EP", Sev::inform) << "Solver type: Incremental Optimization";
276 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
277 SCATTER_REVERSE);
278 CHKERR ep.solveIncrementalOptimizationTAO(ts_elastic, x_elastic,
279 start_step, time);
280 break;
281 case EshelbianCore::SolverType::LoadFactor:
282 MOFEM_LOG("EP", Sev::inform) << "Solver type: Load Factor";
283 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
284 SCATTER_REVERSE);
285 CHKERR ep.solveLoadFactor(ts_elastic, x_elastic, start_step, time);
286 break;
287 case EshelbianCore::SolverType::ShapeOptimization:
288 MOFEM_LOG("EP", Sev::inform) << "Solver type: Shape Optimization";
289 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
290 SCATTER_REVERSE);
291 CHKERR ep.solveSchapeOptimisation(ts_elastic, x_elastic, start_step,
292 time);
293 break;
294 case EshelbianCore::SolverType::TestTopologicalDerivative:
295 MOFEM_LOG("EP", Sev::inform)
296 << "Solver type: Test Topological Derivative";
297 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
298 SCATTER_REVERSE);
299 CHKERR ep.solveTestTopologicalDerivative(ts_elastic, x_elastic,
300 start_step, time);
301 break;
302 case EshelbianCore::SolverType::TestEquilibratedMechanicalValue:
303 MOFEM_LOG("EP", Sev::inform)
304 << "Solver type: Test Equilibrated Mechanical Value";
305 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
306 SCATTER_REVERSE);
307 CHKERR ep.solveTestEquilibratedMechanicalValue(ts_elastic, x_elastic,
308 start_step, time);
309 break;
310 case EshelbianCore::SolverType::TestIncrementalOptimizationLayout:
311 MOFEM_LOG("EP", Sev::inform)
312 << "Solver type: Test Incremental Optimization Layout";
313 CHKERR ep.solveTestIncrementalOptimizationLayout(ts_elastic, x_elastic,
314 start_step, time);
315 break;
316 case EshelbianCore::SolverType::TestIncrementalOptimizationTransaction:
317 MOFEM_LOG("EP", Sev::inform)
318 << "Solver type: Test Incremental Optimization Transaction";
319 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
320 SCATTER_REVERSE);
321 CHKERR ep.solveTestIncrementalOptimizationTransaction(
322 ts_elastic, x_elastic, start_step, time);
323 break;
324 case EshelbianCore::SolverType::
325 TestIncrementalOptimizationObjectiveDerivative:
326 MOFEM_LOG("EP", Sev::inform)
327 << "Solver type: Test Incremental Optimization Objective Derivative";
328 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
329 SCATTER_REVERSE);
330 CHKERR ep.solveTestIncrementalOptimizationObjectiveDerivative(
331 ts_elastic, x_elastic, start_step, time);
332 break;
333 case EshelbianCore::SolverType::
334 TestIncrementalOptimizationConstraintDerivative:
335 MOFEM_LOG("EP", Sev::inform)
336 << "Solver type: Test Incremental Optimization Constraint Derivative";
337 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
338 SCATTER_REVERSE);
339 CHKERR ep.solveTestIncrementalOptimizationConstraintDerivative(
340 ts_elastic, x_elastic, start_step, time);
341 break;
342 default:
343 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
344 "Unknown solver type");
345 break;
346 }
347 }
348 CATCH_ERRORS;
349
350 // finish work cleaning memory, getting statistics, etc.
351 MoFEM::Core::Finalize();
352
353#ifdef ENABLE_PYTHON_BINDING
354 if (Py_FinalizeEx() < 0) {
355 exit(120);
356 }
357#endif
358}
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:49
adolc_tag
constexpr int adolc_tag
Definition ep.cpp:23
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode, RowColData rc=RowColData::COL)
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, RowColData rc=RowColData::COL)
Get smart vector from DM.
Definition DMMoFEM.hpp:1237
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
EshelbianPlasticity::validateRestartLayout
MoFEM::MoFEMErrorCode validateRestartLayout(const EshelbianCore &ep, Vec state, const std::string &binary_file)
Definition EshelbianRestart.cpp:284
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::createTS
auto createTS(MPI_Comm comm)
Definition PetscSmartObj.hpp:257
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:2424
MoFEM::TSAdaptCreateMoFEM
PetscErrorCode TSAdaptCreateMoFEM(TSAdapt adapt)
Craete MOFEM adapt.
Definition TsCtx.cpp:829
EshelbianCore::solverType
static enum SolverType solverType
Definition EshelbianCore.hpp:45
EshelbianCore::TestIncrementalOptimizationConstraintDerivative
@ TestIncrementalOptimizationConstraintDerivative
Definition EshelbianCore.hpp:38
EshelbianCore::TestIncrementalOptimizationLayout
@ TestIncrementalOptimizationLayout
Definition EshelbianCore.hpp:35
EshelbianCore::TestIncrementalOptimizationTransaction
@ TestIncrementalOptimizationTransaction
Definition EshelbianCore.hpp:36
EshelbianCore::TestIncrementalOptimizationObjectiveDerivative
@ TestIncrementalOptimizationObjectiveDerivative
Definition EshelbianCore.hpp:37
EshelbianCore::TestEquilibratedMechanicalValue
@ TestEquilibratedMechanicalValue
Definition EshelbianCore.hpp:34
EshelbianPlasticity::PhysicalEquations::create
static MoFEMErrorCode create(EshelbianCore &ep)
Select and configure the material before registering fields and DMs.
Definition EshelbianMat.cpp:664
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:83
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:68
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:123
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::JsonConfigManager::setMeshFileFromJson
static MoFEMErrorCode setMeshFileFromJson()
Set -file_name from JSON before Core is available.
Definition JsonConfigManager.cpp:451
MoFEM::JsonConfigManager::addCanonicalAttributeNames
static MoFEMErrorCode addCanonicalAttributeNames(const std::vector< std::string > &names)
Definition JsonConfigManager.cpp:419
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:252
restart_flg
PetscBool restart_flg
Definition thermoplastic.cpp:253

Variable Documentation

◆ adolc_tag

constexpr int adolc_tag = 1
constexpr

Definition at line 23 of file ep.cpp.

◆ help

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

Definition at line 49 of file ep.cpp.

Generated by Doxygen 1.9.8 and hosted at University of Glasgow