v0.15.0
Loading...
Searching...
No Matches
ep.cpp
Go to the documentation of this file.
1/**
2 * \file ep.cpp
3 * \example ep.cpp
4 *
5 * \brief Implementation of mix-element for Large strains
6 *
7 * \todo Implementation of plasticity
8 *
9 */
10
11/* This file is part of MoFEM.
12 * MoFEM is free software: you can redistribute it and/or modify it under
13 * the terms of the GNU Lesser General Public License as published by the
14 * Free Software Foundation, either version 3 of the License, or (at your
15 * option) any later version.
16 *
17 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
18 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
19 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
20 * License for more details.
21 *
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
24
25constexpr int adolc_tag = 1;
26
27#define SINGULARITY
28
29#include <MoFEM.hpp>
30using namespace MoFEM;
31using namespace boost::multi_index;
32using namespace boost::multiprecision;
33using namespace boost::numeric;
34#include <cmath>
35
36#include <cholesky.hpp>
37#ifdef ENABLE_PYTHON_BINDING
38 #include <boost/python.hpp>
39 #include <boost/python/def.hpp>
40 #include <boost/python/numpy.hpp>
41namespace bp = boost::python;
42namespace np = boost::python::numpy;
43#endif
44
46using namespace EshelbianPlasticity;
47
48static char help[] = "...\n\n";
49
50int main(int argc, char *argv[]) {
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"));
60 MOFEM_LOG_TAG("EP", "ep");
61 core_log->add_sink(
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(
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];
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 // Register DM Manager
96 DMType dm_name = "DMMOFEM";
97 CHKERR DMRegister_MoFEM(dm_name);
98 DMType dm_name_mg = "DMMOFEM_MG";
100
101 // Create MoAB database
102 moab::Core moab_core;
103 moab::Interface &moab = moab_core;
104
105 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
106 if (pcomm == NULL)
107 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
108 // Read mesh to MOAB
109 PetscBool fully_distributed = PETSC_FALSE;
110 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-fully_distributed",
111 &fully_distributed, PETSC_NULLPTR);
112 if (fully_distributed) {
113 const char *option;
114 if (pcomm->proc_config().proc_size() == 1)
115 option = "";
116 else
117 option = "PARALLEL=READ_PART;"
118 "PARALLEL_RESOLVE_SHARED_ENTS;"
119 "PARTITION=PARALLEL_PARTITION";
120 CHKERR moab.load_file(mesh_file_name, 0, option);
121 } else {
124 }
125
126 // Create MoFEM database and link it to MoAB
127 MOFEM_LOG("EP", Sev::inform) << "Initialise MoFEM database";
128 MoFEM::Core mofem_core(moab);
129 MoFEM::Interface &m_field = mofem_core;
130 MOFEM_LOG("EP", Sev::inform) << "Initialise MoFEM database <- done";
131
132 CHKERR m_field.getInterface<MeshsetsManager>()->setMeshsetFromFile();
133
134 // Register mofem DM
135 CHKERR DMRegister_MoFEM("DMMOFEM");
136
137 BitRefLevel bit_level0 = BitRefLevel().set(0);
138 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
139 0, 3, bit_level0);
140
141 // Data stuctures
142 EshelbianCore ep(m_field);
143
144 auto meshset_ptr = get_temp_meshset_ptr(moab);
145 CHKERR moab.add_entities(
146 *meshset_ptr, CommInterface::getPartEntities(moab, pcomm->rank()));
147 auto get_adj = [&](Range ents, int dim) {
148 Range adj;
149 CHKERR moab.get_adjacencies(ents, dim, false, adj,
150 moab::Interface::UNION);
151 return adj;
152 };
153 CHKERR moab.add_entities(
154 *meshset_ptr,
155 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
156 .subset_by_dimension(SPACE_DIM),
157 2));
158 CHKERR moab.add_entities(
159 *meshset_ptr,
160 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
161 .subset_by_dimension(SPACE_DIM),
162 1));
163 CHKERR moab.add_entities(
164 *meshset_ptr,
165 get_adj(CommInterface::getPartEntities(moab, pcomm->rank())
166 .subset_by_dimension(SPACE_DIM),
167 0));
168
169 // create growing crack surface meshset
171
178
179 CHKERR ep.createExchangeVectors(Sev::inform);
180
181 CHKERR ep.addFields(*meshset_ptr);
182 CHKERR ep.projectGeometry(*meshset_ptr);
183 CHKERR ep.addVolumeFiniteElement(*meshset_ptr);
184 CHKERR ep.addBoundaryFiniteElement(*meshset_ptr);
185 CHKERR ep.addDMs();
186
187 enum MaterialModel {
188 StVenantKirchhoff,
189 MooneyRivlin,
190 Hencky,
191 Neohookean,
192 LastMaterial
193 };
194 const char *list_materials[LastMaterial] = {
195 "stvenant_kirchhoff", "mooney_rivlin", "hencky", "neo_hookean"};
196 PetscInt choice_material = MooneyRivlin;
197 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-material", list_materials,
198 LastMaterial, &choice_material, PETSC_NULLPTR);
199
200 switch (choice_material) {
201 case StVenantKirchhoff:
202 MOFEM_LOG("EP", Sev::inform) << "StVenantKirchhoff material model";
204 MU(1, 0.25), 0);
205 break;
206 case MooneyRivlin:
207 MOFEM_LOG("EP", Sev::inform) << "MooneyRivlin material model";
208 MOFEM_LOG("EP", Sev::inform) << "MU(1, 0)/2 = " << MU(1, 0) / 2;
209 MOFEM_LOG("EP", Sev::inform) << "LAMBDA(1, 0) = " << LAMBDA(1, 0);
211 LAMBDA(1, 0), 0);
212 break;
213 case Hencky:
214 MOFEM_LOG("EP", Sev::inform) << "Hencky material model";
215 CHKERR ep.addMaterial_Hencky(5., 0.25);
216 break;
217 case Neohookean:
218 MOFEM_LOG("EP", Sev::inform) << "Neo-Hookean material model";
220 break;
221 default:
222 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown material");
223 break;
224 }
225
228
229 SmartPetscObj<Vec> x_elastic, f_elastic;
231 f_elastic = vectorDuplicate(x_elastic);
232 CHKERR VecSetOption(f_elastic, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
233
234 PetscInt start_step = 0;
235 if (restart_flg) {
236 PetscViewer viewer;
237 CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, restart_file,
238 FILE_MODE_READ, &viewer);
239 CHKERR VecLoad(x_elastic, viewer);
240 CHKERR PetscViewerDestroy(&viewer);
241 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
242 SCATTER_REVERSE);
243
244 std::string restart_file_str(restart_file);
245 std::regex restart_pattern(R"(restart_([1-9]\d*)\.dat)");
246 std::smatch match;
247 if (std::regex_match(restart_file_str, match, restart_pattern)) {
248 start_step = std::stoi(match[1]);
249 } else {
250 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
251 "Restart file name must be in the format restart_##.dat");
252 }
253 }
254
255 auto ts_elastic = createTS(PETSC_COMM_WORLD);
256 CHKERR TSSetType(ts_elastic, TSBEULER);
257 CHKERR TSAdaptRegister(TSADAPTMOFEM, TSAdaptCreateMoFEM);
258 TSAdapt adapt;
259 CHKERR TSGetAdapt(ts_elastic, &adapt);
260 CHKERR TSAdaptSetType(adapt, TSADAPTNONE);
261 CHKERR TSSetTime(ts_elastic, time);
262 CHKERR TSSetStepNumber(ts_elastic, start_step);
263
264 if(ep.dynamicRelaxation) {
265 CHKERR ep.solveDynamicRelaxation(ts_elastic, x_elastic);
266 } else {
267 CHKERR ep.solveElastic(ts_elastic, x_elastic);
268 }
269
270 }
272
273 // finish work cleaning memory, getting statistics, etc.
275
276#ifdef ENABLE_PYTHON_BINDING
277 if (Py_FinalizeEx() < 0) {
278 exit(120);
279 }
280#endif
281}
Eshelbian plasticity interface.
#define TSADAPTMOFEM
Definition TsCtx.hpp:10
int main()
constexpr int SPACE_DIM
cholesky decomposition
#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:48
constexpr int adolc_tag
Definition ep.cpp:25
#define MU(E, NU)
Definition fem_tools.h:23
#define LAMBDA(E, NU)
Definition fem_tools.h:22
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
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.
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition DMMoFEM.cpp:1157
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.
char mesh_file_name[255]
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
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
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.
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:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
Deprecated interface functions.
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.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.