v0.9.2
lesson7_plastic.cpp
Go to the documentation of this file.
1 /**
2  * \file lesson7_plastic.cpp
3  * \example lesson7_plastic.cpp
4  *
5  * Plane stress elastic problem
6  *
7  */
8 
9 /* This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #include <MoFEM.hpp>
24 
25 using namespace MoFEM;
26 
32 
33 constexpr double young_modulus = 1e1;
34 constexpr double poisson_ratio = 0.25;
35 constexpr double sigmaY = 1;
36 constexpr double H = 1e-2;
37 constexpr double cn = H;
38 constexpr int order = 2;
39 
40 #include <ElasticOps.hpp>
41 #include <PlasticOps.hpp>
42 using namespace OpElasticTools;
43 using namespace OpPlasticTools;
44 
45 struct Example {
46 
47  Example(MoFEM::Interface &m_field) : mField(m_field) {}
48 
49  MoFEMErrorCode runProblem();
50 
51 private:
52  MoFEM::Interface &mField;
53 
54  MoFEMErrorCode setUP();
55  MoFEMErrorCode createCommonData();
56  MoFEMErrorCode bC();
57  MoFEMErrorCode OPs();
58  MoFEMErrorCode tsSolve();
59  MoFEMErrorCode postProcess();
60  MoFEMErrorCode checkResults();
61 
63  boost::shared_ptr<OpPlasticTools::CommonData> commonDataPtr;
64  boost::shared_ptr<PostProcFaceOnRefinedMesh> postProcFe;
65  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
66  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
67 };
68 
69 //! [Run problem]
72  CHKERR setUP();
73  CHKERR createCommonData();
74  CHKERR bC();
75  CHKERR OPs();
76  CHKERR tsSolve();
77  CHKERR postProcess();
78  CHKERR checkResults();
80 }
81 //! [Run problem]
82 
83 //! [Set up problem]
86  Simple *simple = mField.getInterface<Simple>();
87  // Add field
88  CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 2);
89  CHKERR simple->addDomainField("TAU", L2, AINSWORTH_LEGENDRE_BASE, 1);
90  CHKERR simple->addDomainField("EP", L2, AINSWORTH_LEGENDRE_BASE, 3);
91  CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, 2);
92  CHKERR simple->setFieldOrder("U", order);
93  CHKERR simple->setFieldOrder("TAU", order-1);
94  CHKERR simple->setFieldOrder("EP", order-1);
95  CHKERR simple->setUp();
97 }
98 //! [Set up problem]
99 
100 //! [Create common data]
103 
104  auto get_matrial_stiffens = [&](FTensor::Ddg<double, 2, 2> &t_D) {
106 
111  t_D(i, j, k, l) = 0;
112 
113  constexpr double c = young_modulus / (1 - poisson_ratio * poisson_ratio);
114  constexpr double o = poisson_ratio * c;
115 
116  t_D(0, 0, 0, 0) = c;
117  t_D(0, 0, 1, 1) = o;
118 
119  t_D(1, 1, 0, 0) = o;
120  t_D(1, 1, 1, 1) = c;
121 
122  t_D(0, 1, 0, 1) = (1 - poisson_ratio) * c;
123 
125  };
126 
127  commonDataPtr = boost::make_shared<OpPlasticTools::CommonData>();
128  CHKERR get_matrial_stiffens(commonDataPtr->tD);
129  commonDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
130  commonDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
131  commonDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
132  commonDataPtr->plasticSurfacePtr = boost::make_shared<VectorDouble>();
133  commonDataPtr->plasticFlowPtr = boost::make_shared<MatrixDouble>();
134  commonDataPtr->plasticTauPtr = boost::make_shared<VectorDouble>();
135  commonDataPtr->plasticTauDotPtr = boost::make_shared<VectorDouble>();
136  commonDataPtr->plasticStrainPtr = boost::make_shared<MatrixDouble>();
137  commonDataPtr->plasticStrainDotPtr = boost::make_shared<MatrixDouble>();
139 }
140 //! [Create common data]
141 
142 //! [Boundary condition]
145 
146  auto fix_disp = [&](const std::string blockset_name) {
147  Range fix_ents;
149  if (it->getName().compare(0, blockset_name.length(), blockset_name) ==
150  0) {
151  CHKERR mField.get_moab().get_entities_by_handle(it->meshset, fix_ents,
152  true);
153  }
154  }
155  return fix_ents;
156  };
157 
158  auto remove_ents = [&](const Range &&ents, const bool fix_x,
159  const bool fix_y) {
160  auto prb_mng = mField.getInterface<ProblemsManager>();
161  auto simple = mField.getInterface<Simple>();
163  Range verts;
164  CHKERR mField.get_moab().get_connectivity(ents, verts, true);
165  verts.merge(ents);
166  const int lo_coeff = fix_x ? 0 : 1;
167  const int hi_coeff = fix_y ? 1 : 0;
168  CHKERR prb_mng->removeDofsOnEntities(simple->getProblemName(), "U", verts,
169  lo_coeff, hi_coeff);
171  };
172 
173  CHKERR remove_ents(fix_disp("FIX_X"), true, false);
174  CHKERR remove_ents(fix_disp("FIX_Y"), false, true);
175  CHKERR remove_ents(fix_disp("FIX_ALL"), true, true);
176 
178 }
179 //! [Boundary condition]
180 
181 //! [Push operators to pipeline]
184  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
185 
186  auto add_domain_base_ops = [&](auto &pipeline) {
187  pipeline.push_back(new OpCalculateInvJacForFace(invJac));
188  pipeline.push_back(new OpSetInvJacH1ForFace(invJac));
189 
190  pipeline.push_back(
191  new OpCalculateVectorFieldGradient<2, 2>("U", commonDataPtr->mGradPtr));
192  pipeline.push_back(new OpStrain("U", commonDataPtr));
193 
194  pipeline.push_back(new OpCalculateScalarFieldValues(
195  "TAU", commonDataPtr->plasticTauPtr, MBTRI));
196  pipeline.push_back(
198  "TAU", commonDataPtr->plasticTauDotPtr, MBTRI));
199 
200  pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<2>(
201  "EP", commonDataPtr->plasticStrainPtr, MBTRI));
202  pipeline.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<2>(
203  "EP", commonDataPtr->plasticStrainDotPtr, MBTRI));
204 
205  pipeline.push_back(new OpPlasticStress("U", commonDataPtr));
206  pipeline.push_back(new OpCalculatePlasticSurface("U", commonDataPtr));
207  };
208 
209  auto add_domain_ops_lhs = [&](auto &pipeline) {
210  pipeline.push_back(new OpStiffnessMatrixLhs("U", "U", commonDataPtr));
211  pipeline.push_back(
212  new OpCalculatePlasticInternalForceLhs_dEP("U", "EP", commonDataPtr));
213 
214  pipeline.push_back(
215  new OpCalculatePlasticFlowLhs_dU("EP", "U", commonDataPtr));
216  pipeline.push_back(
217  new OpCalculatePlasticFlowLhs_dEP("EP", "EP", commonDataPtr));
218  pipeline.push_back(
219  new OpCalculatePlasticFlowLhs_dTAU("EP", "TAU", commonDataPtr));
220 
221  pipeline.push_back(
222  new OpCalculateContrainsLhs_dU("TAU", "U", commonDataPtr));
223  pipeline.push_back(
224  new OpCalculateContrainsLhs_dEP("TAU", "EP", commonDataPtr));
225  pipeline.push_back(
226  new OpCalculateContrainsLhs_dTAU("TAU", "TAU", commonDataPtr));
227  };
228 
229  auto add_domain_ops_rhs = [&](auto &pipeline) {
230  auto gravity = [](double x, double y) {
231  return FTensor::Tensor1<double, 2>{0., 1.};
232  };
233  pipeline.push_back(new OpForceRhs("U", commonDataPtr, gravity));
234  pipeline.push_back(new OpCalculatePlasticFlowRhs("EP", commonDataPtr));
235  pipeline.push_back(new OpCalculateContrainsRhs("TAU", commonDataPtr));
236  pipeline.push_back(new OpInternalForceRhs("U", commonDataPtr));
237  };
238 
239  auto add_boundary_ops_rhs = [&](auto &pipeline) {
241 
243  if (it->getName().compare(0, 5, "FORCE") == 0) {
244  Range my_edges;
245  std::vector<double> force_vec;
246  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
247  my_edges, true);
248  it->getAttributes(force_vec);
249  pipeline.push_back(new OpEdgeForceRhs("U", my_edges, force_vec));
250  }
251  }
252 
254  };
255 
256  add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
257  add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
258  add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
259  add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
260  add_boundary_ops_rhs(pipeline_mng->getOpBoundaryRhsPipeline());
261 
262  auto integration_rule = [](int, int, int approx_order) {
263  return 2 * approx_order;
264  };
265  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
266  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
267  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
268 
270 }
271 //! [Push operators to pipeline]
272 
273 //! [Solve]
276 
277  Simple *simple = mField.getInterface<Simple>();
278  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
279  ISManager *is_manager = mField.getInterface<ISManager>();
280 
281  auto solver = pipeline_mng->createTS();
282 
283  auto dm = simple->getDM();
284  auto D = smartCreateDMVector(dm);
285 
286  CHKERR TSSetSolution(solver, D);
287  CHKERR TSSetFromOptions(solver);
288  CHKERR TSSetUp(solver);
289 
290  auto set_section_monitor = [&]() {
292  SNES snes;
293  CHKERR TSGetSNES(solver, &snes);
294  PetscViewerAndFormat *vf;
295  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
296  PETSC_VIEWER_DEFAULT, &vf);
297  CHKERR SNESMonitorSet(
298  snes,
299  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
300  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
302  };
303 
304  auto create_post_process_element = [&]() {
306  postProcFe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
307  postProcFe->generateReferenceElementMesh();
308 
309  postProcFe->getOpPtrVector().push_back(
311  postProcFe->getOpPtrVector().push_back(new OpSetInvJacH1ForFace(invJac));
312  postProcFe->getOpPtrVector().push_back(
313  new OpCalculateVectorFieldGradient<2, 2>("U", commonDataPtr->mGradPtr));
314  postProcFe->getOpPtrVector().push_back(new OpStrain("U", commonDataPtr));
315 
316  postProcFe->getOpPtrVector().push_back(
317  new OpCalculateVectorFieldGradient<2, 2>("U", commonDataPtr->mGradPtr));
318  postProcFe->getOpPtrVector().push_back(new OpStrain("U", commonDataPtr));
319  postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
320  "TAU", commonDataPtr->plasticTauPtr, MBTRI));
321  postProcFe->getOpPtrVector().push_back(
323  "EP", commonDataPtr->plasticStrainPtr, MBTRI));
324  postProcFe->getOpPtrVector().push_back(
325  new OpPlasticStress("U", commonDataPtr));
326  postProcFe->getOpPtrVector().push_back(
327  new OpCalculatePlasticSurface("U", commonDataPtr));
328 
329  postProcFe->getOpPtrVector().push_back(new OpPostProcElastic(
330  "U", postProcFe->postProcMesh, postProcFe->mapGaussPts, commonDataPtr));
331  postProcFe->getOpPtrVector().push_back(new OpPostProcPlastic(
332  "U", postProcFe->postProcMesh, postProcFe->mapGaussPts, commonDataPtr));
333  postProcFe->addFieldValuesPostProc("U");
335  };
336 
337  auto scatter_create = [&](auto coeff) {
339  CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
340  ROW, "U", coeff, coeff, is);
341  int loc_size;
342  CHKERR ISGetLocalSize(is, &loc_size);
343  Vec v;
344  CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
345  VecScatter scatter;
346  CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
347  return std::make_tuple(SmartPetscObj<Vec>(v),
348  SmartPetscObj<VecScatter>(scatter));
349  };
350 
351  auto set_time_monitor = [&]() {
353  boost::shared_ptr<Monitor> monitor_ptr(
354  new Monitor(dm, postProcFe, uXScatter, uYScatter));
355  boost::shared_ptr<ForcesAndSourcesCore> null;
356  CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
357  monitor_ptr, null, null);
359  };
360 
361  CHKERR set_section_monitor();
362  CHKERR create_post_process_element();
363  uXScatter = scatter_create(0);
364  uYScatter = scatter_create(1);
365  CHKERR set_time_monitor();
366 
367  CHKERR TSSolve(solver, D);
368 
369  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
370  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
371  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
372 
374 }
375 //! [Solve]
376 
377 //! [Postprocess results]
380 
381 
383 }
384 //! [Postprocess results]
385 
386 //! [Check]
390 }
391 //! [Check]
392 
393 static char help[] = "...\n\n";
394 
395 int main(int argc, char *argv[]) {
396 
397  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
398 
399  try {
400 
401  //! [Register MoFEM discrete manager in PETSc]
402  DMType dm_name = "DMMOFEM";
403  CHKERR DMRegister_MoFEM(dm_name);
404  //! [Register MoFEM discrete manager in PETSc
405 
406  //! [Create MoAB]
407  moab::Core mb_instance; ///< mesh database
408  moab::Interface &moab = mb_instance; ///< mesh database interface
409  //! [Create MoAB]
410 
411  //! [Create MoFEM]
412  MoFEM::Core core(moab); ///< finite element database
413  MoFEM::Interface &m_field = core; ///< finite element database insterface
414  //! [Create MoFEM]
415 
416  //! [Load mesh]
417  Simple *simple = m_field.getInterface<Simple>();
418  CHKERR simple->getOptions();
419  CHKERR simple->loadFile("");
420  //! [Load mesh]
421 
422  //! [Example]
423  Example ex(m_field);
424  CHKERR ex.runProblem();
425  //! [Example]
426  }
427  CATCH_ERRORS;
428 
430 }
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
int main(int argc, char *argv[])
Deprecated interface functions.
Problem manager is used to build and partition problems.
MoFEMErrorCode createCommonData()
[Set up problem]
Section manager is used to create indexes and sectionsFIXME: ISManager is not properly testsed by ato...
Definition: ISManager.hpp:36
constexpr double poisson_ratio
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:445
static constexpr int approx_order
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
MatrixDouble invJac
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
Example(MoFEM::Interface &m_field)
Core (interface) class.
Definition: Core.hpp:50
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProcFe
constexpr double H
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
DataForcesAndSourcesCore::EntData EntData
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:939
MoFEMErrorCode runProblem()
[Operator]
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEMErrorCode OPs()
[Boundary condition]
constexpr int order
MoFEMErrorCode checkResults()
[Print results]
Calculate inverse of jacobian for face element.
MoFEMErrorCode bC()
[Create common data]
Definition: PlasticOps.hpp:203
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
static char help[]
[Check]
constexpr double cn
constexpr double young_modulus
PipelineManager interface.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
#define CHKERR
Inline error check.
Definition: definitions.h:602
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMMoFEM.cpp:915
MoFEMErrorCode setUP()
[Set up problem]
boost::shared_ptr< OpPlasticTools::CommonData > commonDataPtr
Transform local reference derivatives of shape functions to global derivatives.
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
MoFEMErrorCode postProcess()
[Integrate]
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Definition: ElasticOps.hpp:86
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Get rate of scalar field at integration points.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:439
continuous field
Definition: definitions.h:177
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Calculate symmetric tensor field rates ant integratio pts.
intrusive_ptr for managing petsc objects
Definition: AuxPETSc.hpp:128
constexpr double sigmaY
Calculate symmetric tensor field values at integration pts.
Get value at integration points for scalar field.
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
FTensor::Index< 'l', 2 > l
Definition: ContactOps.hpp:29
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:69
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
field with C-1 continuity
Definition: definitions.h:180