v0.15.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 */
213 : public VolumeElementForcesAndSourcesCore::UserDataOperator {
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
269 FTensor::Index<'i', 3> i;
270 FTensor::Index<'j', 3> j;
271 FTensor::Index<'k', 3> k;
272 FTensor::Index<'l', 3> l;
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 */
398 MoFEMErrorCode aSsemble(EntData &row_data);
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:
407 MoFEMErrorCode iNtegrate(EntData &row_data);
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 */
422 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
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:
431 MoFEMErrorCode iNtegrate(EntData &row_data);
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 */
446 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
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 */
461 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
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 */
485 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
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 */
509 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
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:
518 MoFEMErrorCode iNtegrate(EntData &row_data);
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 */
533 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
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 */
560 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
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:
582 MoFEMErrorCode iNtegrate(EntData &row_data);
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:
606 MoFEMErrorCode iNtegrate(EntData &row_data);
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:
631 MoFEMErrorCode iNtegrate(EntData &row_data);
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
653 MoFEMErrorCode doWork(int side, EntityType type,
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) {
700 FTensor::Index<'i', 3> i;
701 FTensor::Index<'j', 3> j;
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
744 FTensor::Index<'i', 3> i;
745 FTensor::Index<'j', 3> j;
746 FTensor::Index<'k', 3> k;
747 FTensor::Index<'l', 3> l;
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
852 FTensor::Index<'i', 3> i;
853 FTensor::Index<'j', 3> j;
854 FTensor::Index<'k', 3> k;
855 FTensor::Index<'l', 3> l;
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
904 FTensor::Index<'i', 3> i;
905 FTensor::Index<'j', 3> j;
906 FTensor::Index<'k', 3> k;
907 FTensor::Index<'l', 3> l;
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
942 FTensor::Index<'i', 3> i;
943 FTensor::Index<'j', 3> j;
944 FTensor::Index<'k', 3> k;
945 FTensor::Index<'l', 3> l;
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
1032 FTensor::Index<'i', 3> i;
1033 FTensor::Index<'j', 3> j;
1034 FTensor::Index<'k', 3> k;
1035 FTensor::Index<'l', 3> l;
1036 FTensor::Index<'m', 3> m;
1037 FTensor::Index<'n', 3> n;
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
1161 FTensor::Index<'i', 3> i;
1162 FTensor::Index<'j', 3> j;
1163 FTensor::Index<'k', 3> k;
1164 FTensor::Index<'l', 3> l;
1165 FTensor::Index<'m', 3> m;
1166 FTensor::Index<'n', 3> n;
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
1336 FTensor::Index<'i', 3> i;
1337 FTensor::Index<'j', 3> j;
1338 FTensor::Index<'k', 3> k;
1339 FTensor::Index<'l', 3> l;
1340 FTensor::Index<'m', 3> m;
1341 FTensor::Index<'n', 3> n;
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
1429 FTensor::Index<'i', 3> i;
1430 FTensor::Index<'j', 3> j;
1431 FTensor::Index<'k', 3> k;
1432 FTensor::Index<'l', 3> l;
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 int block_id = -1;
1459 for (auto &m : (blockSetsPtr)) {
1460 if (m.second.tEts.find(ent_3d) != m.second.tEts.end()) {
1461 const double young = m.second.E;
1462 const double poisson = m.second.PoissonRatio;
1463 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
1464 block_id = m.second.iD;
1465
1466 t_D(i, j, k, l) = 0.;
1467 t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
1468 t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
1469 0.5 * (1 - 2 * poisson);
1470 t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
1471 t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
1472 t_D(i, j, k, l) *= coefficient;
1473
1474 found_block = true;
1475 break;
1476 }
1477 }
1478 if (!found_block)
1479 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1480 "Element not found in any of material blocksets");
1481
1482 int def_val_int = 0;
1483 Tag tag_mat;
1484 CHKERR postProcMesh.tag_get_handle("MAT_ELASTIC", 1, MB_TYPE_INTEGER, tag_mat,
1485 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
1486 double detH = 0.;
1490 FTensor::Tensor2<double, 3, 3> t_small_strain;
1492 FTensor::Tensor2_symmetric<double, 3> t_small_strain_symm;
1493
1494 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1495
1496 if (isFieldDisp) {
1497 t_h(0, 0) += 1;
1498 t_h(1, 1) += 1;
1499 t_h(2, 2) += 1;
1500 }
1501
1502 if (!isALE) {
1503 t_small_strain_symm(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
1504 } else {
1505 CHKERR determinantTensor3by3(t_H, detH);
1506 CHKERR invertTensor3by3(t_H, detH, t_invH);
1507 t_F(i, j) = t_h(i, k) * t_invH(k, j);
1508 t_small_strain_symm(i, j) = (t_F(i, j) || t_F(j, i)) / 2.;
1509 ++t_H;
1510 }
1511
1512 t_small_strain_symm(0, 0) -= 1;
1513 t_small_strain_symm(1, 1) -= 1;
1514 t_small_strain_symm(2, 2) -= 1;
1515
1516 // symmetric tensors need improvement
1517 t_stress_symm(i, j) = t_D(i, j, k, l) * t_small_strain_symm(k, l);
1518 tensor_to_tensor(t_stress_symm, t_stress);
1519
1520 const double psi = 0.5 * t_stress_symm(i, j) * t_small_strain_symm(i, j);
1521
1522 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
1523 CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
1524 &t_stress(0, 0));
1525 CHKERR postProcMesh.tag_set_data(tag_mat, &mapGaussPts[gg], 1,
1526 &block_id);
1527
1528 ++t_h;
1529 }
1530
1532}
1533
1534template <int S>
1535HookeElement::OpAnalyticalInternalStrain_dx<S>::OpAnalyticalInternalStrain_dx(
1536 const std::string row_field,
1537 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1538 StrainFunction strain_fun)
1539 : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1540 strainFun(strain_fun) {}
1541
1542template <int S>
1544HookeElement::OpAnalyticalInternalStrain_dx<S>::iNtegrate(EntData &row_data) {
1545 FTensor::Index<'i', 3> i;
1546 FTensor::Index<'j', 3> j;
1547 FTensor::Index<'k', 3> k;
1548 FTensor::Index<'l', 3> l;
1550
1551 auto get_tensor1 = [](VectorDouble &v, const int r) {
1553 &v(r + 0), &v(r + 1), &v(r + 2));
1554 };
1555
1556 const int nb_integration_pts = getGaussPts().size2();
1557 auto t_coords = getFTensor1CoordsAtGaussPts();
1558
1559 // get element volume
1560 double vol = getVolume();
1561 auto t_w = getFTensor0IntegrationWeight();
1562
1563 nF.resize(nbRows, false);
1564 nF.clear();
1565
1566 // elastic stiffness tensor (4th rank tensor with minor and major
1567 // symmetry)
1569 MAT_TO_DDG(dataAtPts->stiffnessMat));
1570
1571 // get derivatives of base functions on rows
1572 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1573 const int row_nb_base_fun = row_data.getN().size2();
1574
1575 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1576
1577 auto t_fun_strain = strainFun(t_coords);
1579 t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1580
1581 // calculate scalar weight times element volume
1582 double a = t_w * vol;
1583
1584 auto t_nf = get_tensor1(nF, 0);
1585
1586 int rr = 0;
1587 for (; rr != nbRows / 3; ++rr) {
1588 t_nf(i) += a * t_row_diff_base(j) * t_stress(i, j);
1589 ++t_row_diff_base;
1590 ++t_nf;
1591 }
1592
1593 for (; rr != row_nb_base_fun; ++rr)
1594 ++t_row_diff_base;
1595
1596 ++t_w;
1597 ++t_coords;
1598 ++t_D;
1599 }
1600
1602}
1603
1604template <int S>
1605HookeElement::OpAnalyticalInternalAleStrain_dX<S>::
1606 OpAnalyticalInternalAleStrain_dX(
1607 const std::string row_field,
1608 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1609 StrainFunction strain_fun,
1610 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1611 : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1612 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1613
1614template <int S>
1615MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dX<S>::iNtegrate(
1616 EntData &row_data) {
1617 FTensor::Index<'i', 3> i;
1618 FTensor::Index<'j', 3> j;
1619 FTensor::Index<'k', 3> k;
1620 FTensor::Index<'l', 3> l;
1622
1623 auto get_tensor1 = [](VectorDouble &v, const int r) {
1625 &v(r + 0), &v(r + 1), &v(r + 2));
1626 };
1627
1628 const int nb_integration_pts = getGaussPts().size2();
1629
1630 auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1631 auto t_coords = get_coords();
1632
1633 // get element volume
1634 double vol = getVolume();
1635 auto t_w = getFTensor0IntegrationWeight();
1636
1637 nF.resize(nbRows, false);
1638 nF.clear();
1639
1640 // elastic stiffness tensor (4th rank tensor with minor and major
1641 // symmetry)
1643 MAT_TO_DDG(dataAtPts->stiffnessMat));
1644 auto t_F = getFTensor2FromMat<3, 3>(*(dataAtPts->FMat));
1645 auto &det_H = *dataAtPts->detHVec;
1646 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1647
1648 // get derivatives of base functions on rows
1649 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1650 const int row_nb_base_fun = row_data.getN().size2();
1651
1652 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1653
1654 auto t_fun_strain = strainFun(t_coords);
1656 t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1657 FTensor::Tensor2<double, 3, 3> t_eshelby_stress;
1658 t_eshelby_stress(i, j) = -t_F(k, i) * t_stress(k, j);
1659
1660 // calculate scalar weight times element volume
1661 double a = t_w * vol * det_H[gg];
1662
1663 auto t_nf = get_tensor1(nF, 0);
1664
1665 int rr = 0;
1666 for (; rr != nbRows / 3; ++rr) {
1667 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1668 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1669 t_nf(i) += a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1670 ++t_row_diff_base;
1671 ++t_nf;
1672 }
1673
1674 for (; rr != row_nb_base_fun; ++rr)
1675 ++t_row_diff_base;
1676
1677 ++t_w;
1678 ++t_coords;
1679 ++t_F;
1680 ++t_invH;
1681 ++t_D;
1682 }
1683
1685}
1686
1687template <int S>
1688HookeElement::OpAnalyticalInternalAleStrain_dx<S>::
1689 OpAnalyticalInternalAleStrain_dx(
1690 const std::string row_field,
1691 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1692 StrainFunction strain_fun,
1693 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1694 : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1695 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1696
1697template <int S>
1698MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dx<S>::iNtegrate(
1699 EntData &row_data) {
1700 FTensor::Index<'i', 3> i;
1701 FTensor::Index<'j', 3> j;
1702 FTensor::Index<'k', 3> k;
1703 FTensor::Index<'l', 3> l;
1705
1706 auto get_tensor1 = [](VectorDouble &v, const int r) {
1708 &v(r + 0), &v(r + 1), &v(r + 2));
1709 };
1710
1711 const int nb_integration_pts = getGaussPts().size2();
1712
1713 auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1714 auto t_coords = get_coords();
1715
1716 // get element volume
1717 double vol = getVolume();
1718 auto t_w = getFTensor0IntegrationWeight();
1719
1720 nF.resize(nbRows, false);
1721 nF.clear();
1722
1723 // elastic stiffness tensor (4th rank tensor with minor and major
1724 // symmetry)
1726 MAT_TO_DDG(dataAtPts->stiffnessMat));
1727 auto &det_H = *dataAtPts->detHVec;
1728 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1729
1730 // get derivatives of base functions on rows
1731 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1732 const int row_nb_base_fun = row_data.getN().size2();
1733
1734 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1735
1736 auto t_fun_strain = strainFun(t_coords);
1738 t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1739
1740 // calculate scalar weight times element volume
1741 double a = t_w * vol * det_H[gg];
1742
1743 auto t_nf = get_tensor1(nF, 0);
1744
1745 int rr = 0;
1746 for (; rr != nbRows / 3; ++rr) {
1747 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1748 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1749 t_nf(i) += a * t_row_diff_base_pulled(j) * t_stress(i, j);
1750 ++t_row_diff_base;
1751 ++t_nf;
1752 }
1753
1754 for (; rr != row_nb_base_fun; ++rr)
1755 ++t_row_diff_base;
1756
1757 ++t_w;
1758 ++t_coords;
1759 ++t_invH;
1760 ++t_D;
1761 }
1762
1764}
1765
1766#endif // __HOOKE_ELEMENT_HPP
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 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 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()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
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)
const double n
refractive index of diffusive medium
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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
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.
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
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.
FTensor::Index< 'm', 3 > m
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
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
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
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
OpAleLhs_dX_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunction
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunction
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > &t_coords) > StrainFunction
int nbCols
number if dof on column
bool isDiag
true if this block is on diagonal
VectorDouble nF
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
int nbIntegrationPts
number of integration points
int nbRows
number of dofs on rows
MatrixDouble transK
VectorInt colIndices
VectorInt rowIndices
MatrixDouble K
SmartPetscObj< Vec > ghostVec
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
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
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
const double rhoN
exponent n in E(p) = E * (p / p_0)^n
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
std::vector< EntityHandle > & mapGaussPts
map< int, BlockData > & blockSetsPtr
moab::Interface & postProcMesh
boost::shared_ptr< DataAtIntegrationPts > dataAtPts