v0.10.0
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 setupProblem();
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 setupProblem();
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]
MoFEMErrorCode setupProblem()
[Set up problem]
CoreTmp< 0 > Core
Definition: Core.hpp:1129
const double D
diffusivity
[Example]
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:509
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:76
MatrixDouble invJac
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
Example(MoFEM::Interface &m_field)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
static Index< 'l', 3 > l
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]
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()
[Run problem]
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()
[Postprocess results]
Calculate inverse of jacobian for face element.
MoFEMErrorCode bC()
[Create common data]
static Index< 'i', 3 > i
Definition: PlasticOps.hpp:203
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
static Index< 'j', 3 > j
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
static Index< 'o', 3 > o
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:604
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
boost::shared_ptr< OpPlasticTools::CommonData > commonDataPtr
Transform local reference derivatives of shape functions to global derivatives.
MoFEMErrorCode postProcess()
[Solve]
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
static Index< 'k', 3 > k
Definition: ElasticOps.hpp:86
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Get rate of scalar field at integration points.
Core (interface) class.
Definition: Core.hpp:77
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
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.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...
field with C-1 continuity
Definition: definitions.h:180