v0.14.0
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 
25 constexpr int adolc_tag = 1;
26 
27 #define SINGULARITY
28 
29 #include <MoFEM.hpp>
30 using namespace MoFEM;
31 using namespace boost::multi_index;
32 using namespace boost::multiprecision;
33 using namespace boost::numeric;
34 
35 #include <cholesky.hpp>
36 #ifdef PYTHON_SDF
37  #include <boost/python.hpp>
38  #include <boost/python/def.hpp>
39  #include <boost/python/numpy.hpp>
40 namespace bp = boost::python;
41 namespace np = boost::python::numpy;
42 #endif
43 
44 #include <EshelbianPlasticity.hpp>
45 #include <TSEPAdapt.hpp>
46 using namespace EshelbianPlasticity;
47 
48 static char help[] = "...\n\n";
49 
50 int 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"));
59  LogManager::setLog("EP");
60  MOFEM_LOG_TAG("EP", "ep");
61 
62 #ifdef PYTHON_SDF
63  Py_Initialize();
64  np::initialize();
65  MOFEM_LOG("EP", Sev::inform) << "Python initialised";
66 #else
67  MOFEM_LOG("EP", Sev::inform) << "Python NOT initialised";
68 #endif
69 
70  core_log->add_sink(
72  LogManager::setLog("EPSYNC");
73  MOFEM_LOG_TAG("EPSYNC", "ep");
74 
75  try {
76 
77  // Get mesh file
78  PetscBool flg = PETSC_TRUE;
79  char mesh_file_name[255];
80  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
81  255, &flg);
82  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-file_name", mesh_file_name,
83  255, &flg);
84  char restart_file[255];
85  PetscBool restart_flg = PETSC_TRUE;
86  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-restart", restart_file, 255,
87  &restart_flg);
88  double time = 0;
89  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-time", &time, PETSC_NULL);
90 
91  // Register DM Manager
92  DMType dm_name = "DMMOFEM";
93  CHKERR DMRegister_MoFEM(dm_name);
94  DMType dm_name_mg = "DMMOFEM_MG";
96 
97  // Create MoAB database
98  moab::Core moab_core;
99  moab::Interface &moab = moab_core;
100 
101  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
102  if (pcomm == NULL)
103  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
104  // Read mesh to MOAB
105  const char *option;
106  if (pcomm->proc_config().proc_size() == 1)
107  option = "";
108  else
109  option = "PARALLEL=READ_PART;"
110  "PARALLEL_RESOLVE_SHARED_ENTS;"
111  "PARTITION=PARALLEL_PARTITION";
112  CHKERR moab.load_file(mesh_file_name, 0, option);
113 
114  // Create MoFEM database and link it to MoAB
115  MoFEM::Core mofem_core(moab);
116  MoFEM::Interface &m_field = mofem_core;
117 
118  // Register mofem DM
119  CHKERR DMRegister_MoFEM("DMMOFEM");
120 
121  BitRefLevel bit_level0 = BitRefLevel().set(0);
122  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
123  0, 3, bit_level0);
124 
125  // Data stuctures
126  EshelbianCore ep(m_field);
127 
132 
133  CHKERR ep.addFields();
134  CHKERR ep.projectGeometry();
137  CHKERR ep.addDMs();
138 
139  enum MaterialModel {
140  StVenantKirchhoff,
141  MooneyRivlin,
142  Hencky,
143  Neohookean,
144  LastMaterial
145  };
146  const char *list_materials[LastMaterial] = {
147  "stvenant_kirchhoff", "mooney_rivlin", "hencky", "neo_hookean"};
148  PetscInt choice_material = MooneyRivlin;
149  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-material", list_materials,
150  LastMaterial, &choice_material, PETSC_NULL);
151 
152  switch (choice_material) {
153  case StVenantKirchhoff:
154  MOFEM_LOG("EP", Sev::inform) << "StVenantKirchhoff material model";
156  MU(1, 0.25), 0);
157  break;
158  case MooneyRivlin:
159  MOFEM_LOG("EP", Sev::inform) << "MooneyRivlin material model";
160  MOFEM_LOG("EP", Sev::inform) << "MU(1, 0)/2 = " << MU(1, 0) / 2;
161  MOFEM_LOG("EP", Sev::inform) << "LAMBDA(1, 0) = " << LAMBDA(1, 0);
163  LAMBDA(1, 0), 0);
164  break;
165  case Hencky:
166  MOFEM_LOG("EP", Sev::inform) << "Hencky material model";
167  CHKERR ep.addMaterial_Hencky(5., 0.25);
168  break;
169  case Neohookean:
170  MOFEM_LOG("EP", Sev::inform) << "Neohookean material model";
172  break;
173  default:
174  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown material");
175  break;
176  }
177 
180 
181  SmartPetscObj<Vec> x_elastic, f_elastic;
183  f_elastic = vectorDuplicate(x_elastic);
184  CHKERR VecSetOption(f_elastic, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
185 
186  if (restart_flg) {
187  PetscViewer viewer;
188  CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, restart_file,
189  FILE_MODE_READ, &viewer);
190  CHKERR VecLoad(x_elastic, viewer);
191  CHKERR PetscViewerDestroy(&viewer);
192  CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x_elastic, INSERT_VALUES,
193  SCATTER_REVERSE);
194  }
195 
196  auto ts_elastic = createTS(PETSC_COMM_WORLD);
197  CHKERR TSSetType(ts_elastic, TSBEULER);
198  CHKERR TSAdaptRegister(TSADAPTEP, TSAdaptCreate_EP);
199  TSAdapt adapt;
200  CHKERR TSGetAdapt(ts_elastic, &adapt);
201  CHKERR TSAdaptSetType(adapt, TSADAPTEP);
202  // CHKERR TSSetFromOptions(ts_elastic);
203 
204  CHKERR TSSetTime(ts_elastic, time);
205  CHKERR ep.solveElastic(ts_elastic, x_elastic);
206  }
207  CATCH_ERRORS;
208 
209  // finish work cleaning memory, getting statistics, etc.
211 
212 #ifdef PYTHON_SDF
213  if (Py_FinalizeEx() < 0) {
214  exit(120);
215  }
216 #endif
217 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
EshelbianPlasticity::EshelbianCore::addDMs
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0), const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:1101
EshelbianPlasticity::EshelbianCore::projectGeometry
MoFEMErrorCode projectGeometry(const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:707
MoFEM::createTS
auto createTS(MPI_Comm comm)
Definition: PetscSmartObj.hpp:249
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
EshelbianPlasticity::EshelbianCore::addMaterial_Hencky
MoFEMErrorCode addMaterial_Hencky(double E, double nu)
Definition: EshelbianADOL-C.cpp:491
EshelbianPlasticity::EshelbianCore::addBoundaryFiniteElement
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:912
EshelbianPlasticity::EshelbianCore::setElasticElementToTs
MoFEMErrorCode setElasticElementToTs(DM dm)
Definition: EshelbianPlasticity.cpp:2192
EshelbianPlasticity::EshelbianCore::setElasticElementOps
MoFEMErrorCode setElasticElementOps(const int tag)
Definition: EshelbianPlasticity.cpp:2169
LAMBDA
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
EshelbianPlasticity
Definition: CGGTonsorialBubbleBase.hpp:11
MoFEM.hpp
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:523
EshelbianPlasticity::EshelbianCore::getSpatialTractionFreeBc
MoFEMErrorCode getSpatialTractionFreeBc(const EntityHandle meshset=0)
Definition: EshelbianPlasticity.hpp:1061
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
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:298
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
adolc_tag
constexpr int adolc_tag
Definition: ep.cpp:25
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1167
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
TSEPAdapt.hpp
Step adaptation.
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
EshelbianPlasticity::EshelbianCore::addVolumeFiniteElement
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:850
EshelbianPlasticity::EshelbianCore::dmElastic
SmartPetscObj< DM > dmElastic
Elastic problem.
Definition: EshelbianPlasticity.hpp:964
MoFEM::DMRegister_MGViaApproxOrders
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:302
TSAdaptCreate_EP
PetscErrorCode TSAdaptCreate_EP(TSAdapt adapt)
Definition: TSEPAdapt.hpp:87
EshelbianPlasticity::EshelbianCore::addMaterial_HMHHStVenantKirchhoff
MoFEMErrorCode addMaterial_HMHHStVenantKirchhoff(const int tape, const double lambda, const double mu, const double sigma_y)
Definition: EshelbianADOL-C.cpp:464
MoFEM::LogManager::getStrmSync
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition: LogManager.cpp:348
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
EshelbianPlasticity::EshelbianCore::addMaterial_HMHNeohookean
MoFEMErrorCode addMaterial_HMHNeohookean(const int tape, const double c10, const double K)
Definition: EshelbianADOL-C.cpp:482
EshelbianPlasticity::EshelbianCore::solveElastic
MoFEMErrorCode solveElastic(TS ts, Vec x)
Definition: EshelbianPlasticity.cpp:2225
EshelbianPlasticity.hpp
Eshelbian plasticity interface.
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
cholesky.hpp
cholesky decomposition
EshelbianPlasticity::EshelbianCore
Definition: EshelbianPlasticity.hpp:894
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
help
static char help[]
Definition: ep.cpp:48
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
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
EshelbianPlasticity::EshelbianCore::addFields
MoFEMErrorCode addFields(const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:483
TSADAPTEP
#define TSADAPTEP
Definition: TSEPAdapt.hpp:15
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
main
int main(int argc, char *argv[])
Definition: ep.cpp:50
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::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< Vec >
EshelbianPlasticity::EshelbianCore::getSpatialTractionBc
MoFEMErrorCode getSpatialTractionBc()
Definition: EshelbianPlasticity.cpp:2599
EshelbianPlasticity::EshelbianCore::addMaterial_HMHMooneyRivlin
MoFEMErrorCode addMaterial_HMHMooneyRivlin(const int tape, const double alpha, const double beta, const double lambda, const double sigma_y)
Definition: EshelbianADOL-C.cpp:473
EshelbianPlasticity::EshelbianCore::getSpatialRotationBc
MoFEMErrorCode getSpatialRotationBc()
Definition: EshelbianPlasticity.hpp:1039
EshelbianPlasticity::EshelbianCore::getSpatialDispBc
MoFEMErrorCode getSpatialDispBc()
[Getting norms]
Definition: EshelbianPlasticity.cpp:2561
MU
#define MU(E, NU)
Definition: fem_tools.h:23