v0.14.0
Loading...
Searching...
No Matches
HookeElement.hpp
Go to the documentation of this file.
1/**
2 * \file HookeElement.hpp
3 * \example HookeElement.hpp
4 *
5 * \brief Operators and data structures for linear elastic
6 * analysis
7 *
8 * Implemention of operators for Hooke material. Implementation is extended to
9 * the case when the mesh is moving as results of topological changes, also the
10 * calculation of material forces and associated tangent matrices are added to
11 * implementation.
12 *
13 * In other words, spatial deformation is small but topological changes large.
14 */
15
16
17
18#ifndef __HOOKE_ELEMENT_HPP
19#define __HOOKE_ELEMENT_HPP
20
21#ifndef __BASICFINITEELEMENTS_HPP__
23#endif // __BASICFINITEELEMENTS_HPP__
24
25#ifndef __NONLINEAR_ELASTIC_HPP
26
28
29 /** \brief data for calculation heat conductivity and heat capacity elements
30 * \ingroup nonlinear_elastic_elem
31 */
32 struct BlockData {
33 int iD;
34 double E;
36 Range tEts; ///< constrains elements in block set
39 };
40};
41
42#endif // __NONLINEAR_ELASTIC_HPP
43
44/** \brief structure grouping operators and data used for calculation of
45 * nonlinear elastic element \ingroup nonlinear_elastic_elem
46 *
47 * In order to assemble matrices and right hand vectors, the loops over
48 * elements, entities over that elements and finally loop over integration
49 * points are executed.
50 *
51 * Following implementation separate those three categories of loops and to each
52 * loop attach operator.
53 *
54 */
55
56#ifndef __CONVECTIVE_MASS_ELEMENT_HPP
58 /** \brief data for calculation inertia forces
59 * \ingroup user_modules
60 */
61 struct BlockData {
62 double rho0; ///< reference density
63 VectorDouble a0; ///< constant acceleration
64 Range tEts; ///< elements in block set
65 };
66}
67
68#endif //__CONVECTIVE_MASS_ELEMENT_HPP
69struct HookeElement {
70
73
75 using UserDataOperator = ForcesAndSourcesCore::UserDataOperator;
77 VolumeElementForcesAndSourcesCore::UserDataOperator;
78
80
81 boost::shared_ptr<MatrixDouble> smallStrainMat;
82 boost::shared_ptr<MatrixDouble> hMat;
83 boost::shared_ptr<MatrixDouble> FMat;
84
85 boost::shared_ptr<MatrixDouble> HMat;
86 boost::shared_ptr<VectorDouble> detHVec;
87 boost::shared_ptr<MatrixDouble> invHMat;
88
89 boost::shared_ptr<MatrixDouble> cauchyStressMat;
90 boost::shared_ptr<MatrixDouble> stiffnessMat;
91 boost::shared_ptr<VectorDouble> energyVec;
92 boost::shared_ptr<MatrixDouble> eshelbyStressMat;
93
94 boost::shared_ptr<MatrixDouble> eshelbyStress_dx;
95
97
98 smallStrainMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
99 hMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
100 FMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
101
102 HMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
103 detHVec = boost::shared_ptr<VectorDouble>(new VectorDouble());
104 invHMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
105
106 cauchyStressMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
107 stiffnessMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
108 energyVec = boost::shared_ptr<VectorDouble>(new VectorDouble());
109 eshelbyStressMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
110
111 eshelbyStress_dx = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
112 }
113
116 };
117
118 template <bool D = false>
120
121 OpCalculateStrain(const std::string row_field, const std::string col_field,
122 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
123
124 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
125
126 private:
127 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
128 };
129
131
132 OpCalculateStrainAle(const std::string row_field,
133 const std::string col_field,
134 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
135
136 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
137
138 private:
139 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
140 };
141
142#define MAT_TO_DDG(SM) \
143 &(*SM)(0, 0), &(*SM)(1, 0), &(*SM)(2, 0), &(*SM)(3, 0), &(*SM)(4, 0), \
144 &(*SM)(5, 0), &(*SM)(6, 0), &(*SM)(7, 0), &(*SM)(8, 0), &(*SM)(9, 0), \
145 &(*SM)(10, 0), &(*SM)(11, 0), &(*SM)(12, 0), &(*SM)(13, 0), \
146 &(*SM)(14, 0), &(*SM)(15, 0), &(*SM)(16, 0), &(*SM)(17, 0), \
147 &(*SM)(18, 0), &(*SM)(19, 0), &(*SM)(20, 0), &(*SM)(21, 0), \
148 &(*SM)(22, 0), &(*SM)(23, 0), &(*SM)(24, 0), &(*SM)(25, 0), \
149 &(*SM)(26, 0), &(*SM)(27, 0), &(*SM)(28, 0), &(*SM)(29, 0), \
150 &(*SM)(30, 0), &(*SM)(31, 0), &(*SM)(32, 0), &(*SM)(33, 0), \
151 &(*SM)(34, 0), &(*SM)(35, 0)
152
153 template <int S = 0> struct OpCalculateStress : public VolUserDataOperator {
154
155 OpCalculateStress(const std::string row_field, const std::string col_field,
156 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
157
158 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
159
160 protected:
161 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
162 };
163
165
166 OpCalculateEnergy(const std::string row_field, const std::string col_field,
167 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
169
170 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
171
172 protected:
173 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
175 };
176
178
180 const std::string row_field, const std::string col_field,
181 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
182
183 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
184
185 protected:
186 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
187 };
188
189 template <int S = 0>
191
193 const std::string row_field, const std::string col_field,
194 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
195 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
196
197 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
198
199 protected:
200 boost::shared_ptr<map<int, BlockData>>
201 blockSetsPtr; ///< Structure keeping data about problem, like
202 ///< material parameters
203 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
204 };
205
206 /** * @brief Assemble mass matrix for elastic element TODO: CHANGE FORMULA *
207 * \f[
208 * {\bf{M}} = \int\limits_\Omega
209 * \f]
210 *
211 */
214
219
220 boost::shared_ptr<DataAtIntegrationPts> commonData;
221
222 OpCalculateMassMatrix(const std::string row_field,
223 const std::string col_field, BlockData &data,
224 MassBlockData &mass_data,
225 boost::shared_ptr<DataAtIntegrationPts> &common_data,
226 bool symm = true)
228 row_field, col_field, OPROWCOL, symm),
229 commonData(common_data), dAta(data), massData(mass_data) {}
230
231 PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
232 EntityType col_type,
234 EntitiesFieldData::EntData &col_data) {
236
237 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
239 &m(3 * r + 0, 3 * c + 0), &m(3 * r + 0, 3 * c + 1),
240 &m(3 * r + 0, 3 * c + 2), &m(3 * r + 1, 3 * c + 0),
241 &m(3 * r + 1, 3 * c + 1), &m(3 * r + 1, 3 * c + 2),
242 &m(3 * r + 2, 3 * c + 0), &m(3 * r + 2, 3 * c + 1),
243 &m(3 * r + 2, 3 * c + 2));
244 };
245
246 const int row_nb_dofs = row_data.getIndices().size();
247 if (!row_nb_dofs)
249 const int col_nb_dofs = col_data.getIndices().size();
250 if (!col_nb_dofs)
252 if (dAta.tEts.find(getFEEntityHandle()) == dAta.tEts.end()) {
254 }
255 if (massData.tEts.find(getFEEntityHandle()) == massData.tEts.end()) {
257 }
258
259 const bool diagonal_block =
260 (row_type == col_type) && (row_side == col_side);
261 // get number of integration points
262 // Set size can clear local tangent matrix
263 locK.resize(row_nb_dofs, col_nb_dofs, false);
264 locK.clear();
265
266 const int row_nb_gauss_pts = row_data.getN().size1();
267 const int row_nb_base_functions = row_data.getN().size2();
268
273
274 double density = massData.rho0;
275
276 // get integration weights
277 auto t_w = getFTensor0IntegrationWeight();
278
279 // integrate local matrix for entity block
280 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
281
282 auto t_row_base_func = row_data.getFTensor0N(gg, 0);
283
284 // Get volume and integration weight
285 double w = getVolume() * t_w;
286
287 for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
288 auto t_col_base_func = col_data.getFTensor0N(gg, 0);
289 for (int col_bb = 0; col_bb != col_nb_dofs / 3; col_bb++) {
290 auto t_assemble = get_tensor2(locK, row_bb, col_bb);
291 t_assemble(i, j) += density * t_row_base_func * t_col_base_func * w;
292 // Next base function for column
293 ++t_col_base_func;
294 }
295 // Next base function for row
296 ++t_row_base_func;
297 }
298 // Next integration point for getting weight
299 ++t_w;
300 }
301
302 CHKERR MatSetValues(getKSPB(), row_data, col_data, &locK(0, 0),
303 ADD_VALUES);
304
305 // is symmetric
306 if (row_type != col_type || row_side != col_side) {
307 translocK.resize(col_nb_dofs, row_nb_dofs, false);
308 noalias(translocK) = trans(locK);
309
310 CHKERR MatSetValues(getKSPB(), col_data, row_data, &translocK(0, 0),
311 ADD_VALUES);
312 }
313
315 }
316 };
317
319 protected:
320 boost::shared_ptr<map<int, BlockData>>
321 blockSetsPtr; ///< Structure keeping data about problem, like
322 ///< material parameters
323 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
324
325 boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
326 const double rhoN; ///< exponent n in E(p) = E * (p / p_0)^n
327 const double rHo0; ///< p_0 reference density in E(p) = E * (p / p_0)^n
328 // // where p is density, E - youngs modulus
329 public:
331 const std::string row_field, const std::string col_field,
332 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
333 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
334 boost::shared_ptr<VectorDouble> rho_at_gauss_pts, const double rho_n,
335 const double rho_0);
336
337 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
338 };
339
341
342 OpAssemble(const std::string row_field, const std::string col_field,
343 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
344 const char type, bool symm = false);
345
346 /**
347 * \brief Do calculations for give operator
348 * @param row_side row side number (local number) of entity on element
349 * @param col_side column side number (local number) of entity on element
350 * @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
351 * @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
352 * @param row_data data for row
353 * @param col_data data for column
354 * @return error code
355 */
356 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
357 EntityType col_type, EntData &row_data,
358 EntData &col_data);
359
360 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
361
362 protected:
363 // Finite element stiffness sub-matrix K_ij
367
368 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
369
372
373 int nbRows; ///< number of dofs on rows
374 int nbCols; ///< number if dof on column
375 int nbIntegrationPts; ///< number of integration points
376 bool isDiag; ///< true if this block is on diagonal
377
378 virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
379
380 virtual MoFEMErrorCode iNtegrate(EntData &row_data);
381
382 /**
383 * \brief Assemble local entity block matrix
384 * @param row_data row data (consist base functions on row entity)
385 * @param col_data column data (consist base functions on column
386 * entity)
387 * @return error code
388 */
389 MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data);
390
391 /**
392 * \brief Assemble local entity right-hand vector
393 * @param row_data row data (consist base functions on row entity)
394 * @param col_data column data (consist base functions on column
395 * entity)
396 * @return error code
397 */
399 };
400
401 struct OpRhs_dx : public OpAssemble {
402
403 OpRhs_dx(const std::string row_field, const std::string col_field,
404 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
405
406 protected:
408 };
409
410 template <int S = 0> struct OpLhs_dx_dx : public OpAssemble {
411
412 OpLhs_dx_dx(const std::string row_field, const std::string col_field,
413 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
414
415 protected:
416 /**
417 * \brief Integrate B^T D B operator
418 * @param row_data row data (consist base functions on row entity)
419 * @param col_data column data (consist base functions on column entity)
420 * @return error code
421 */
423 };
424
425 struct OpAleRhs_dx : public OpAssemble {
426
427 OpAleRhs_dx(const std::string row_field, const std::string col_field,
428 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
429
430 protected:
432 };
433
434 template <int S = 0> struct OpAleLhs_dx_dx : public OpAssemble {
435
436 OpAleLhs_dx_dx(const std::string row_field, const std::string col_field,
437 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
438
439 protected:
440 /**
441 * \brief Integrate B^T D B operator
442 * @param row_data row data (consist base functions on row entity)
443 * @param col_data column data (consist base functions on column entity)
444 * @return error code
445 */
447 };
448
449 template <int S = 0> struct OpAleLhs_dx_dX : public OpAssemble {
450
451 OpAleLhs_dx_dX(const std::string row_field, const std::string col_field,
452 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
453
454 protected:
455 /**
456 * \brief Integrate tangent stiffness for spatial momentum
457 * @param row_data row data (consist base functions on row entity)
458 * @param col_data column data (consist base functions on column entity)
459 * @return error code
460 */
462 };
463
465
466 boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
467 boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
468 const double rhoN;
469 const double rHo0;
470
472 const std::string row_field, const std::string col_field,
473 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
474 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
475 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
476 const double rho_n, const double rho_0);
477
478 protected:
479 /**
480 * \brief Integrate tangent stiffness for spatial momentum
481 * @param row_data row data (consist base functions on row entity)
482 * @param col_data column data (consist base functions on column entity)
483 * @return error code
484 */
486 };
487
489
490 boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
491 boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
492 const double rhoN;
493 const double rHo0;
494
496 const std::string row_field, const std::string col_field,
497 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
498 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
499 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
500 const double rho_n, const double rho_0);
501
502 protected:
503 /**
504 * \brief Integrate tangent stiffness for material momentum
505 * @param row_data row data (consist base functions on row entity)
506 * @param col_data column data (consist base functions on column entity)
507 * @return error code
508 */
510 };
511
512 struct OpAleRhs_dX : public OpAssemble {
513
514 OpAleRhs_dX(const std::string row_field, const std::string col_field,
515 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
516
517 protected:
519 };
520
521 template <int S = 0> struct OpAleLhs_dX_dX : public OpAssemble {
522
523 OpAleLhs_dX_dX(const std::string row_field, const std::string col_field,
524 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
525
526 protected:
527 /**
528 * \brief Integrate tangent stiffness for material momentum
529 * @param row_data row data (consist base functions on row entity)
530 * @param col_data column data (consist base functions on column entity)
531 * @return error code
532 */
534 };
535
536 template <int S = 0> struct OpAleLhsPre_dX_dx : public VolUserDataOperator {
537
538 OpAleLhsPre_dX_dx(const std::string row_field, const std::string col_field,
539 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
540
541 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
542
543 private:
544 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
545 };
546
547 struct OpAleLhs_dX_dx : public OpAssemble {
548
549 OpAleLhs_dX_dx(const std::string row_field, const std::string col_field,
550 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
551 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
552
553 protected:
554 /**
555 * \brief Integrate tangent stiffness for material momentum
556 * @param row_data row data (consist base functions on row entity)
557 * @param col_data column data (consist base functions on column entity)
558 * @return error code
559 */
561 };
562
563 template <int S> struct OpAnalyticalInternalStrain_dx : public OpAssemble {
564
565 typedef boost::function<
566
568
570
571 )
572
573 >
575
577 const std::string row_field,
578 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
579 StrainFunction strain_fun);
580
581 protected:
584 };
585
586 template <int S> struct OpAnalyticalInternalAleStrain_dX : public OpAssemble {
587
588 typedef boost::function<
589
591
593
594 )
595
596 >
598
600 const std::string row_field,
601 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
602 StrainFunction strain_fun,
603 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
604
605 protected:
608 boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
609 };
610
611 template <int S> struct OpAnalyticalInternalAleStrain_dx : public OpAssemble {
612
613 typedef boost::function<
614
616
618
619 )
620
621 >
623
625 const std::string row_field,
626 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
627 StrainFunction strain_fun,
628 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
629
630 protected:
633 boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
634 };
635
636 template <class ELEMENT>
637 struct OpPostProcHookeElement : public ELEMENT::UserDataOperator {
638 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
639 map<int, BlockData>
640 &blockSetsPtr; // FIXME: (works only with the first block)
641 moab::Interface &postProcMesh;
642 std::vector<EntityHandle> &mapGaussPts;
643 bool isALE;
645
646 OpPostProcHookeElement(const string row_field,
647 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
648 map<int, BlockData> &block_sets_ptr,
649 moab::Interface &post_proc_mesh,
650 std::vector<EntityHandle> &map_gauss_pts,
651 bool is_ale = false, bool is_field_disp = true);
652
655 };
656
657 static MoFEMErrorCode
659 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr);
660
661 static MoFEMErrorCode
663 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
664 const std::string element_name, const std::string x_field,
665 const std::string X_field, const bool ale);
666
667 static MoFEMErrorCode
668 setOperators(boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr,
669 boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr,
670 boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
671 const std::string x_field, const std::string X_field,
672 const bool ale, const bool field_disp,
673 const EntityType type = MBTET,
674 boost::shared_ptr<DataAtIntegrationPts> data_at_pts = nullptr);
675
676 static MoFEMErrorCode
677 calculateEnergy(DM dm, boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
678 const std::string x_field, const std::string X_field,
679 const bool ale, const bool field_disp,
680 SmartPetscObj<Vec> &v_energy_ptr);
681
682private:
684};
685
686template <bool D>
687HookeElement::OpCalculateStrain<D>::OpCalculateStrain(
688 const std::string row_field, const std::string col_field,
689 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
690 : VolUserDataOperator(row_field, col_field, OPROW, true),
691 dataAtPts(data_at_pts) {
692 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
693}
694
695template <bool D>
696MoFEMErrorCode HookeElement::OpCalculateStrain<D>::doWork(int row_side,
697 EntityType row_type,
698 EntData &row_data) {
702 // get number of integration points
703 const int nb_integration_pts = getGaussPts().size2();
704 dataAtPts->smallStrainMat->resize(6, nb_integration_pts, false);
705 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
706 auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hMat));
707
708 for (int gg = 0; gg != nb_integration_pts; ++gg) {
709 t_strain(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
710
711 // If displacement field, not field o spatial positons is given
712 if (!D) {
713 t_strain(0, 0) -= 1;
714 t_strain(1, 1) -= 1;
715 t_strain(2, 2) -= 1;
716 }
717
718 ++t_strain;
719 ++t_h;
720 }
722}
723
724template <int S>
725HookeElement::OpAleLhs_dx_dx<S>::OpAleLhs_dx_dx(
726 const std::string row_field, const std::string col_field,
727 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
728 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
729
730template <int S>
731MoFEMErrorCode HookeElement::OpAleLhs_dx_dx<S>::iNtegrate(EntData &row_data,
732 EntData &col_data) {
734
735 // get sub-block (3x3) of local stiffens matrix, here represented by
736 // second order tensor
737 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
739 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
740 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
741 &m(r + 2, c + 2));
742 };
743
748
749 // get element volume
750 double vol = getVolume();
751
752 // get intergrayion weights
753 auto t_w = getFTensor0IntegrationWeight();
754
755 // get derivatives of base functions on rows
756 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
757 const int row_nb_base_fun = row_data.getN().size2();
758
759 // Elastic stiffness tensor (4th rank tensor with minor and major
760 // symmetry)
762 MAT_TO_DDG(dataAtPts->stiffnessMat));
763
764 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
765 auto &det_H = *dataAtPts->detHVec;
766
767 // iterate over integration points
768 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
769
770 // calculate scalar weight times element volume
771 double a = t_w * vol * det_H[gg];
772
773 // iterate over row base functions
774 int rr = 0;
775 for (; rr != nbRows / 3; ++rr) {
776
777 // get sub matrix for the row
778 auto t_m = get_tensor2(K, 3 * rr, 0);
779
780 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
781 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
782
784 // I mix up the indices here so that it behaves like a
785 // Dg. That way I don't have to have a separate wrapper
786 // class Christof_Expr, which simplifies things.
787 t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base_pulled(i));
788
789 // get derivatives of base functions for columns
790 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
791
792 // iterate column base functions
793 for (int cc = 0; cc != nbCols / 3; ++cc) {
794
795 FTensor::Tensor1<double, 3> t_col_diff_base_pulled;
796 t_col_diff_base_pulled(j) = t_col_diff_base(i) * t_invH(i, j);
797
798 // integrate block local stiffens matrix
799 t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base_pulled(k);
800
801 // move to next column base function
802 ++t_col_diff_base;
803
804 // move to next block of local stiffens matrix
805 ++t_m;
806 }
807
808 // move to next row base function
809 ++t_row_diff_base;
810 }
811
812 for (; rr != row_nb_base_fun; ++rr)
813 ++t_row_diff_base;
814
815 // move to next integration weight
816 ++t_w;
817 ++t_D;
818 ++t_invH;
819 }
820
822}
823
824template <int S>
825HookeElement::OpCalculateHomogeneousStiffness<S>::
826 OpCalculateHomogeneousStiffness(
827 const std::string row_field, const std::string col_field,
828 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
829 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
830 : VolUserDataOperator(row_field, col_field, OPROW, true),
831 blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts) {
832 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
833}
834
835template <int S>
836MoFEMErrorCode HookeElement::OpCalculateHomogeneousStiffness<S>::doWork(
837 int row_side, EntityType row_type, EntData &row_data) {
839
840 for (auto &m : (*blockSetsPtr)) {
841 if (m.second.tEts.find(getFEEntityHandle()) != m.second.tEts.end()) {
842
843 dataAtPts->stiffnessMat->resize(36, 1, false);
845 MAT_TO_DDG(dataAtPts->stiffnessMat));
846 const double young = m.second.E;
847 const double poisson = m.second.PoissonRatio;
848
849 // coefficient used in intermediate calculation
850 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
851
856
857 t_D(i, j, k, l) = 0.;
858
859 t_D(0, 0, 0, 0) = 1 - poisson;
860 t_D(1, 1, 1, 1) = 1 - poisson;
861 t_D(2, 2, 2, 2) = 1 - poisson;
862
863 t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * poisson);
864 t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * poisson);
865 t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * poisson);
866
867 t_D(0, 0, 1, 1) = poisson;
868 t_D(1, 1, 0, 0) = poisson;
869 t_D(0, 0, 2, 2) = poisson;
870 t_D(2, 2, 0, 0) = poisson;
871 t_D(1, 1, 2, 2) = poisson;
872 t_D(2, 2, 1, 1) = poisson;
873
874 t_D(i, j, k, l) *= coefficient;
875
876 break;
877 }
878 }
879
881}
882
883template <int S>
884HookeElement::OpCalculateStress<S>::OpCalculateStress(
885 const std::string row_field, const std::string col_field,
886 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
887 : VolUserDataOperator(row_field, col_field, OPROW, true),
888 dataAtPts(data_at_pts) {
889 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
890}
891
892template <int S>
893MoFEMErrorCode HookeElement::OpCalculateStress<S>::doWork(int row_side,
894 EntityType row_type,
895 EntData &row_data) {
897 // get number of integration points
898 const int nb_integration_pts = getGaussPts().size2();
899 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
900 dataAtPts->cauchyStressMat->resize(6, nb_integration_pts, false);
901 auto t_cauchy_stress =
902 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
903
908
909 // elastic stiffness tensor (4th rank tensor with minor and major
910 // symmetry)
912 MAT_TO_DDG(dataAtPts->stiffnessMat));
913 for (int gg = 0; gg != nb_integration_pts; ++gg) {
914 t_cauchy_stress(i, j) = t_D(i, j, k, l) * t_strain(k, l);
915 ++t_strain;
916 ++t_cauchy_stress;
917 ++t_D;
918 }
920}
921
922template <int S>
923HookeElement::OpLhs_dx_dx<S>::OpLhs_dx_dx(
924 const std::string row_field, const std::string col_field,
925 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
926 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
927
928template <int S>
929MoFEMErrorCode HookeElement::OpLhs_dx_dx<S>::iNtegrate(EntData &row_data,
930 EntData &col_data) {
932
933 // get sub-block (3x3) of local stiffens matrix, here represented by
934 // second order tensor
935 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
937 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
938 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
939 &m(r + 2, c + 2));
940 };
941
946
947 // get element volume
948 double vol = getVolume();
949
950 // get intergrayion weights
951 auto t_w = getFTensor0IntegrationWeight();
952
953 // get derivatives of base functions on rows
954 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
955 const int row_nb_base_fun = row_data.getN().size2();
956
957 // Elastic stiffness tensor (4th rank tensor with minor and major
958 // symmetry)
960 MAT_TO_DDG(dataAtPts->stiffnessMat));
961
962 // iterate over integration points
963 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
964
965 // calculate scalar weight times element volume
966 double a = t_w * vol;
967
968 // iterate over row base functions
969 int rr = 0;
970 for (; rr != nbRows / 3; ++rr) {
971
972 // get sub matrix for the row
973 auto t_m = get_tensor2(K, 3 * rr, 0);
974
975 // get derivatives of base functions for columns
976 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
977
979 // I mix up the indices here so that it behaves like a
980 // Dg. That way I don't have to have a separate wrapper
981 // class Christof_Expr, which simplifies things.
982 t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
983
984 // iterate column base functions
985 for (int cc = 0; cc != nbCols / 3; ++cc) {
986
987 // integrate block local stiffens matrix
988 t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
989
990 // move to next column base function
991 ++t_col_diff_base;
992
993 // move to next block of local stiffens matrix
994 ++t_m;
995 }
996
997 // move to next row base function
998 ++t_row_diff_base;
999 }
1000
1001 for (; rr != row_nb_base_fun; ++rr)
1002 ++t_row_diff_base;
1003
1004 // move to next integration weight
1005 ++t_w;
1006 ++t_D;
1007 }
1008
1010}
1011
1012template <int S>
1013HookeElement::OpAleLhs_dx_dX<S>::OpAleLhs_dx_dX(
1014 const std::string row_field, const std::string col_field,
1015 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1016 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
1017
1018template <int S>
1019MoFEMErrorCode HookeElement::OpAleLhs_dx_dX<S>::iNtegrate(EntData &row_data,
1020 EntData &col_data) {
1022
1023 // get sub-block (3x3) of local stiffens matrix, here represented by
1024 // second order tensor
1025 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1027 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1028 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1029 &m(r + 2, c + 2));
1030 };
1031
1038
1039 // get element volume
1040 double vol = getVolume();
1041
1042 // get intergrayion weights
1043 auto t_w = getFTensor0IntegrationWeight();
1044
1045 // get derivatives of base functions on rows
1046 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1047 const int row_nb_base_fun = row_data.getN().size2();
1048
1049 // Elastic stiffness tensor (4th rank tensor with minor and major
1050 // symmetry)
1052 MAT_TO_DDG(dataAtPts->stiffnessMat));
1053
1054 auto t_cauchy_stress =
1055 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1056 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1057 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1058 auto &det_H = *dataAtPts->detHVec;
1059
1060 // iterate over integration points
1061 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1062
1063 // calculate scalar weight times element volume
1064 double a = t_w * vol * det_H[gg];
1065
1067 t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1068
1069 // iterate over row base functions
1070 int rr = 0;
1071 for (; rr != nbRows / 3; ++rr) {
1072
1073 // get sub matrix for the row
1074 auto t_m = get_tensor2(K, 3 * rr, 0);
1075
1076 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1077 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1078
1079 FTensor::Tensor1<double, 3> t_row_stress;
1080 t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_cauchy_stress(i, j);
1081
1082 FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
1083 t_row_diff_base_pulled_dX(j, k, l) =
1084 -(t_invH(i, k) * t_row_diff_base(i)) * t_invH(l, j);
1085
1086 FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
1087 t_row_dX_stress(i, k, l) =
1088 a * (t_row_diff_base_pulled_dX(j, k, l) * t_cauchy_stress(j, i));
1089
1091 t_row_D(l, j, k) = (a * t_row_diff_base_pulled(i)) * t_D(i, j, k, l);
1092
1093 FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1094 // FIXME: This operator is not implemented, doing operation by hand
1095 // t_row_stress_dX(i, m, n) = t_row_D(i, k, l) * t_F_dX(k, l, m, n);
1096 t_row_stress_dX(i, j, k) = 0;
1097 for (int ii = 0; ii != 3; ++ii)
1098 for (int mm = 0; mm != 3; ++mm)
1099 for (int nn = 0; nn != 3; ++nn) {
1100 auto &v = t_row_stress_dX(ii, mm, nn);
1101 for (int kk = 0; kk != 3; ++kk)
1102 for (int ll = 0; ll != 3; ++ll)
1103 v += t_row_D(ii, kk, ll) * t_F_dX(kk, ll, mm, nn);
1104 }
1105
1106 // get derivatives of base functions for columns
1107 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1108
1109 // iterate column base functions
1110 for (int cc = 0; cc != nbCols / 3; ++cc) {
1111
1112 t_m(i, k) += t_row_stress(i) * (t_invH(j, k) * t_col_diff_base(j));
1113 t_m(i, k) += t_row_dX_stress(i, k, l) * t_col_diff_base(l);
1114 t_m(i, k) += t_row_stress_dX(i, k, l) * t_col_diff_base(l);
1115
1116 // move to next column base function
1117 ++t_col_diff_base;
1118
1119 // move to next block of local stiffens matrix
1120 ++t_m;
1121 }
1122
1123 // move to next row base function
1124 ++t_row_diff_base;
1125 }
1126
1127 for (; rr != row_nb_base_fun; ++rr)
1128 ++t_row_diff_base;
1129
1130 // move to next integration weight
1131 ++t_w;
1132 ++t_D;
1133 ++t_cauchy_stress;
1134 ++t_invH;
1135 ++t_h;
1136 }
1137
1139}
1140
1141template <int S>
1142HookeElement::OpAleLhs_dX_dX<S>::OpAleLhs_dX_dX(
1143 const std::string row_field, const std::string col_field,
1144 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1145 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
1146
1147template <int S>
1148MoFEMErrorCode HookeElement::OpAleLhs_dX_dX<S>::iNtegrate(EntData &row_data,
1149 EntData &col_data) {
1151
1152 // get sub-block (3x3) of local stiffens matrix, here represented by
1153 // second order tensor
1154 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1156 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1157 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1158 &m(r + 2, c + 2));
1159 };
1160
1167
1168 // get element volume
1169 double vol = getVolume();
1170
1171 // get intergrayion weights
1172 auto t_w = getFTensor0IntegrationWeight();
1173
1174 // get derivatives of base functions on rows
1175 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1176 const int row_nb_base_fun = row_data.getN().size2();
1177
1178 // Elastic stiffness tensor (4th rank tensor with minor and major
1179 // symmetry)
1181 MAT_TO_DDG(dataAtPts->stiffnessMat));
1182 auto t_cauchy_stress =
1183 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1184 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
1185 auto t_eshelby_stress =
1186 getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
1187 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1188 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1189 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1190 auto &det_H = *dataAtPts->detHVec;
1191
1192 // iterate over integration points
1193 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1194
1195 // calculate scalar weight times element volume
1196 double a = t_w * vol * det_H[gg];
1197
1199 t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1200
1202 t_D_strain_dX(i, j, m, n) = 0.;
1203 for (int ii = 0; ii != 3; ++ii)
1204 for (int jj = 0; jj != 3; ++jj)
1205 for (int ll = 0; ll != 3; ++ll)
1206 for (int kk = 0; kk != 3; ++kk) {
1207 auto &v = t_D_strain_dX(ii, jj, kk, ll);
1208 for (int mm = 0; mm != 3; ++mm)
1209 for (int nn = 0; nn != 3; ++nn)
1210 v += t_D(ii, jj, mm, nn) * t_F_dX(mm, nn, kk, ll);
1211 }
1212
1213 FTensor::Tensor4<double, 3, 3, 3, 3> t_eshelby_stress_dX;
1214 t_eshelby_stress_dX(i, j, m, n) = t_F(k, i) * t_D_strain_dX(k, j, m, n);
1215
1216 for (int ii = 0; ii != 3; ++ii)
1217 for (int jj = 0; jj != 3; ++jj)
1218 for (int mm = 0; mm != 3; ++mm)
1219 for (int nn = 0; nn != 3; ++nn) {
1220 auto &v = t_eshelby_stress_dX(ii, jj, mm, nn);
1221 for (int kk = 0; kk != 3; ++kk)
1222 v += t_F_dX(kk, ii, mm, nn) * t_cauchy_stress(kk, jj);
1223 }
1224
1225 t_eshelby_stress_dX(i, j, k, l) *= -1;
1226
1228 t_energy_dX(k, l) = t_F_dX(i, j, k, l) * t_cauchy_stress(i, j);
1229 t_energy_dX(k, l) +=
1230 (t_strain(m, n) * t_D(m, n, i, j)) * t_F_dX(i, j, k, l);
1231 t_energy_dX(k, l) /= 2.;
1232
1233 for (int kk = 0; kk != 3; ++kk)
1234 for (int ll = 0; ll != 3; ++ll) {
1235 auto v = t_energy_dX(kk, ll);
1236 for (int ii = 0; ii != 3; ++ii)
1237 t_eshelby_stress_dX(ii, ii, kk, ll) += v;
1238 }
1239
1240 // iterate over row base functions
1241 int rr = 0;
1242 for (; rr != nbRows / 3; ++rr) {
1243
1244 // get sub matrix for the row
1245 auto t_m = get_tensor2(K, 3 * rr, 0);
1246
1247 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1248 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1249
1250 FTensor::Tensor1<double, 3> t_row_stress;
1251 t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1252
1253 FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
1254 t_row_diff_base_pulled_dX(j, k, l) =
1255 -(t_row_diff_base(i) * t_invH(i, k)) * t_invH(l, j);
1256
1257 FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
1258 t_row_dX_stress(i, k, l) =
1259 a * (t_row_diff_base_pulled_dX(j, k, l) * t_eshelby_stress(i, j));
1260
1261 FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1262 t_row_stress_dX(i, m, n) =
1263 a * t_row_diff_base_pulled(j) * t_eshelby_stress_dX(i, j, m, n);
1264
1265 // get derivatives of base functions for columns
1266 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1267
1268 // iterate column base functions
1269 for (int cc = 0; cc != nbCols / 3; ++cc) {
1270
1271 t_m(i, k) += t_row_stress(i) * (t_invH(j, k) * t_col_diff_base(j));
1272 t_m(i, k) += t_row_dX_stress(i, k, l) * t_col_diff_base(l);
1273 t_m(i, k) += t_row_stress_dX(i, k, l) * t_col_diff_base(l);
1274
1275 // move to next column base function
1276 ++t_col_diff_base;
1277
1278 // move to next block of local stiffens matrix
1279 ++t_m;
1280 }
1281
1282 // move to next row base function
1283 ++t_row_diff_base;
1284 }
1285
1286 for (; rr != row_nb_base_fun; ++rr)
1287 ++t_row_diff_base;
1288
1289 // move to next integration weight
1290 ++t_w;
1291 ++t_D;
1292 ++t_cauchy_stress;
1293 ++t_strain;
1294 ++t_eshelby_stress;
1295 ++t_h;
1296 ++t_invH;
1297 ++t_F;
1298 }
1299
1301}
1302
1303template <int S>
1304HookeElement::OpAleLhsPre_dX_dx<S>::OpAleLhsPre_dX_dx(
1305 const std::string row_field, const std::string col_field,
1306 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1307 : VolUserDataOperator(row_field, col_field, OPROW, true),
1308 dataAtPts(data_at_pts) {
1309 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1310}
1311
1312template <int S>
1313MoFEMErrorCode HookeElement::OpAleLhsPre_dX_dx<S>::doWork(int row_side,
1314 EntityType row_type,
1315 EntData &row_data) {
1317
1318 const int nb_integration_pts = row_data.getN().size1();
1319
1320 auto get_eshelby_stress_dx = [this, nb_integration_pts]() {
1322 t_eshelby_stress_dx;
1323 dataAtPts->eshelbyStress_dx->resize(81, nb_integration_pts, false);
1324 int mm = 0;
1325 for (int ii = 0; ii != 3; ++ii)
1326 for (int jj = 0; jj != 3; ++jj)
1327 for (int kk = 0; kk != 3; ++kk)
1328 for (int ll = 0; ll != 3; ++ll)
1329 t_eshelby_stress_dx.ptr(ii, jj, kk, ll) =
1330 &(*dataAtPts->eshelbyStress_dx)(mm++, 0);
1331 return t_eshelby_stress_dx;
1332 };
1333
1334 auto t_eshelby_stress_dx = get_eshelby_stress_dx();
1335
1342
1343 // Elastic stiffness tensor (4th rank tensor with minor and major
1344 // symmetry)
1346 MAT_TO_DDG(dataAtPts->stiffnessMat));
1347 auto t_cauchy_stress =
1348 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1349 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1350 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1351
1352 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1353
1354 t_eshelby_stress_dx(i, j, m, n) =
1355 (t_F(k, i) * t_D(k, j, m, l)) * t_invH(n, l);
1356 for (int ii = 0; ii != 3; ++ii)
1357 for (int jj = 0; jj != 3; ++jj)
1358 for (int mm = 0; mm != 3; ++mm)
1359 for (int nn = 0; nn != 3; ++nn) {
1360 auto &v = t_eshelby_stress_dx(ii, jj, mm, nn);
1361 v += t_invH(nn, ii) * t_cauchy_stress(mm, jj);
1362 }
1363 t_eshelby_stress_dx(i, j, k, l) *= -1;
1364
1366 t_energy_dx(m, n) = t_invH(n, j) * t_cauchy_stress(m, j);
1367
1368 for (int mm = 0; mm != 3; ++mm)
1369 for (int nn = 0; nn != 3; ++nn) {
1370 auto v = t_energy_dx(mm, nn);
1371 for (int ii = 0; ii != 3; ++ii)
1372 t_eshelby_stress_dx(ii, ii, mm, nn) += v;
1373 }
1374
1375 ++t_D;
1376 ++t_invH;
1377 ++t_cauchy_stress;
1378 ++t_eshelby_stress_dx;
1379 ++t_F;
1380 }
1381
1383}
1384
1385template <class ELEMENT>
1386HookeElement::OpPostProcHookeElement<ELEMENT>::OpPostProcHookeElement(
1387 const string row_field, boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1388 map<int, BlockData> &block_sets_ptr, moab::Interface &post_proc_mesh,
1389 std::vector<EntityHandle> &map_gauss_pts, bool is_ale, bool is_field_disp)
1390 : ELEMENT::UserDataOperator(row_field, UserDataOperator::OPROW),
1391 dataAtPts(data_at_pts), blockSetsPtr(block_sets_ptr),
1392 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), isALE(is_ale),
1393 isFieldDisp(is_field_disp) {}
1394
1395template <class ELEMENT>
1396MoFEMErrorCode HookeElement::OpPostProcHookeElement<ELEMENT>::doWork(
1397 int side, EntityType type, EntitiesFieldData::EntData &data) {
1399
1400 if (type != MBVERTEX) {
1402 }
1403
1404 auto tensor_to_tensor = [](const auto &t1, auto &t2) {
1405 t2(0, 0) = t1(0, 0);
1406 t2(1, 1) = t1(1, 1);
1407 t2(2, 2) = t1(2, 2);
1408 t2(0, 1) = t2(1, 0) = t1(1, 0);
1409 t2(0, 2) = t2(2, 0) = t1(2, 0);
1410 t2(1, 2) = t2(2, 1) = t1(2, 1);
1411 };
1412
1413 std::array<double, 9> def_val;
1414 def_val.fill(0);
1415
1416 auto make_tag = [&](auto name, auto size) {
1417 Tag th;
1418 CHKERR postProcMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
1419 MB_TAG_CREAT | MB_TAG_SPARSE,
1420 def_val.data());
1421 return th;
1422 };
1423
1424 auto th_stress = make_tag("STRESS", 9);
1425 auto th_psi = make_tag("ENERGY", 1);
1426
1427 const int nb_integration_pts = mapGaussPts.size();
1428
1433
1434 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1435 auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
1436
1437 dataAtPts->stiffnessMat->resize(36, 1, false);
1439 MAT_TO_DDG(dataAtPts->stiffnessMat));
1440
1441 EntityHandle ent = this->getFEEntityHandle();
1442 auto type = type_from_handle(ent);
1443 EntityHandle ent_3d = ent;
1444 if (type == MBTRI || type == MBQUAD) {
1445 Range ents;
1446 auto &m_field = this->getPtrFE()->mField;
1447 CHKERR m_field.get_moab().get_adjacencies(&ent, 1, 3, false, ents,
1448 moab::Interface::UNION);
1449#ifndef NDEBUG
1450 if (ents.empty())
1451 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1452 "Could not find a 3D element adjacent to a given face element");
1453#endif
1454 ent_3d = ents.front();
1455 }
1456
1457 bool found_block = false;
1458 for (auto &m : (blockSetsPtr)) {
1459 if (m.second.tEts.find(ent_3d) != m.second.tEts.end()) {
1460 const double young = m.second.E;
1461 const double poisson = m.second.PoissonRatio;
1462 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
1463
1464 t_D(i, j, k, l) = 0.;
1465 t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
1466 t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
1467 0.5 * (1 - 2 * poisson);
1468 t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
1469 t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
1470 t_D(i, j, k, l) *= coefficient;
1471
1472 found_block = true;
1473 break;
1474 }
1475 }
1476 if (!found_block)
1477 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1478 "Element not found in any of material blocksets");
1479
1480 double detH = 0.;
1484 FTensor::Tensor2<double, 3, 3> t_small_strain;
1486 FTensor::Tensor2_symmetric<double, 3> t_small_strain_symm;
1487
1488 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1489
1490 if (isFieldDisp) {
1491 t_h(0, 0) += 1;
1492 t_h(1, 1) += 1;
1493 t_h(2, 2) += 1;
1494 }
1495
1496 if (!isALE) {
1497 t_small_strain_symm(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
1498 } else {
1499 CHKERR determinantTensor3by3(t_H, detH);
1500 CHKERR invertTensor3by3(t_H, detH, t_invH);
1501 t_F(i, j) = t_h(i, k) * t_invH(k, j);
1502 t_small_strain_symm(i, j) = (t_F(i, j) || t_F(j, i)) / 2.;
1503 ++t_H;
1504 }
1505
1506 t_small_strain_symm(0, 0) -= 1;
1507 t_small_strain_symm(1, 1) -= 1;
1508 t_small_strain_symm(2, 2) -= 1;
1509
1510 // symmetric tensors need improvement
1511 t_stress_symm(i, j) = t_D(i, j, k, l) * t_small_strain_symm(k, l);
1512 tensor_to_tensor(t_stress_symm, t_stress);
1513
1514 const double psi = 0.5 * t_stress_symm(i, j) * t_small_strain_symm(i, j);
1515
1516 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
1517 CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
1518 &t_stress(0, 0));
1519
1520 ++t_h;
1521 }
1522
1524}
1525
1526template <int S>
1527HookeElement::OpAnalyticalInternalStrain_dx<S>::OpAnalyticalInternalStrain_dx(
1528 const std::string row_field,
1529 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1530 StrainFunction strain_fun)
1531 : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1532 strainFun(strain_fun) {}
1533
1534template <int S>
1536HookeElement::OpAnalyticalInternalStrain_dx<S>::iNtegrate(EntData &row_data) {
1542
1543 auto get_tensor1 = [](VectorDouble &v, const int r) {
1545 &v(r + 0), &v(r + 1), &v(r + 2));
1546 };
1547
1548 const int nb_integration_pts = getGaussPts().size2();
1549 auto t_coords = getFTensor1CoordsAtGaussPts();
1550
1551 // get element volume
1552 double vol = getVolume();
1553 auto t_w = getFTensor0IntegrationWeight();
1554
1555 nF.resize(nbRows, false);
1556 nF.clear();
1557
1558 // elastic stiffness tensor (4th rank tensor with minor and major
1559 // symmetry)
1561 MAT_TO_DDG(dataAtPts->stiffnessMat));
1562
1563 // get derivatives of base functions on rows
1564 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1565 const int row_nb_base_fun = row_data.getN().size2();
1566
1567 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1568
1569 auto t_fun_strain = strainFun(t_coords);
1571 t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1572
1573 // calculate scalar weight times element volume
1574 double a = t_w * vol;
1575
1576 auto t_nf = get_tensor1(nF, 0);
1577
1578 int rr = 0;
1579 for (; rr != nbRows / 3; ++rr) {
1580 t_nf(i) += a * t_row_diff_base(j) * t_stress(i, j);
1581 ++t_row_diff_base;
1582 ++t_nf;
1583 }
1584
1585 for (; rr != row_nb_base_fun; ++rr)
1586 ++t_row_diff_base;
1587
1588 ++t_w;
1589 ++t_coords;
1590 ++t_D;
1591 }
1592
1594}
1595
1596template <int S>
1597HookeElement::OpAnalyticalInternalAleStrain_dX<S>::
1598 OpAnalyticalInternalAleStrain_dX(
1599 const std::string row_field,
1600 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1601 StrainFunction strain_fun,
1602 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1603 : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1604 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1605
1606template <int S>
1607MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dX<S>::iNtegrate(
1608 EntData &row_data) {
1614
1615 auto get_tensor1 = [](VectorDouble &v, const int r) {
1617 &v(r + 0), &v(r + 1), &v(r + 2));
1618 };
1619
1620 const int nb_integration_pts = getGaussPts().size2();
1621
1622 auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1623 auto t_coords = get_coords();
1624
1625 // get element volume
1626 double vol = getVolume();
1627 auto t_w = getFTensor0IntegrationWeight();
1628
1629 nF.resize(nbRows, false);
1630 nF.clear();
1631
1632 // elastic stiffness tensor (4th rank tensor with minor and major
1633 // symmetry)
1635 MAT_TO_DDG(dataAtPts->stiffnessMat));
1636 auto t_F = getFTensor2FromMat<3, 3>(*(dataAtPts->FMat));
1637 auto &det_H = *dataAtPts->detHVec;
1638 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1639
1640 // get derivatives of base functions on rows
1641 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1642 const int row_nb_base_fun = row_data.getN().size2();
1643
1644 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1645
1646 auto t_fun_strain = strainFun(t_coords);
1648 t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1649 FTensor::Tensor2<double, 3, 3> t_eshelby_stress;
1650 t_eshelby_stress(i, j) = -t_F(k, i) * t_stress(k, j);
1651
1652 // calculate scalar weight times element volume
1653 double a = t_w * vol * det_H[gg];
1654
1655 auto t_nf = get_tensor1(nF, 0);
1656
1657 int rr = 0;
1658 for (; rr != nbRows / 3; ++rr) {
1659 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1660 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1661 t_nf(i) += a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1662 ++t_row_diff_base;
1663 ++t_nf;
1664 }
1665
1666 for (; rr != row_nb_base_fun; ++rr)
1667 ++t_row_diff_base;
1668
1669 ++t_w;
1670 ++t_coords;
1671 ++t_F;
1672 ++t_invH;
1673 ++t_D;
1674 }
1675
1677}
1678
1679template <int S>
1680HookeElement::OpAnalyticalInternalAleStrain_dx<S>::
1681 OpAnalyticalInternalAleStrain_dx(
1682 const std::string row_field,
1683 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1684 StrainFunction strain_fun,
1685 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1686 : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1687 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1688
1689template <int S>
1690MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dx<S>::iNtegrate(
1691 EntData &row_data) {
1697
1698 auto get_tensor1 = [](VectorDouble &v, const int r) {
1700 &v(r + 0), &v(r + 1), &v(r + 2));
1701 };
1702
1703 const int nb_integration_pts = getGaussPts().size2();
1704
1705 auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1706 auto t_coords = get_coords();
1707
1708 // get element volume
1709 double vol = getVolume();
1710 auto t_w = getFTensor0IntegrationWeight();
1711
1712 nF.resize(nbRows, false);
1713 nF.clear();
1714
1715 // elastic stiffness tensor (4th rank tensor with minor and major
1716 // symmetry)
1718 MAT_TO_DDG(dataAtPts->stiffnessMat));
1719 auto &det_H = *dataAtPts->detHVec;
1720 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1721
1722 // get derivatives of base functions on rows
1723 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1724 const int row_nb_base_fun = row_data.getN().size2();
1725
1726 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1727
1728 auto t_fun_strain = strainFun(t_coords);
1730 t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1731
1732 // calculate scalar weight times element volume
1733 double a = t_w * vol * det_H[gg];
1734
1735 auto t_nf = get_tensor1(nF, 0);
1736
1737 int rr = 0;
1738 for (; rr != nbRows / 3; ++rr) {
1739 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1740 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1741 t_nf(i) += a * t_row_diff_base_pulled(j) * t_stress(i, j);
1742 ++t_row_diff_base;
1743 ++t_nf;
1744 }
1745
1746 for (; rr != row_nb_base_fun; ++rr)
1747 ++t_row_diff_base;
1748
1749 ++t_w;
1750 ++t_coords;
1751 ++t_invH;
1752 ++t_D;
1753 }
1754
1756}
1757
1758#endif // __HOOKE_ELEMENT_HPP
static MoFEMErrorCode setOperators(boost::shared_ptr< ForcesAndSourcesCore > fe_lhs_ptr, boost::shared_ptr< ForcesAndSourcesCore > fe_rhs_ptr, boost::shared_ptr< map< int, BlockData > > block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, const EntityType type=MBTET, boost::shared_ptr< DataAtIntegrationPts > data_at_pts=nullptr)
#define MAT_TO_DDG(SM)
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr)
MatrixDouble invJac
static MoFEMErrorCode addElasticElement(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr, const std::string element_name, const std::string x_field, const std::string X_field, const bool ale)
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData > > block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, SmartPetscObj< Vec > &v_energy_ptr)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
double D
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
UBlasVector< int > VectorInt
Definition: Types.hpp:67
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1761
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1486
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1475
int r
Definition: sdf.py:8
data for calculation inertia forces
VectorDouble a0
constant acceleration
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
boost::shared_ptr< MatrixDouble > HMat
boost::shared_ptr< MatrixDouble > invHMat
boost::shared_ptr< MatrixDouble > smallStrainMat
boost::shared_ptr< MatrixDouble > cauchyStressMat
boost::shared_ptr< MatrixDouble > eshelbyStress_dx
boost::shared_ptr< MatrixDouble > FMat
boost::shared_ptr< MatrixDouble > eshelbyStressMat
boost::shared_ptr< MatrixDouble > stiffnessMat
boost::shared_ptr< MatrixDouble > hMat
boost::shared_ptr< VectorDouble > energyVec
boost::shared_ptr< VectorDouble > detHVec
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
intrusive_ptr for managing petsc objects
data for calculation heat conductivity and heat capacity elements
Range tEts
constrains elements in block set
structure grouping operators and data used for calculation of nonlinear elastic element
OpAleLhs_dX_dX(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for material momentum.
OpAleLhs_dX_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for material momentum.
OpAleLhs_dx_dX(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for spatial momentum.
OpAleLhs_dx_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate B^T D B operator.
OpAleLhsPre_dX_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for material momentum.
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
OpAleLhsWithDensity_dX_dX(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > rho_grad_at_gauss_pts, const double rho_n, const double rho_0)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for spatial momentum.
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
OpAleLhsWithDensity_dx_dX(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > rho_grad_at_gauss_pts, const double rho_n, const double rho_0)
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
MoFEMErrorCode iNtegrate(EntData &row_data)
OpAleRhs_dX(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
OpAleRhs_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode iNtegrate(EntData &row_data)
MoFEMErrorCode iNtegrate(EntData &row_data)
OpAnalyticalInternalAleStrain_dX(const std::string row_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts, StrainFunction strain_fun, boost::shared_ptr< MatrixDouble > mat_pos_at_pts_ptr)
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunction
OpAnalyticalInternalAleStrain_dx(const std::string row_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts, StrainFunction strain_fun, boost::shared_ptr< MatrixDouble > mat_pos_at_pts_ptr)
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
MoFEMErrorCode iNtegrate(EntData &row_data)
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunction
OpAnalyticalInternalStrain_dx(const std::string row_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, StrainFunction strain_fun)
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > &t_coords) > StrainFunction
MoFEMErrorCode iNtegrate(EntData &row_data)
int nbCols
number if dof on column
bool isDiag
true if this block is on diagonal
VectorDouble nF
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Assemble local entity block matrix.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Do calculations for give operator.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
int nbIntegrationPts
number of integration points
int nbRows
number of dofs on rows
MatrixDouble transK
VectorInt colIndices
virtual MoFEMErrorCode iNtegrate(EntData &row_data)
VectorInt rowIndices
MoFEMErrorCode aSsemble(EntData &row_data)
Assemble local entity right-hand vector.
OpAssemble(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, const char type, bool symm=false)
MatrixDouble K
OpCalculateEnergy(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, SmartPetscObj< Vec > ghost_vec=SmartPetscObj< Vec >())
SmartPetscObj< Vec > ghostVec
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
OpCalculateEshelbyStress(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
OpCalculateHomogeneousStiffness(const std::string row_field, const std::string col_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr, boost::shared_ptr< DataAtIntegrationPts > data_at_pts)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Assemble mass matrix for elastic element TODO: CHANGE FORMULA *.
OpCalculateMassMatrix(const std::string row_field, const std::string col_field, BlockData &data, MassBlockData &mass_data, boost::shared_ptr< DataAtIntegrationPts > &common_data, bool symm=true)
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MassBlockData & massData
boost::shared_ptr< DataAtIntegrationPts > commonData
const double rHo0
p_0 reference density in E(p) = E * (p / p_0)^n
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
OpCalculateStiffnessScaledByDensityField(const std::string row_field, const std::string col_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, const double rho_n, const double rho_0)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
const double rhoN
exponent n in E(p) = E * (p / p_0)^n
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateStrainAle(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
OpCalculateStrain(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateStress(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
OpLhs_dx_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate B^T D B operator.
bool isALE
std::vector< EntityHandle > & mapGaussPts
map< int, BlockData > & blockSetsPtr
bool isFieldDisp
moab::Interface & postProcMesh
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
OpPostProcHookeElement(const string row_field, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, map< int, BlockData > &block_sets_ptr, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, bool is_ale=false, bool is_field_disp=true)
OpRhs_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
MoFEMErrorCode iNtegrate(EntData &row_data)