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"));
62 MOFEM_LOG_TAG("EP", "ep");
63 core_log->add_sink(
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(
78 LogManager::setLog("EPSYNC");
79 MOFEM_LOG_TAG("EPSYNC", "ep");
80
81 try {
82
84
85 // Get mesh file
86 PetscBool flg = PETSC_FALSE;
87 char mesh_file_name[255] = "";
88 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-my_file", mesh_file_name,
89 255, &flg);
90 if (!flg) {
91 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-file_name",
92 mesh_file_name, 255, &flg);
93 }
94 if (!flg) {
95 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
96 "Missing required mesh option: use -file_name or -my_file");
97 }
98 char restart_file[255] = "no_restart_file";
99 PetscBool restart_flg = PETSC_TRUE;
100 CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-restart", restart_file, 255,
101 &restart_flg);
102 double time = 0;
103 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-time", &time, PETSC_NULLPTR);
104
105
106 MOFEM_LOG("EP", Sev::inform) << "Mesh file: " << mesh_file_name;
107 MOFEM_LOG("EP", Sev::inform) << "Restart file: " << restart_file;
108 MOFEM_LOG("EP", Sev::inform) << "Time: " << time;
109
110 // Register DM Manager
111 DMType dm_name = "DMMOFEM";
112 CHKERR DMRegister_MoFEM(dm_name);
113 DMType dm_name_mg = "DMMOFEM_MG";
115
116 // Create MoAB database
117 moab::Core moab_core;
118 moab::Interface &moab = moab_core;
119
120 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
121 if (pcomm == NULL)
122 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
123 // Read mesh to MOAB
124 PetscBool fully_distributed = PETSC_FALSE;
125 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-fully_distributed",
126 &fully_distributed, PETSC_NULLPTR);
127 if (fully_distributed) {
128 const char *option;
129 if (pcomm->proc_config().proc_size() == 1)
130 option = "";
131 else
132 option = "PARALLEL=READ_PART;"
133 "PARALLEL_RESOLVE_SHARED_ENTS;"
134 "PARTITION=PARALLEL_PARTITION";
135 CHKERR moab.load_file(mesh_file_name, 0, option);
136 } else {
138 moab, mesh_file_name, SPACE_DIM);
139 }
140
141 // Create MoFEM database and link it to MoAB
142 MOFEM_LOG("EP", Sev::inform) << "Initialise MoFEM database";
143 MoFEM::Core mofem_core(moab);
144 MoFEM::Interface &m_field = mofem_core;
145 MOFEM_LOG("EP", Sev::inform) << "Initialise MoFEM database <- done";
146 CHKERR JsonConfigManager::addCanonicalAttributeNames({"kn", "kt"}); //TODO: remove this line and add kn kt to the core JsonConfigManagerCanonicalAttributes.hpp
147 CHKERR m_field.getInterface<MeshsetsManager>()->setMeshsetFromFile();
148
149 // Register mofem DM
150 CHKERR DMRegister_MoFEM("DMMOFEM");
151
152 BitRefLevel bit_level0 = BitRefLevel().set(0);
153 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
154 0, 3, bit_level0);
155
156 // Data stuctures
157 EshelbianCore ep(m_field);
158
159 auto meshset_ptr = get_temp_meshset_ptr(moab);
160 CHKERR moab.add_entities(
161 *meshset_ptr, CommInterface::getPartEntities(moab, pcomm->rank()));
162 auto get_adj = [&](Range ents, int dim) {
163 Range adj;
164 CHKERR moab.get_adjacencies(ents, dim, false, adj,
165 moab::Interface::UNION);
166 return adj;
167 };
168 CHKERR moab.add_entities(
169 *meshset_ptr,
170 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
171 .subset_by_dimension(SPACE_DIM),
172 2));
173 CHKERR moab.add_entities(
174 *meshset_ptr,
175 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
176 .subset_by_dimension(SPACE_DIM),
177 1));
178 CHKERR moab.add_entities(
179 *meshset_ptr,
180 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
181 .subset_by_dimension(SPACE_DIM),
182 0));
183
184 // create growing crack surface meshset
185 CHKERR ep.createCrackSurfaceMeshset();
186
187 CHKERR ep.getSpatialDispBc();
188 CHKERR ep.getSpatialRotationBc();
189 CHKERR ep.getSpatialTractionBc();
190 CHKERR ep.getSpatialTractionFreeBc();
191 CHKERR ep.getExternalStrain();
192 CHKERR ep.setBlockTagsOnSkin();
193
194 CHKERR ep.createExchangeVectors(Sev::inform);
195
196 CHKERR ep.projectMaterialTags();
197 CHKERR ep.addFields(*meshset_ptr);
198 CHKERR ep.projectGeometry(*meshset_ptr, 0.0);
199 CHKERR ep.addVolumeFiniteElement(*meshset_ptr);
200 CHKERR ep.addBoundaryFiniteElement(*meshset_ptr);
201 CHKERR ep.addDMs();
202
203 const char *list_materials[EshelbianCore::MaterialModel::LastMaterial] = {
204 "stvenant_kirchhoff", "mooney_rivlin", "hencky",
205 "neo_hookean", "storakers", "meta"};
206 PetscInt choice_material = EshelbianCore::materialModel;
207 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-material",
208 list_materials,
210 &choice_material, PETSC_NULLPTR);
212 static_cast<EshelbianCore::MaterialModel>(choice_material);
213
214 auto material_model_names_dictionary = std::map<std::string, std::string>{
215 {"grad", "stretchH1AtPts"}, {"P", "PAtPts"}, {"P_dF", "PAtPts_du"}};
217 true; // by default use deformation gradient as common data
218
219 switch (choice_material) {
221 MOFEM_LOG("EP", Sev::inform) << "St. Venant-Kirchhoff material model";
222 MOFEM_LOG("EP", Sev::inform) << "Meta Material model";
223 ep.dataAtPts = boost::make_shared<DataAtIntegrationPts>();
224 auto physical_equations_ptr = MatOps::createMatOpsPhysicalEquationsPtr<
226 MatOps::createMatOpsDataPtr(ep.dataAtPts->matOpsData,
227 material_model_names_dictionary),
229 CHKERR physical_equations_ptr->getOptions(&m_field);
230 CHKERR physical_equations_ptr->recordTape();
231 CHKERR ep.addMaterial_Core(physical_equations_ptr);
232 };
234 MOFEM_LOG("EP", Sev::inform) << "Mooney-Rivlin material model";
235 ep.dataAtPts = boost::make_shared<DataAtIntegrationPts>();
236 auto physical_equations_ptr = MatOps::createMatOpsPhysicalEquationsPtr<
238 MatOps::createMatOpsDataPtr(ep.dataAtPts->matOpsData,
239 material_model_names_dictionary),
241 CHKERR physical_equations_ptr->getOptions(&m_field);
242 CHKERR physical_equations_ptr->recordTape();
243 CHKERR ep.addMaterial_Core(physical_equations_ptr);
244 }; break;
246 MOFEM_LOG("EP", Sev::inform) << "Hencky material model";
247 CHKERR ep.addMaterial_Hencky(5., 0.25);
248 break;
250 MOFEM_LOG("EP", Sev::inform) << "Neo-Hookean material model";
251 // The power-law Neo-Hookean model requires K > 4*c10/3. These values
252 // are only defaults; -neo_hookean_c10 and -neo_hookean_K still override
253 // them.
254 CHKERR ep.addMaterial_HMHNeohookean(1.0, 2.0);
255 break;
257 MOFEM_LOG("EP", Sev::inform) << "Storakers material model";
258 CHKERR ep.addMaterial_HMHStorakers(0.5, 1.0, 1.0);
259 break;
261 MOFEM_LOG("EP", Sev::inform) << "Meta Material model";
262 ep.dataAtPts = boost::make_shared<DataAtIntegrationPts>();
263 auto physical_equations_ptr =
266 MatOps::createMatOpsDataPtr(ep.dataAtPts->matOpsData,
267 material_model_names_dictionary),
269 CHKERR physical_equations_ptr->getOptions(&m_field);
270 CHKERR physical_equations_ptr->recordTape();
271 CHKERR ep.addMaterial_Core(physical_equations_ptr);
272 }; break;
273 default:
274 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown material");
275 break;
276 }
277
278 CHKERR ep.setElasticElementOps(adolc_tag);
279 CHKERR ep.setElasticElementToTs(ep.dmElastic);
280
281 auto x_elastic = createDMVector(ep.dmElastic);
282 auto f_elastic = vectorDuplicate(x_elastic);
283 CHKERR VecSetOption(f_elastic, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
284
285 PetscInt start_step = 0;
286 if (restart_flg) {
287 PetscViewer viewer;
288 CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, restart_file,
289 FILE_MODE_READ, &viewer);
290 CHKERR VecLoad(x_elastic, viewer);
291 CHKERR PetscViewerDestroy(&viewer);
292 CHKERR VecGhostUpdateBegin(x_elastic, INSERT_VALUES, SCATTER_FORWARD);
293 CHKERR VecGhostUpdateEnd(x_elastic, INSERT_VALUES, SCATTER_FORWARD);
294
295 std::string restart_file_str(restart_file);
296 std::regex restart_pattern(R"(restart_([1-9]\d*)\.dat)");
297 std::smatch match;
298 if (std::regex_search(restart_file_str, match, restart_pattern)) {
299 start_step = std::stoi(match[1]);
300 } else {
301 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
302 "Restart file name must be in the format restart_##.dat");
303 }
304 }
305
306 auto ts_elastic = createTS(PETSC_COMM_WORLD);
307 CHKERR TSSetType(ts_elastic, TSBEULER);
308 CHKERR TSAdaptRegister(TSADAPTMOFEM, TSAdaptCreateMoFEM);
309 TSAdapt adapt;
310 CHKERR TSGetAdapt(ts_elastic, &adapt);
311 CHKERR TSAdaptSetType(adapt, TSADAPTNONE);
312
315 MOFEM_LOG("EP", Sev::inform) << "Solver type: TS";
316 CHKERR TSSetTime(ts_elastic, time);
317 CHKERR TSSetStepNumber(ts_elastic, start_step);
318 CHKERR ep.solveElastic(ts_elastic, x_elastic);
319 break;
321 MOFEM_LOG("EP", Sev::inform) << "Solver type: Dynamic Relaxation";
322 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
323 SCATTER_REVERSE);
324 CHKERR ep.solveDynamicRelaxation(ts_elastic, x_elastic, start_step, time);
325 break;
327 MOFEM_LOG("EP", Sev::inform) << "Solver type: Cohesive";
328 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
329 SCATTER_REVERSE);
330 CHKERR ep.solveCohesiveCrackGrowth(ts_elastic, x_elastic, start_step,
331 time);
332 break;
334 MOFEM_LOG("EP", Sev::inform) << "Solver type: Load Factor";
335 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
336 SCATTER_REVERSE);
337 CHKERR ep.solveLoadFactor(ts_elastic, x_elastic, start_step, time);
338 break;
340 MOFEM_LOG("EP", Sev::inform) << "Solver type: Shape Optimization";
341 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
342 SCATTER_REVERSE);
343 CHKERR ep.solveSchapeOptimisation(ts_elastic, x_elastic, start_step,
344 time);
345 break;
347 MOFEM_LOG("EP", Sev::inform)
348 << "Solver type: Test Topological Derivative";
349 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
350 SCATTER_REVERSE);
351 CHKERR ep.solveTestTopologicalDerivative(ts_elastic, x_elastic,
352 start_step, time);
353 break;
354 default:
355 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
356 "Unknown solver type");
357 break;
358 }
359
360 }
362
363 // finish work cleaning memory, getting statistics, etc.
365
366#ifdef ENABLE_PYTHON_BINDING
367 if (Py_FinalizeEx() < 0) {
368 exit(120);
369 }
370#endif
371}
#define TSADAPTMOFEM
Definition TsCtx.hpp:10
constexpr int SPACE_DIM
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define CHKERR
Inline error check.
static char help[]
Definition ep.cpp:50
constexpr int adolc_tag
Definition ep.cpp:25
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
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
auto createDMVector(DM dm, RowColData rc=RowColData::COL)
Get smart vector from DM.
Definition DMMoFEM.hpp:1237
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.
@ MODEL_3D
Definition MatOps.hpp:181
boost::shared_ptr< PhysicalEquations > createMatOpsPhysicalEquationsPtr(boost::shared_ptr< MatOpsData > mat_ops_data_ptr, int tag)
boost::shared_ptr< MatOpsData > createMatOpsDataPtr()
Definition MatOps.cpp:671
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
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.
Definition TsCtx.cpp:829
static enum SolverType solverType
static enum MaterialModel materialModel
static bool useDeformationGradient
static int setTagName(std::string name, int tag=-1)
Definition MatOps.cpp:30
Managing BitRefLevels.
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
Core (interface) class.
Definition Core.hpp:83
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:123
Deprecated interface functions.
static MoFEMErrorCode setMeshFileFromJson()
Set -file_name from JSON before Core is available.
static MoFEMErrorCode addCanonicalAttributeNames(const std::vector< std::string > &names)
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.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
char restart_file[255]
PetscBool restart_flg

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.