v0.15.0
Loading...
Searching...
No Matches
EshelbianPlasticity.hpp
Go to the documentation of this file.
1/**
2 * \file EshelbianPlasticity.hpp
3 * \brief Eshelbian plasticity interface
4 *
5 * \brief Problem implementation for mix element for large-strain elasticity
6 *
7 * For reference on mixed formulation see: \cite gopalakrishnan2012second and
8 * \cite cockburn2010new
9 *
10 * \todo Implementation of plasticity
11 */
12
13// DO NOT DELETE MIGHT BE USEFULL
14// #include <boost/multiprecision/mpfr.hpp>
15// using mpfr50 = boost::multiprecision::number<
16// boost::multiprecision::mpfr_float_backend<50>>;
17
18// inline mpfr50 exp_hi(const mpfr50 x) { return boost::multiprecision::exp(x); }
19// inline mpfr50 log_hi(const mpfr50 x) { return boost::multiprecision::log(x); }
20
21constexpr int SPACE_DIM = 3;
22// constexpr auto A = AssemblyType::BLOCK_SCHUR;
23constexpr auto A = AssemblyType::BLOCK_MAT;
24
25#ifndef __ESHELBIAN_PLASTICITY_HPP__
26 #define __ESHELBIAN_PLASTICITY_HPP__
27
28 #ifdef ENABLE_PYTHON_BINDING
29 #include <boost/python.hpp>
30 #include <boost/python/def.hpp>
31 #include <boost/python/numpy.hpp>
32namespace bp = boost::python;
33namespace np = boost::python::numpy;
34 #endif
35namespace EshelbianPlasticity {
36
37struct ContactTree;
38
43 LOG /*linar extenion*/,
44 LOG_QUADRATIC /*quadratic extension*/,
46};
48
49using MatrixPtr = boost::shared_ptr<MatrixDouble>;
50using VectorPtr = boost::shared_ptr<VectorDouble>;
51
52using EntData = EntitiesFieldData::EntData;
53using UserDataOperator = ForcesAndSourcesCore::UserDataOperator;
54using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator;
55using FaceUserDataOperator = FaceElementForcesAndSourcesCore::UserDataOperator;
56using EleOnSide = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
57using SideEleOp = EleOnSide::UserDataOperator;
58
60
61struct AnalyticalExprPython;
62
64 : public boost::enable_shared_from_this<DataAtIntegrationPts> {
65
66 MatrixDouble approxPAtPts;
67 MatrixDouble approxSigmaAtPts;
68 MatrixDouble divPAtPts;
69 MatrixDouble gradPAtPts;
70 MatrixDouble divSigmaAtPts;
71 MatrixDouble wL2AtPts;
72 MatrixDouble wL2DotAtPts;
73 MatrixDouble wL2DotDotAtPts;
75 MatrixDouble stretchTensorAtPts;
76 VectorDouble jacobianAtPts;
77 MatrixDouble tractionAtPts;
78
84 MatrixDouble rotAxisAtPts;
85 MatrixDouble rotAxisDotAtPts;
86 MatrixDouble rotAxisGradDotAtPts;
87 MatrixDouble rotAxis0AtPts;
88 MatrixDouble WAtPts;
89 MatrixDouble W0AtPts;
90 MatrixDouble GAtPts;
91 MatrixDouble G0AtPts;
92 MatrixDouble wH1AtPts;
93 MatrixDouble XH1AtPts;
94 MatrixDouble contactL2AtPts;
95 MatrixDouble wGradH1AtPts;
96 MatrixDouble logStretch2H1AtPts;
98
99 MatrixDouble hAtPts;
100 MatrixDouble hdOmegaAtPts;
101 MatrixDouble hdLogStretchAtPts;
102 MatrixDouble leviKirchhoffAtPts;
107 MatrixDouble adjointPdUAtPts;
109 MatrixDouble adjointPdUdPAtPts;
110 MatrixDouble rotMatAtPts;
111 MatrixDouble PAtPts;
112 MatrixDouble SigmaAtPts;
113 VectorDouble energyAtPts; //< this is density of energy at integration points
114 MatrixDouble flowL2AtPts;
115
116 MatrixDouble varRotAxis;
117 MatrixDouble varLogStreach;
118 MatrixDouble varPiola;
119 MatrixDouble varDivPiola;
120 MatrixDouble varWL2;
121
122 MatrixDouble P_du;
123
124 MatrixDouble eigenVals;
125 MatrixDouble eigenVecs;
126 VectorInt nbUniq;
127 MatrixDouble eigenValsC;
128 MatrixDouble eigenVecsC;
129 VectorInt nbUniqC;
130
131 MatrixDouble matD;
132 MatrixDouble matAxiatorD;
133 MatrixDouble matDeviatorD;
134 MatrixDouble matInvD;
135
137
138 MatrixDouble facePiolaAtPts;
139 MatrixDouble hybridDispAtPts;
143
144 double mu;
145 double lambda;
146 double piolaScale = 1.;
147 double faceEnergy = 0.;
148
149 inline auto getPiolaScalePtr() {
150 return boost::shared_ptr<double>(shared_from_this(), &piolaScale);
151 }
152
154 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
156 }
158 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &approxPAtPts);
159 }
160
162 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &divPAtPts);
163 }
164
166 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &gradPAtPts);
167 }
168
170 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &divSigmaAtPts);
171 }
172
174 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wL2AtPts);
175 }
176
178 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wL2DotAtPts);
179 }
180
182 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wL2DotDotAtPts);
183 }
184
186 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
188 }
189
191 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
193 }
194
196 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
198 }
199
201 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
203 }
204
206 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &rotAxisAtPts);
207 }
208
210 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &rotAxis0AtPts);
211 }
212
214 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
216 }
217
219 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
221 }
222
224 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &GAtPts);
225 }
226
228 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &GAtPts);
229 }
230
232 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matD);
233 }
234
236 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matInvD);
237 }
238
240 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matAxiatorD);
241 }
242
244 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matDeviatorD);
245 }
246
248 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wH1AtPts);
249 }
250
252 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &XH1AtPts);
253 }
254
256 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactL2AtPts);
257 }
258
260 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wGradH1AtPts);
261 }
262
264 return boost::shared_ptr<VectorDouble>(shared_from_this(), &jacobianAtPts);
265 }
266
268 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
270 };
271
273 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varRotAxis);
274 }
275
277 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varLogStreach);
278 }
279
281 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varPiola);
282 }
283
285 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varDivPiola);
286 }
287
289 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varWL2);
290 }
291
293 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
295 }
296
298 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &facePiolaAtPts);
299 }
300
302 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
304 }
305
307 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
309 }
310
312 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &eigenVals);
313 }
314
315 boost::shared_ptr<PhysicalEquations> physicsPtr;
316};
317
318struct OpJacobian;
319
320// Forward declarations
321struct ExternalStrain;
322using ExternalStrainVec = std::vector<ExternalStrain>;
324
331
333
336
338 PhysicalEquations(const int size_active, const int size_dependent)
339 : activeVariables(size_active, 0),
340 dependentVariablesPiola(size_dependent, 0),
341 dependentVariablesPiolaDirevatives(size_dependent * size_active, 0) {}
342 virtual ~PhysicalEquations() = default;
343
344 virtual MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h) = 0;
345
346 virtual UserDataOperator *
347 returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
348 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
349 boost::shared_ptr<PhysicalEquations> physics_ptr);
350
351 virtual VolUserDataOperator *
352 returnOpSpatialPhysical(const std::string &field_name,
353 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
354 const double alpha_u);
355
356 virtual VolUserDataOperator *
358 const std::string &field_name,
359 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
360 boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
361 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
362
364 std::string row_field, std::string col_field,
365 boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha);
366
367 virtual VolUserDataOperator *
368 returnOpCalculateEnergy(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
369 boost::shared_ptr<double> total_energy_ptr);
370
372 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
373 boost::shared_ptr<PhysicalEquations> physics_ptr);
374
376 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
377 boost::shared_ptr<PhysicalEquations> physics_ptr);
378
379 virtual VolUserDataOperator *
380 returnOpSetScale(boost::shared_ptr<double> scale_ptr,
381 boost::shared_ptr<PhysicalEquations> physics_ptr);
382
383 std::vector<double> activeVariables;
384 std::vector<double> dependentVariablesPiola;
386
387 /** \name Active variables */
388
389 /**@{*/
390
391 template <int S>
392 inline static DTensor2Ptr get_VecTensor2(std::vector<double> &v) {
393 return DTensor2Ptr(&v[S + 0], &v[S + 1], &v[S + 2], &v[S + 3], &v[S + 4],
394 &v[S + 5], &v[S + 6], &v[S + 7], &v[S + 8]);
395 }
396
397 template <int S>
398 inline static DTensor0Ptr get_VecTensor0(std::vector<double> &v) {
399 return DTensor0Ptr(&v[S + 0]);
400 }
401
402 template <int S0>
403 inline static DTensor3Ptr get_vecTensor3(std::vector<double> &v,
404 const int nba) {
405
406 const int A00 = nba * 0 + S0;
407 const int A01 = nba * 1 + S0;
408 const int A02 = nba * 2 + S0;
409 const int A10 = nba * 3 + S0;
410 const int A11 = nba * 4 + S0;
411 const int A12 = nba * 5 + S0;
412 const int A20 = nba * 6 + S0;
413 const int A21 = nba * 7 + S0;
414 const int A22 = nba * 8 + S0;
415
416 return DTensor3Ptr(
417
418 &v[A00 + 0], &v[A00 + 1], &v[A00 + 2], &v[A01 + 0], &v[A01 + 1],
419 &v[A01 + 2], &v[A02 + 0], &v[A02 + 1], &v[A02 + 2],
420
421 &v[A10 + 0], &v[A10 + 1], &v[A10 + 2], &v[A11 + 0], &v[A11 + 1],
422 &v[A11 + 2], &v[A12 + 0], &v[A12 + 1], &v[A12 + 2],
423
424 &v[A20 + 0], &v[A20 + 1], &v[A20 + 2], &v[A21 + 0], &v[A21 + 1],
425 &v[A21 + 2], &v[A22 + 0], &v[A22 + 1], &v[A22 + 2]
426
427 );
428 }
429
430 /**@}*/
431
432 /** \name Active variables */
433
434 /**@{*/
435
436 inline DTensor2Ptr get_h() { return get_VecTensor2<0>(activeVariables); }
437
438 /**@}*/
439
440 /** \name Dependent variables */
441
442 /**@{*/
443
445 return get_VecTensor2<0>(dependentVariablesPiola);
446 }
447
448 /**@}*/
449
450 /** \name Derivatives of dependent variables */
451
452 /**@{*/
453
455 return get_vecTensor3<0>(dependentVariablesPiolaDirevatives,
456 activeVariables.size());
457 }
459 return get_vecTensor3<3>(dependentVariablesPiolaDirevatives,
460 activeVariables.size());
461 }
463 return get_vecTensor3<6>(dependentVariablesPiolaDirevatives,
464 activeVariables.size());
465 }
466
467 /**@}*/
468};
469
470struct BcDisp {
471 BcDisp(std::string name, std::vector<double> attr, Range faces);
472 std::string blockName;
474 VectorDouble3 vals;
475 VectorInt3 flags;
476};
477using BcDispVec = std::vector<BcDisp>;
478
479struct BcRot {
480 BcRot(std::string name, std::vector<double> attr, Range faces);
481 std::string blockName;
483 VectorDouble vals;
484 double theta;
485};
486using BcRotVec = std::vector<BcRot>;
487
488typedef std::vector<Range> TractionFreeBc;
489
491 TractionBc(std::string name, std::vector<double> attr, Range faces);
492 std::string blockName;
494 VectorDouble3 vals;
495 VectorInt3 flags;
496};
497using TractionBcVec = std::vector<TractionBc>;
498
500 NormalDisplacementBc(std::string name, std::vector<double> attr, Range faces);
501 std::string blockName;
503 double val;
504};
505using NormalDisplacementBcVec = std::vector<NormalDisplacementBc>;
506
508 AnalyticalDisplacementBc(std::string name, std::vector<double> attr,
509 Range faces);
510 std::string blockName;
512 VectorInt3 flags;
513};
514using AnalyticalDisplacementBcVec = std::vector<AnalyticalDisplacementBc>;
515
517 AnalyticalTractionBc(std::string name, std::vector<double> attr, Range faces);
518 std::string blockName;
520 VectorInt3 flags;
521};
522using AnalyticalTractionBcVec = std::vector<AnalyticalTractionBc>;
523
525 PressureBc(std::string name, std::vector<double> attr, Range faces);
526 std::string blockName;
528 double val;
529};
530using PressureBcVec = std::vector<PressureBc>;
531
533 ExternalStrain(std::string name, std::vector<double> attr, Range ents);
534 std::string blockName;
536 double val;
538};
539
540 #ifdef ENABLE_PYTHON_BINDING
541struct AnalyticalExprPython {
542 AnalyticalExprPython() = default;
543 virtual ~AnalyticalExprPython() = default;
544
545 MoFEMErrorCode analyticalExprInit(const std::string py_file);
546 MoFEMErrorCode evalAnalyticalDisp(double delta_t, double t, np::ndarray x,
547 np::ndarray y, np::ndarray z,
548 np::ndarray nx, np::ndarray ny,
549 np::ndarray nz,
550 const std::string &block_name,
551 np::ndarray &analytical_expr);
552 MoFEMErrorCode evalAnalyticalTraction(double delta_t, double t, np::ndarray x,
553 np::ndarray y, np::ndarray z,
554 np::ndarray nx, np::ndarray ny,
555 np::ndarray nz,
556 const std::string &block_name,
557 np::ndarray &analytical_expr);
558
559 template <typename T>
560 inline std::vector<T>
561 py_list_to_std_vector(const boost::python::object &iterable) {
562 return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
563 boost::python::stl_input_iterator<T>());
564 }
565
566private:
567 bp::object mainNamespace;
568 bp::object analyticalDispFun;
569 bp::object analyticalTractionFun;
570};
571
572extern boost::weak_ptr<AnalyticalExprPython> AnalyticalExprPythonWeakPtr;
573 #endif
574
576 #include "EshelbianCore.hpp"
577 #include "EshelbianOperators.hpp"
578
579} // namespace EshelbianPlasticity
580
581#endif //__ESHELBIAN_PLASTICITY_HPP__
constexpr int SPACE_DIM
constexpr auto A
const double v
phase velocity of light in medium (cm/ns)
boost::shared_ptr< MatrixDouble > MatrixPtr
std::vector< AnalyticalTractionBc > AnalyticalTractionBcVec
std::vector< TractionBc > TractionBcVec
std::vector< AnalyticalDisplacementBc > AnalyticalDisplacementBcVec
std::vector< PressureBc > PressureBcVec
std::vector< Range > TractionFreeBc
std::vector< ExternalStrain > ExternalStrainVec
ForcesAndSourcesCore::UserDataOperator UserDataOperator
std::vector< BcRot > BcRotVec
std::vector< NormalDisplacementBc > NormalDisplacementBcVec
FaceElementForcesAndSourcesCore::UserDataOperator FaceUserDataOperator
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
std::vector< BcDisp > BcDispVec
boost::shared_ptr< VectorDouble > VectorPtr
constexpr double t
plate stiffness
Definition plate.cpp:58
constexpr auto field_name
boost::shared_ptr< PhysicalEquations > physicsPtr
virtual VolUserDataOperator * returnOpCalculateVarStretchFromStress(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
virtual VolUserDataOperator * returnOpCalculateEnergy(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< double > total_energy_ptr)
FTensor::Tensor3< double, 3, 3, 3 > DTensor3
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > DTensor3Ptr
virtual VolUserDataOperator * returnOpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
virtual VolUserDataOperator * returnOpSpatialPhysicalExternalStrain(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< ExternalStrainVec > external_strain_vec_ptr, std::map< std::string, boost::shared_ptr< ScalingMethod > > smv)
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, 3, 3 > DTensor2Ptr
virtual VolUserDataOperator * returnOpCalculateStretchFromStress(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
static DTensor3Ptr get_vecTensor3(std::vector< double > &v, const int nba)
FTensor::Tensor2< adouble, 3, 3 > ATensor2
FTensor::Tensor2< double, 3, 3 > DTensor2
PhysicalEquations(const int size_active, const int size_dependent)
virtual VolUserDataOperator * returnOpSetScale(boost::shared_ptr< double > scale_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
static DTensor2Ptr get_VecTensor2(std::vector< double > &v)
FTensor::Tensor1< adouble, 3 > ATensor1
FTensor::Tensor3< adouble, 3, 3, 3 > ATensor3
virtual MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h)=0
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > DTensor0Ptr
virtual VolUserDataOperator * returnOpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
virtual UserDataOperator * returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
static DTensor0Ptr get_VecTensor0(std::vector< double > &v)
Data on single entity (This is passed as argument to DataOperator::doWork)