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>>;
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); }
21constexpr int SPACE_DIM = 3;
22// constexpr auto A = AssemblyType::BLOCK_SCHUR;
23constexpr auto A = AssemblyType::BLOCK_MAT;
25#ifndef __ESHELBIAN_PLASTICITY_HPP__
26 #define __ESHELBIAN_PLASTICITY_HPP__
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 {
39 using CachePhi = boost::tuple<int, int, MatrixDouble>;
41 CGGUserPolynomialBase(boost::shared_ptr<CachePhi> cache_phi = nullptr);
44 MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
45 BaseFunctionUnknownInterface **iface) const;
46 MoFEMErrorCode getValue(MatrixDouble &pts,
47 boost::shared_ptr<BaseFunctionCtx> ctx_ptr);
49private:
50 MatrixDouble shapeFun;
51 boost::shared_ptr<CachePhi> cachePhiPtr;
53 MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts);
54};
56struct ContactTree;
57
62 LOG /*linar extenion*/,
63 LOG_QUADRATIC /*quadratic extension*/,
65};
67
68using MatrixPtr = boost::shared_ptr<MatrixDouble>;
69using VectorPtr = boost::shared_ptr<VectorDouble>;
71using EntData = EntitiesFieldData::EntData;
72using UserDataOperator = ForcesAndSourcesCore::UserDataOperator;
73using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator;
74using FaceUserDataOperator = FaceElementForcesAndSourcesCore::UserDataOperator;
75using EleOnSide = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
76using SideEleOp = EleOnSide::UserDataOperator;
77
79
80struct AnalyticalExprPython;
81
83 : public boost::enable_shared_from_this<DataAtIntegrationPts> {
84
85 MatrixDouble approxPAtPts;
86 MatrixDouble approxSigmaAtPts;
87 MatrixDouble divPAtPts;
88 MatrixDouble divSigmaAtPts;
89 MatrixDouble wL2AtPts;
90 MatrixDouble wL2DotAtPts;
91 MatrixDouble wL2DotDotAtPts;
93 MatrixDouble stretchTensorAtPts;
94 VectorDouble jacobianAtPts;
95 MatrixDouble tractionAtPts;
106 MatrixDouble WAtPts;
107 MatrixDouble W0AtPts;
108 MatrixDouble GAtPts;
109 MatrixDouble G0AtPts;
110 MatrixDouble wH1AtPts;
111 MatrixDouble XH1AtPts;
112 MatrixDouble contactL2AtPts;
113 MatrixDouble wGradH1AtPts;
114 MatrixDouble logStretch2H1AtPts;
125 MatrixDouble adjointPdUAtPts;
127 MatrixDouble adjointPdUdPAtPts;
128 MatrixDouble rotMatAtPts;
129 MatrixDouble PAtPts;
130 MatrixDouble SigmaAtPts;
131 VectorDouble energyAtPts; //< this is density of energy at integration points
132 MatrixDouble flowL2AtPts;
133
134 MatrixDouble varRotAxis;
135 MatrixDouble varLogStreach;
136 MatrixDouble varPiola;
137 MatrixDouble varDivPiola;
138 MatrixDouble varWL2;
140 MatrixDouble P_du;
142 MatrixDouble eigenVals;
143 MatrixDouble eigenVecs;
144 VectorInt nbUniq;
145 MatrixDouble eigenValsC;
146 MatrixDouble eigenVecsC;
147 VectorInt nbUniqC;
148
149 MatrixDouble matD;
150 MatrixDouble matAxiatorD;
151 MatrixDouble matDeviatorD;
152 MatrixDouble matInvD;
153
156 MatrixDouble facePiolaAtPts;
161 double mu;
162 double lambda;
163 double piolaScale = 1.;
164 double faceEnergy = 0.;
166 inline auto getPiolaScalePtr() {
167 return boost::shared_ptr<double>(shared_from_this(), &piolaScale);
168 }
169
171 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
173 }
175 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &approxPAtPts);
177
179 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &divPAtPts);
183 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &divSigmaAtPts);
187 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wL2AtPts);
191 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wL2DotAtPts);
192 }
195 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wL2DotDotAtPts);
196 }
199 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
204 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
209 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
211 }
212
214 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
216 }
219 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &rotAxisAtPts);
223 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &rotAxis0AtPts);
224 }
227 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
229 }
232 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
235
237 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &GAtPts);
241 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &GAtPts);
242 }
243
245 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matD);
246 }
247
249 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matInvD);
250 }
251
253 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matAxiatorD);
254 }
255
257 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matDeviatorD);
258 }
259
261 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wH1AtPts);
262 }
263
265 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &XH1AtPts);
266 }
267
269 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactL2AtPts);
270 }
271
273 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wGradH1AtPts);
274 }
275
277 return boost::shared_ptr<VectorDouble>(shared_from_this(), &jacobianAtPts);
278 }
279
281 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
286 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varRotAxis);
287 }
288
290 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varLogStreach);
291 }
292
294 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varPiola);
296
298 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varDivPiola);
299 }
300
302 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varWL2);
306 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
309
311 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &facePiolaAtPts);
312 }
313
315 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
317 }
318
320 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &eigenVals);
321 }
322
323 boost::shared_ptr<PhysicalEquations> physicsPtr;
326struct OpJacobian;
327
328// Forward declarations
329struct ExternalStrain;
330using ExternalStrainVec = std::vector<ExternalStrain>;
332
339
344
346 PhysicalEquations(const int size_active, const int size_dependent)
347 : activeVariables(size_active, 0),
348 dependentVariablesPiola(size_dependent, 0),
349 dependentVariablesPiolaDirevatives(size_dependent * size_active, 0) {}
350 virtual ~PhysicalEquations() = default;
352 virtual MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h) = 0;
353
354 virtual UserDataOperator *
355 returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
356 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
357 boost::shared_ptr<PhysicalEquations> physics_ptr);
358
359 virtual VolUserDataOperator *
360 returnOpSpatialPhysical(const std::string &field_name,
361 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
362 const double alpha_u);
366 const std::string &field_name,
367 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
368 boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
369 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
370
372 std::string row_field, std::string col_field,
373 boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha);
374
375 virtual VolUserDataOperator *
376 returnOpCalculateEnergy(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
377 boost::shared_ptr<double> total_energy_ptr);
378
380 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
381 boost::shared_ptr<PhysicalEquations> physics_ptr);
382
384 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
385 boost::shared_ptr<PhysicalEquations> physics_ptr);
387 virtual VolUserDataOperator *
388 returnOpSetScale(boost::shared_ptr<double> scale_ptr,
389 boost::shared_ptr<PhysicalEquations> physics_ptr);
390
391 std::vector<double> activeVariables;
392 std::vector<double> dependentVariablesPiola;
395 /** \name Active variables */
396
397 /**@{*/
398
399 template <int S>
400 inline static DTensor2Ptr get_VecTensor2(std::vector<double> &v) {
401 return DTensor2Ptr(&v[S + 0], &v[S + 1], &v[S + 2], &v[S + 3], &v[S + 4],
402 &v[S + 5], &v[S + 6], &v[S + 7], &v[S + 8]);
405 template <int S>
406 inline static DTensor0Ptr get_VecTensor0(std::vector<double> &v) {
407 return DTensor0Ptr(&v[S + 0]);
410 template <int S0>
411 inline static DTensor3Ptr get_vecTensor3(std::vector<double> &v,
412 const int nba) {
414 const int A00 = nba * 0 + S0;
415 const int A01 = nba * 1 + S0;
416 const int A02 = nba * 2 + S0;
417 const int A10 = nba * 3 + S0;
418 const int A11 = nba * 4 + S0;
419 const int A12 = nba * 5 + S0;
420 const int A20 = nba * 6 + S0;
421 const int A21 = nba * 7 + S0;
422 const int A22 = nba * 8 + S0;
424 return DTensor3Ptr(
426 &v[A00 + 0], &v[A00 + 1], &v[A00 + 2], &v[A01 + 0], &v[A01 + 1],
427 &v[A01 + 2], &v[A02 + 0], &v[A02 + 1], &v[A02 + 2],
429 &v[A10 + 0], &v[A10 + 1], &v[A10 + 2], &v[A11 + 0], &v[A11 + 1],
430 &v[A11 + 2], &v[A12 + 0], &v[A12 + 1], &v[A12 + 2],
432 &v[A20 + 0], &v[A20 + 1], &v[A20 + 2], &v[A21 + 0], &v[A21 + 1],
433 &v[A21 + 2], &v[A22 + 0], &v[A22 + 1], &v[A22 + 2]
435 );
438 /**@}*/
440 /** \name Active variables */
442 /**@{*/
443
446 /**@}*/
448 /** \name Dependent variables */
450 /**@{*/
455
456 /**@}*/
458 /** \name Derivatives of dependent variables */
459
460 /**@{*/
461
474
475 /**@}*/
476};
477
478struct BcDisp {
479 BcDisp(std::string name, std::vector<double> attr, Range faces);
480 std::string blockName;
482 VectorDouble3 vals;
483 VectorInt3 flags;
484};
485using BcDispVec = std::vector<BcDisp>;
486
487struct BcRot {
488 BcRot(std::string name, std::vector<double> attr, Range faces);
489 std::string blockName;
491 VectorDouble vals;
492 double theta;
493};
494using BcRotVec = std::vector<BcRot>;
495
496typedef std::vector<Range> TractionFreeBc;
497
499 TractionBc(std::string name, std::vector<double> attr, Range faces);
500 std::string blockName;
502 VectorDouble3 vals;
503 VectorInt3 flags;
505using TractionBcVec = std::vector<TractionBc>;
506
508 NormalDisplacementBc(std::string name, std::vector<double> attr, Range faces);
509 std::string blockName;
511 double val;
512};
513using NormalDisplacementBcVec = std::vector<NormalDisplacementBc>;
514
516 AnalyticalDisplacementBc(std::string name, std::vector<double> attr,
517 Range faces);
518 std::string blockName;
520 VectorInt3 flags;
521};
522using AnalyticalDisplacementBcVec = std::vector<AnalyticalDisplacementBc>;
525 AnalyticalTractionBc(std::string name, std::vector<double> attr, Range faces);
526 std::string blockName;
528 VectorInt3 flags;
529};
530using AnalyticalTractionBcVec = std::vector<AnalyticalTractionBc>;
533 PressureBc(std::string name, std::vector<double> attr, Range faces);
534 std::string blockName;
536 double val;
538using PressureBcVec = std::vector<PressureBc>;
539
541 ExternalStrain(std::string name, std::vector<double> attr, Range ents);
542 std::string blockName;
544 double val;
548 #ifdef ENABLE_PYTHON_BINDING
549struct AnalyticalExprPython {
550 AnalyticalExprPython() = default;
551 virtual ~AnalyticalExprPython() = default;
553 MoFEMErrorCode analyticalExprInit(const std::string py_file);
554 MoFEMErrorCode evalAnalyticalDisp(double delta_t, double t, np::ndarray x,
555 np::ndarray y, np::ndarray z,
556 np::ndarray nx, np::ndarray ny,
557 np::ndarray nz,
558 const std::string &block_name,
559 np::ndarray &analytical_expr);
560 MoFEMErrorCode evalAnalyticalTraction(double delta_t, double t, np::ndarray x,
561 np::ndarray y, np::ndarray z,
562 np::ndarray nx, np::ndarray ny,
563 np::ndarray nz,
564 const std::string &block_name,
565 np::ndarray &analytical_expr);
566
567 template <typename T>
568 inline std::vector<T>
569 py_list_to_std_vector(const boost::python::object &iterable) {
570 return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
571 boost::python::stl_input_iterator<T>());
574private:
575 bp::object mainNamespace;
576 bp::object analyticalDispFun;
577 bp::object analyticalTractionFun;
580extern boost::weak_ptr<AnalyticalExprPython> AnalyticalExprPythonWeakPtr;
581 #endif
583 #include "EshelbianCore.hpp"
584 #include "EshelbianOperators.hpp"
586} // namespace EshelbianPlasticity
587
588#endif //__ESHELBIAN_PLASTICITY_HPP__
constexpr int SPACE_DIM
ForcesAndSourcesCore::UserDataOperator UserDataOperator
const double v
phase velocity of light in medium (cm/ns)
boost::shared_ptr< VectorDouble > VectorPtr
FaceElementForcesAndSourcesCore::UserDataOperator FaceUserDataOperator
std::vector< AnalyticalTractionBc > AnalyticalTractionBcVec
std::vector< BcDisp > BcDispVec
std::vector< NormalDisplacementBc > NormalDisplacementBcVec
std::vector< BcRot > BcRotVec
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
std::vector< ExternalStrain > ExternalStrainVec
std::vector< AnalyticalDisplacementBc > AnalyticalDisplacementBcVec
std::vector< Range > TractionFreeBc
std::vector< PressureBc > PressureBcVec
std::vector< TractionBc > TractionBcVec
boost::shared_ptr< MatrixDouble > MatrixPtr
constexpr double t
plate stiffness
Definition plate.cpp:58
constexpr auto field_name
AnalyticalDisplacementBc(std::string name, std::vector< double > attr, Range faces)
AnalyticalTractionBc(std::string name, std::vector< double > attr, Range faces)
BcDisp(std::string name, std::vector< double > attr, Range faces)
BcRot(std::string name, std::vector< double > attr, Range faces)
CGGUserPolynomialBase(boost::shared_ptr< CachePhi > cache_phi=nullptr)
MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts)
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
boost::tuple< int, int, MatrixDouble > CachePhi
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, BaseFunctionUnknownInterface **iface) const
boost::shared_ptr< PhysicalEquations > physicsPtr
ExternalStrain(std::string name, std::vector< double > attr, Range ents)
NormalDisplacementBc(std::string name, std::vector< double > attr, Range faces)
FTensor::Tensor1< adouble, 3 > ATensor1
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)
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)
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::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > DTensor3Ptr
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::Tensor3< adouble, 3, 3, 3 > ATensor3
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, 3, 3 > DTensor2Ptr
FTensor::Tensor2< double, 3, 3 > DTensor2
FTensor::Tensor3< double, 3, 3, 3 > DTensor3
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > DTensor0Ptr
virtual MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h)=0
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)
PressureBc(std::string name, std::vector< double > attr, Range faces)
TractionBc(std::string name, std::vector< double > attr, Range faces)
Data on single entity (This is passed as argument to DataOperator::doWork)