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

52 {
53
54 // initialize petsc
55 const char param_file[] = "param_file.petsc";
56 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
57
58 // Add logging channel for example
59 auto core_log = logging::core::get();
60 core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "EP"));
61 LogManager::setLog("EP");
62 MOFEM_LOG_TAG("EP", "ep");
63 core_log->add_sink(
64 LogManager::createSink(LogManager::getStrmSelf(), "EPSELF"));
65 LogManager::setLog("EPSELF");
66 MOFEM_LOG_TAG("EPSELF", "ep");
67
68 #ifdef ENABLE_PYTHON_BINDING
69 Py_Initialize();
70 np::initialize();
71 MOFEM_LOG("EP", Sev::inform) << "Python initialised";
72#else
73 MOFEM_LOG("EP", Sev::inform) << "Python NOT initialised";
74#endif
75
76 core_log->add_sink(
77 LogManager::createSink(LogManager::getStrmSync(), "EPSYNC"));
78 LogManager::setLog("EPSYNC");
79 MOFEM_LOG_TAG("EPSYNC", "ep");
80
81 try {
82
83 // Get mesh file
84 PetscBool flg = PETSC_FALSE;
85 char mesh_file_name[255] = "";
86 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-my_file", mesh_file_name,
87 255, &flg);
88 if (!flg) {
89 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-file_name",
90 mesh_file_name, 255, &flg);
91 }
92 if (!flg) {
93 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
94 "Missing required mesh option: use -file_name or -my_file");
95 }
96 char restart_file[255] = "no_restart_file";
97 PetscBool restart_flg = PETSC_TRUE;
98 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-restart", restart_file, 255,
99 &restart_flg);
100 double time = 0;
101 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-time", &time, PETSC_NULLPTR);
102
103
104 MOFEM_LOG("EP", Sev::inform) << "Mesh file: " << mesh_file_name;
105 MOFEM_LOG("EP", Sev::inform) << "Restart file: " << restart_file;
106 MOFEM_LOG("EP", Sev::inform) << "Time: " << time;
107
108 // Register DM Manager
109 DMType dm_name = "DMMOFEM";
110 CHKERR DMRegister_MoFEM(dm_name);
111 DMType dm_name_mg = "DMMOFEM_MG";
112 CHKERR DMRegister_MGViaApproxOrders(dm_name_mg);
113
114 // Create MoAB database
115 moab::Core moab_core;
116 moab::Interface &moab = moab_core;
117
118 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
119 if (pcomm == NULL)
120 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
121 // Read mesh to MOAB
122 PetscBool fully_distributed = PETSC_FALSE;
123 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-fully_distributed",
124 &fully_distributed, PETSC_NULLPTR);
125 if (fully_distributed) {
126 const char *option;
127 if (pcomm->proc_config().proc_size() == 1)
128 option = "";
129 else
130 option = "PARALLEL=READ_PART;"
131 "PARALLEL_RESOLVE_SHARED_ENTS;"
132 "PARTITION=PARALLEL_PARTITION";
133 CHKERR moab.load_file(mesh_file_name, 0, option);
134 } else {
135 CHKERR CommInterface::loadFileRootProcAllRestDistributed(
136 moab, mesh_file_name, SPACE_DIM);
137 }
138
139 // Create MoFEM database and link it to MoAB
140 MOFEM_LOG("EP", Sev::inform) << "Initialise MoFEM database";
141 MoFEM::Core mofem_core(moab);
142 MoFEM::Interface &m_field = mofem_core;
143 MOFEM_LOG("EP", Sev::inform) << "Initialise MoFEM database <- done";
144
145 CHKERR m_field.getInterface<MeshsetsManager>()->setMeshsetFromFile();
146
147 // Register mofem DM
148 CHKERR DMRegister_MoFEM("DMMOFEM");
149
150 BitRefLevel bit_level0 = BitRefLevel().set(0);
151 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
152 0, 3, bit_level0);
153
154 // Data stuctures
155 EshelbianCore ep(m_field);
156
157 auto meshset_ptr = get_temp_meshset_ptr(moab);
158 CHKERR moab.add_entities(
159 *meshset_ptr, CommInterface::getPartEntities(moab, pcomm->rank()));
160 auto get_adj = [&](Range ents, int dim) {
161 Range adj;
162 CHKERR moab.get_adjacencies(ents, dim, false, adj,
163 moab::Interface::UNION);
164 return adj;
165 };
166 CHKERR moab.add_entities(
167 *meshset_ptr,
168 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
169 .subset_by_dimension(SPACE_DIM),
170 2));
171 CHKERR moab.add_entities(
172 *meshset_ptr,
173 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
174 .subset_by_dimension(SPACE_DIM),
175 1));
176 CHKERR moab.add_entities(
177 *meshset_ptr,
178 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
179 .subset_by_dimension(SPACE_DIM),
180 0));
181
182 // create growing crack surface meshset
183 CHKERR ep.createCrackSurfaceMeshset();
184
185 CHKERR ep.getSpatialDispBc();
186 CHKERR ep.getSpatialRotationBc();
187 CHKERR ep.getSpatialTractionBc();
188 CHKERR ep.getSpatialTractionFreeBc();
189 CHKERR ep.getExternalStrain();
190 CHKERR ep.setBlockTagsOnSkin();
191
192 CHKERR ep.createExchangeVectors(Sev::inform);
193
194 CHKERR ep.projectMaterialTags();
195 CHKERR ep.addFields(*meshset_ptr);
196 CHKERR ep.projectGeometry(*meshset_ptr, 0.0);
197 CHKERR ep.addVolumeFiniteElement(*meshset_ptr);
198 CHKERR ep.addBoundaryFiniteElement(*meshset_ptr);
199 CHKERR ep.addDMs();
200
201 const char *list_materials[EshelbianCore::MaterialModel::LastMaterial] = {
202 "stvenant_kirchhoff", "mooney_rivlin", "hencky", "neo_hookean", "meta"};
203 PetscInt choice_material = EshelbianCore::materialModel;
204 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-material",
205 list_materials,
206 EshelbianCore::MaterialModel::LastMaterial,
207 &choice_material, PETSC_NULLPTR);
208 EshelbianCore::materialModel =
209 static_cast<EshelbianCore::MaterialModel>(choice_material);
210
211 auto material_model_names_dictionary = std::map<std::string, std::string>{
212 {"grad", "stretchH1AtPts"}, {"P", "PAtPts"}, {"P_dF", "PAtPts_du"}};
213 MatOps::MatElastic::useDeformationGradient =
214 true; // by default use deformation gradient as common data
215
216 switch (choice_material) {
217 case EshelbianCore::MaterialModel::StVenantKirchhoff: {
218 MOFEM_LOG("EP", Sev::inform) << "St. Venant-Kirchhoff material model";
219 MOFEM_LOG("EP", Sev::inform) << "Meta Material model";
220 ep.dataAtPts = boost::make_shared<DataAtIntegrationPts>();
221 auto physical_equations_ptr = MatOps::createMatOpsPhysicalEquationsPtr<
222 MatOps::ELASTICITY::STVENANTKIRCHHOFF, MatOps::MODEL_3D>(
223 MatOps::createMatOpsDataPtr(ep.dataAtPts->matOpsData,
224 material_model_names_dictionary),
225 MatOps::MatOpsTagsRegistry::setTagName("MetaElastic", 0));
226 CHKERR physical_equations_ptr->getOptions(&m_field);
227 CHKERR physical_equations_ptr->recordTape();
228 CHKERR ep.addMaterial_Core(physical_equations_ptr);
229 };
230 case EshelbianCore::MaterialModel::MooneyRivlin: {
231 MOFEM_LOG("EP", Sev::inform) << "Mooney-Rivlin material model";
232 ep.dataAtPts = boost::make_shared<DataAtIntegrationPts>();
233 auto physical_equations_ptr = MatOps::createMatOpsPhysicalEquationsPtr<
234 MatOps::ELASTICITY::MOONEYRIVLINWRIGGERSEQ63, MatOps::MODEL_3D>(
235 MatOps::createMatOpsDataPtr(ep.dataAtPts->matOpsData,
236 material_model_names_dictionary),
237 MatOps::MatOpsTagsRegistry::setTagName("MetaElastic", 0));
238 CHKERR physical_equations_ptr->getOptions(&m_field);
239 CHKERR physical_equations_ptr->recordTape();
240 CHKERR ep.addMaterial_Core(physical_equations_ptr);
241 }; break;
242 case EshelbianCore::MaterialModel::Hencky:
243 MOFEM_LOG("EP", Sev::inform) << "Hencky material model";
244 CHKERR ep.addMaterial_Hencky(5., 0.25);
245 break;
246 case EshelbianCore::MaterialModel::Neohookean:
247 MOFEM_LOG("EP", Sev::inform) << "Neo-Hookean material model";
248 CHKERR ep.addMaterial_HMHNeohookean(adolc_tag, 1.0, 1.0);
249 break;
250 case EshelbianCore::MaterialModel::Meta: {
251 MOFEM_LOG("EP", Sev::inform) << "Meta Material model";
252 ep.dataAtPts = boost::make_shared<DataAtIntegrationPts>();
253 auto physical_equations_ptr =
254 MatOps::createMatOpsPhysicalEquationsPtr<MatOps::ELASTICITY::META,
255 MatOps::MODEL_3D>(
256 MatOps::createMatOpsDataPtr(ep.dataAtPts->matOpsData,
257 material_model_names_dictionary),
258 MatOps::MatOpsTagsRegistry::setTagName("MetaElastic", 0));
259 CHKERR physical_equations_ptr->getOptions(&m_field);
260 CHKERR physical_equations_ptr->recordTape();
261 CHKERR ep.addMaterial_Core(physical_equations_ptr);
262 }; break;
263 default:
264 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown material");
265 break;
266 }
267
268 CHKERR ep.setElasticElementOps(adolc_tag);
269 CHKERR ep.setElasticElementToTs(ep.dmElastic);
270
271 auto x_elastic = createDMVector(ep.dmElastic);
272 auto f_elastic = vectorDuplicate(x_elastic);
273 CHKERR VecSetOption(f_elastic, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
274
275 PetscInt start_step = 0;
276 if (restart_flg) {
277 PetscViewer viewer;
278 CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, restart_file,
279 FILE_MODE_READ, &viewer);
280 CHKERR VecLoad(x_elastic, viewer);
281 CHKERR PetscViewerDestroy(&viewer);
282 CHKERR VecGhostUpdateBegin(x_elastic, INSERT_VALUES, SCATTER_FORWARD);
283 CHKERR VecGhostUpdateEnd(x_elastic, INSERT_VALUES, SCATTER_FORWARD);
284
285 std::string restart_file_str(restart_file);
286 std::regex restart_pattern(R"(restart_([1-9]\d*)\.dat)");
287 std::smatch match;
288 if (std::regex_search(restart_file_str, match, restart_pattern)) {
289 start_step = std::stoi(match[1]);
290 } else {
291 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
292 "Restart file name must be in the format restart_##.dat");
293 }
294 }
295
296 auto ts_elastic = createTS(PETSC_COMM_WORLD);
297 CHKERR TSSetType(ts_elastic, TSBEULER);
298 CHKERR TSAdaptRegister(TSADAPTMOFEM, TSAdaptCreateMoFEM);
299 TSAdapt adapt;
300 CHKERR TSGetAdapt(ts_elastic, &adapt);
301 CHKERR TSAdaptSetType(adapt, TSADAPTNONE);
302
303 switch (EshelbianCore::solverType) {
304 case EshelbianCore::SolverType::TimeSolver:
305 MOFEM_LOG("EP", Sev::inform) << "Solver type: TS";
306 CHKERR TSSetTime(ts_elastic, time);
307 CHKERR TSSetStepNumber(ts_elastic, start_step);
308 CHKERR ep.solveElastic(ts_elastic, x_elastic);
309 break;
310 case EshelbianCore::SolverType::DynamicRelaxation:
311 MOFEM_LOG("EP", Sev::inform) << "Solver type: Dynamic Relaxation";
312 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
313 SCATTER_REVERSE);
314 CHKERR ep.solveDynamicRelaxation(ts_elastic, x_elastic, start_step, time);
315 break;
316 case EshelbianCore::SolverType::Cohesive:
317 MOFEM_LOG("EP", Sev::inform) << "Solver type: Cohesive";
318 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
319 SCATTER_REVERSE);
320 CHKERR ep.solveCohesiveCrackGrowth(ts_elastic, x_elastic, start_step,
321 time);
322 break;
323 case EshelbianCore::SolverType::LoadFactor:
324 MOFEM_LOG("EP", Sev::inform) << "Solver type: Load Factor";
325 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
326 SCATTER_REVERSE);
327 CHKERR ep.solveLoadFactor(ts_elastic, x_elastic, start_step, time);
328 break;
329 case EshelbianCore::SolverType::ShapeOptimization:
330 MOFEM_LOG("EP", Sev::inform) << "Solver type: Shape Optimization";
331 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
332 SCATTER_REVERSE);
333 CHKERR ep.solveSchapeOptimisation(ts_elastic, x_elastic, start_step,
334 time);
335 break;
336 case EshelbianCore::SolverType::TestTopologicalDerivative:
337 MOFEM_LOG("EP", Sev::inform)
338 << "Solver type: Test Topological Derivative";
339 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
340 SCATTER_REVERSE);
341 CHKERR ep.solveTestTopologicalDerivative(ts_elastic, x_elastic,
342 start_step, time);
343 break;
344 default:
345 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
346 "Unknown solver type");
347 break;
348 }
349
350 }
351 CATCH_ERRORS;
352
353 // finish work cleaning memory, getting statistics, etc.
354 MoFEM::Core::Finalize();
355
356#ifdef ENABLE_PYTHON_BINDING
357 if (Py_FinalizeEx() < 0) {
358 exit(120);
359 }
360#endif
361}
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:50
adolc_tag
constexpr int adolc_tag
Definition ep.cpp:25
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
MatOps::MODEL_3D
@ MODEL_3D
Definition MatOps.hpp:171
MatOps::createMatOpsPhysicalEquationsPtr
boost::shared_ptr< PhysicalEquations > createMatOpsPhysicalEquationsPtr(boost::shared_ptr< MatOpsData > mat_ops_data_ptr, int tag)
MatOps::createMatOpsDataPtr
boost::shared_ptr< MatOpsData > createMatOpsDataPtr()
Definition MatOps.cpp:671
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: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::materialModel
static enum MaterialModel materialModel
Definition EshelbianCore.hpp:46
MatOps::MatElastic::useDeformationGradient
static bool useDeformationGradient
Definition MatElastic.hpp:18
MatOps::MatOpsTagsRegistry::setTagName
static int setTagName(std::string name, int tag=-1)
Definition MatOps.cpp:30
MoFEM::BitRefManager
Managing BitRefLevels.
Definition BitRefManager.hpp:21
MoFEM::CommInterface::getPartEntities
static Range getPartEntities(moab::Interface &moab, int part)
Definition CommInterface.cpp:1418
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:1205
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::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 25 of file ep.cpp.

◆ help

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

Definition at line 50 of file ep.cpp.

Generated by Doxygen 1.9.8 and hosted at University of Glasgow