v0.13.1
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 * This file is part of MoFEM.
18 * MoFEM is free software: you can redistribute it and/or modify it under
19 * the terms of the GNU Lesser General Public License as published by the
20 * Free Software Foundation, either version 3 of the License, or (at your
21 * option) any later version.
22 *
23 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
24 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
25 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
26 * License for more details.
27 *
28 * You should have received a copy of the GNU Lesser General Public
29 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
30
31#ifndef __HOOKE_ELEMENT_HPP
32#define __HOOKE_ELEMENT_HPP
33
34#ifndef __BASICFINITEELEMENTS_HPP__
36#endif // __BASICFINITEELEMENTS_HPP__
37
38#ifndef __NONLINEAR_ELASTIC_HPP
39
41
42 /** \brief data for calculation heat conductivity and heat capacity elements
43 * \ingroup nonlinear_elastic_elem
44 */
45 struct BlockData {
46 int iD;
47 double E;
49 Range tEts; ///< constrains elements in block set
52 };
53};
54
55#endif // __NONLINEAR_ELASTIC_HPP
56
57/** \brief structure grouping operators and data used for calculation of
58 * nonlinear elastic element \ingroup nonlinear_elastic_elem
59 *
60 * In order to assemble matrices and right hand vectors, the loops over
61 * elements, entities over that elements and finally loop over integration
62 * points are executed.
63 *
64 * Following implementation separate those three categories of loops and to each
65 * loop attach operator.
66 *
67 */
68
69#ifndef __CONVECTIVE_MASS_ELEMENT_HPP
71 /** \brief data for calculation inertia forces
72 * \ingroup user_modules
73 */
74 struct BlockData {
75 double rho0; ///< reference density
76 VectorDouble a0; ///< constant acceleration
77 Range tEts; ///< elements in block set
78 };
79}
80
81#endif //__CONVECTIVE_MASS_ELEMENT_HPP
82struct HookeElement {
83
86
91
93
94 boost::shared_ptr<MatrixDouble> smallStrainMat;
95 boost::shared_ptr<MatrixDouble> hMat;
96 boost::shared_ptr<MatrixDouble> FMat;
97
98 boost::shared_ptr<MatrixDouble> HMat;
99 boost::shared_ptr<VectorDouble> detHVec;
100 boost::shared_ptr<MatrixDouble> invHMat;
101
102 boost::shared_ptr<MatrixDouble> cauchyStressMat;
103 boost::shared_ptr<MatrixDouble> stiffnessMat;
104 boost::shared_ptr<VectorDouble> energyVec;
105 boost::shared_ptr<MatrixDouble> eshelbyStressMat;
106
107 boost::shared_ptr<MatrixDouble> eshelbyStress_dx;
108
110
111 smallStrainMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
112 hMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
113 FMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
114
115 HMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
116 detHVec = boost::shared_ptr<VectorDouble>(new VectorDouble());
117 invHMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
118
119 cauchyStressMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
120 stiffnessMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
121 energyVec = boost::shared_ptr<VectorDouble>(new VectorDouble());
122 eshelbyStressMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
123
124 eshelbyStress_dx = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
125 }
126
129 };
130
131 template <bool D = false>
133
134 OpCalculateStrain(const std::string row_field, const std::string col_field,
135 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
136
137 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
138
139 private:
140 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
141 };
142
144
145 OpCalculateStrainAle(const std::string row_field,
146 const std::string col_field,
147 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
148
149 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
150
151 private:
152 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
153 };
154
155#define MAT_TO_DDG(SM) \
156 &(*SM)(0, 0), &(*SM)(1, 0), &(*SM)(2, 0), &(*SM)(3, 0), &(*SM)(4, 0), \
157 &(*SM)(5, 0), &(*SM)(6, 0), &(*SM)(7, 0), &(*SM)(8, 0), &(*SM)(9, 0), \
158 &(*SM)(10, 0), &(*SM)(11, 0), &(*SM)(12, 0), &(*SM)(13, 0), \
159 &(*SM)(14, 0), &(*SM)(15, 0), &(*SM)(16, 0), &(*SM)(17, 0), \
160 &(*SM)(18, 0), &(*SM)(19, 0), &(*SM)(20, 0), &(*SM)(21, 0), \
161 &(*SM)(22, 0), &(*SM)(23, 0), &(*SM)(24, 0), &(*SM)(25, 0), \
162 &(*SM)(26, 0), &(*SM)(27, 0), &(*SM)(28, 0), &(*SM)(29, 0), \
163 &(*SM)(30, 0), &(*SM)(31, 0), &(*SM)(32, 0), &(*SM)(33, 0), \
164 &(*SM)(34, 0), &(*SM)(35, 0)
165
166 template <int S = 0> struct OpCalculateStress : public VolUserDataOperator {
167
168 OpCalculateStress(const std::string row_field, const std::string col_field,
169 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
170
171 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
172
173 protected:
174 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
175 };
176
178
179 OpCalculateEnergy(const std::string row_field, const std::string col_field,
180 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;
188 };
189
191
193 const std::string row_field, const std::string col_field,
194 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
195
196 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
197
198 protected:
199 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
200 };
201
202 template <int S = 0>
204
206 const std::string row_field, const std::string col_field,
207 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
208 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
209
210 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
211
212 protected:
213 boost::shared_ptr<map<int, BlockData>>
214 blockSetsPtr; ///< Structure keeping data about problem, like
215 ///< material parameters
216 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
217 };
218
219 /** * @brief Assemble mass matrix for elastic element TODO: CHANGE FORMULA *
220 * \f[
221 * {\bf{M}} = \int\limits_\Omega
222 * \f]
223 *
224 */
227
232
233 boost::shared_ptr<DataAtIntegrationPts> commonData;
234
235 OpCalculateMassMatrix(const std::string row_field,
236 const std::string col_field, BlockData &data,
237 MassBlockData &mass_data,
238 boost::shared_ptr<DataAtIntegrationPts> &common_data,
239 bool symm = true)
241 row_field, col_field, OPROWCOL, symm),
242 commonData(common_data), dAta(data), massData(mass_data) {}
243
244 PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
245 EntityType col_type,
247 EntitiesFieldData::EntData &col_data) {
249
250 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
252 &m(3 * r + 0, 3 * c + 0), &m(3 * r + 0, 3 * c + 1),
253 &m(3 * r + 0, 3 * c + 2), &m(3 * r + 1, 3 * c + 0),
254 &m(3 * r + 1, 3 * c + 1), &m(3 * r + 1, 3 * c + 2),
255 &m(3 * r + 2, 3 * c + 0), &m(3 * r + 2, 3 * c + 1),
256 &m(3 * r + 2, 3 * c + 2));
257 };
258
259 const int row_nb_dofs = row_data.getIndices().size();
260 if (!row_nb_dofs)
262 const int col_nb_dofs = col_data.getIndices().size();
263 if (!col_nb_dofs)
265 if (dAta.tEts.find(getFEEntityHandle()) == dAta.tEts.end()) {
267 }
268 if (massData.tEts.find(getFEEntityHandle()) == massData.tEts.end()) {
270 }
271
272 const bool diagonal_block =
273 (row_type == col_type) && (row_side == col_side);
274 // get number of integration points
275 // Set size can clear local tangent matrix
276 locK.resize(row_nb_dofs, col_nb_dofs, false);
277 locK.clear();
278
279 const int row_nb_gauss_pts = row_data.getN().size1();
280 const int row_nb_base_functions = row_data.getN().size2();
281
286
287 double density = massData.rho0;
288
289 // get integration weights
290 auto t_w = getFTensor0IntegrationWeight();
291
292 // integrate local matrix for entity block
293 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
294
295 auto t_row_base_func = row_data.getFTensor0N(gg, 0);
296
297 // Get volume and integration weight
298 double w = getVolume() * t_w;
299
300 for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
301 auto t_col_base_func = col_data.getFTensor0N(gg, 0);
302 for (int col_bb = 0; col_bb != col_nb_dofs / 3; col_bb++) {
303 auto t_assemble = get_tensor2(locK, row_bb, col_bb);
304 t_assemble(i, j) += density * t_row_base_func * t_col_base_func * w;
305 // Next base function for column
306 ++t_col_base_func;
307 }
308 // Next base function for row
309 ++t_row_base_func;
310 }
311 // Next integration point for getting weight
312 ++t_w;
313 }
314
315 CHKERR MatSetValues(getKSPB(), row_data, col_data, &locK(0, 0),
316 ADD_VALUES);
317
318 // is symmetric
319 if (row_type != col_type || row_side != col_side) {
320 translocK.resize(col_nb_dofs, row_nb_dofs, false);
321 noalias(translocK) = trans(locK);
322
323 CHKERR MatSetValues(getKSPB(), col_data, row_data, &translocK(0, 0),
324 ADD_VALUES);
325 }
326
328 }
329 };
330
332 protected:
333 boost::shared_ptr<map<int, BlockData>>
334 blockSetsPtr; ///< Structure keeping data about problem, like
335 ///< material parameters
336 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
337
338 boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
339 const double rhoN; ///< exponent n in E(p) = E * (p / p_0)^n
340 const double rHo0; ///< p_0 reference density in E(p) = E * (p / p_0)^n
341 // // where p is density, E - youngs modulus
342 public:
344 const std::string row_field, const std::string col_field,
345 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
346 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
347 boost::shared_ptr<VectorDouble> rho_at_gauss_pts, const double rho_n,
348 const double rho_0);
349
350 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
351 };
352
354
355 OpAssemble(const std::string row_field, const std::string col_field,
356 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
357 const char type, bool symm = false);
358
359 /**
360 * \brief Do calculations for give operator
361 * @param row_side row side number (local number) of entity on element
362 * @param col_side column side number (local number) of entity on element
363 * @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
364 * @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
365 * @param row_data data for row
366 * @param col_data data for column
367 * @return error code
368 */
369 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
370 EntityType col_type, EntData &row_data,
371 EntData &col_data);
372
373 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
374
375 protected:
376 // Finite element stiffness sub-matrix K_ij
380
381 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
382
385
386 int nbRows; ///< number of dofs on rows
387 int nbCols; ///< number if dof on column
388 int nbIntegrationPts; ///< number of integration points
389 bool isDiag; ///< true if this block is on diagonal
390
391 virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
392
393 virtual MoFEMErrorCode iNtegrate(EntData &row_data);
394
395 /**
396 * \brief Assemble local entity block matrix
397 * @param row_data row data (consist base functions on row entity)
398 * @param col_data column data (consist base functions on column
399 * entity)
400 * @return error code
401 */
402 MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data);
403
404 /**
405 * \brief Assemble local entity right-hand vector
406 * @param row_data row data (consist base functions on row entity)
407 * @param col_data column data (consist base functions on column
408 * entity)
409 * @return error code
410 */
412 };
413
414 struct OpRhs_dx : public OpAssemble {
415
416 OpRhs_dx(const std::string row_field, const std::string col_field,
417 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
418
419 protected:
421 };
422
423 template <int S = 0> struct OpLhs_dx_dx : public OpAssemble {
424
425 OpLhs_dx_dx(const std::string row_field, const std::string col_field,
426 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
427
428 protected:
429 /**
430 * \brief Integrate B^T D B operator
431 * @param row_data row data (consist base functions on row entity)
432 * @param col_data column data (consist base functions on column entity)
433 * @return error code
434 */
436 };
437
438 struct OpAleRhs_dx : public OpAssemble {
439
440 OpAleRhs_dx(const std::string row_field, const std::string col_field,
441 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
442
443 protected:
445 };
446
447 template <int S = 0> struct OpAleLhs_dx_dx : public OpAssemble {
448
449 OpAleLhs_dx_dx(const std::string row_field, const std::string col_field,
450 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
451
452 protected:
453 /**
454 * \brief Integrate B^T D B operator
455 * @param row_data row data (consist base functions on row entity)
456 * @param col_data column data (consist base functions on column entity)
457 * @return error code
458 */
460 };
461
462 template <int S = 0> struct OpAleLhs_dx_dX : public OpAssemble {
463
464 OpAleLhs_dx_dX(const std::string row_field, const std::string col_field,
465 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
466
467 protected:
468 /**
469 * \brief Integrate tangent stiffness for spatial momentum
470 * @param row_data row data (consist base functions on row entity)
471 * @param col_data column data (consist base functions on column entity)
472 * @return error code
473 */
475 };
476
478
479 boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
480 boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
481 const double rhoN;
482 const double rHo0;
483
485 const std::string row_field, const std::string col_field,
486 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
487 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
488 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
489 const double rho_n, const double rho_0);
490
491 protected:
492 /**
493 * \brief Integrate tangent stiffness for spatial momentum
494 * @param row_data row data (consist base functions on row entity)
495 * @param col_data column data (consist base functions on column entity)
496 * @return error code
497 */
499 };
500
502
503 boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
504 boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
505 const double rhoN;
506 const double rHo0;
507
509 const std::string row_field, const std::string col_field,
510 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
511 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
512 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
513 const double rho_n, const double rho_0);
514
515 protected:
516 /**
517 * \brief Integrate tangent stiffness for material momentum
518 * @param row_data row data (consist base functions on row entity)
519 * @param col_data column data (consist base functions on column entity)
520 * @return error code
521 */
523 };
524
525 struct OpAleRhs_dX : public OpAssemble {
526
527 OpAleRhs_dX(const std::string row_field, const std::string col_field,
528 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
529
530 protected:
532 };
533
534 template <int S = 0> struct OpAleLhs_dX_dX : public OpAssemble {
535
536 OpAleLhs_dX_dX(const std::string row_field, const std::string col_field,
537 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
538
539 protected:
540 /**
541 * \brief Integrate tangent stiffness for material momentum
542 * @param row_data row data (consist base functions on row entity)
543 * @param col_data column data (consist base functions on column entity)
544 * @return error code
545 */
547 };
548
549 template <int S = 0> struct OpAleLhsPre_dX_dx : public VolUserDataOperator {
550
551 OpAleLhsPre_dX_dx(const std::string row_field, const std::string col_field,
552 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
553
554 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
555
556 private:
557 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
558 };
559
560 struct OpAleLhs_dX_dx : public OpAssemble {
561
562 OpAleLhs_dX_dx(const std::string row_field, const std::string col_field,
563 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
564 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
565
566 protected:
567 /**
568 * \brief Integrate tangent stiffness for material momentum
569 * @param row_data row data (consist base functions on row entity)
570 * @param col_data column data (consist base functions on column entity)
571 * @return error code
572 */
574 };
575
576 template <int S> struct OpAnalyticalInternalStrain_dx : public OpAssemble {
577
578 typedef boost::function<
579
581
583
584 )
585
586 >
588
590 const std::string row_field,
591 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
592 StrainFunction strain_fun);
593
594 protected:
597 };
598
599 template <int S> struct OpAnalyticalInternalAleStrain_dX : public OpAssemble {
600
601 typedef boost::function<
602
604
606
607 )
608
609 >
611
613 const std::string row_field,
614 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
615 StrainFunction strain_fun,
616 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
617
618 protected:
621 boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
622 };
623
624 template <int S> struct OpAnalyticalInternalAleStrain_dx : public OpAssemble {
625
626 typedef boost::function<
627
629
631
632 )
633
634 >
636
638 const std::string row_field,
639 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
640 StrainFunction strain_fun,
641 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
642
643 protected:
646 boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
647 };
648
649 template <class ELEMENT>
651 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
652 map<int, BlockData>
653 &blockSetsPtr; // FIXME: (works only with the first block)
655 std::vector<EntityHandle> &mapGaussPts;
656 bool isALE;
658
659 OpPostProcHookeElement(const string row_field,
660 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
661 map<int, BlockData> &block_sets_ptr,
662 moab::Interface &post_proc_mesh,
663 std::vector<EntityHandle> &map_gauss_pts,
664 bool is_ale = false, bool is_field_disp = true);
665
668 };
669
670 static MoFEMErrorCode
672 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr);
673
674 static MoFEMErrorCode
676 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
677 const std::string element_name, const std::string x_field,
678 const std::string X_field, const bool ale);
679
680 static MoFEMErrorCode
681 setOperators(boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr,
682 boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr,
683 boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
684 const std::string x_field, const std::string X_field,
685 const bool ale, const bool field_disp,
686 const EntityType type = MBTET,
687 boost::shared_ptr<DataAtIntegrationPts> data_at_pts = nullptr);
688
689 static MoFEMErrorCode
690 calculateEnergy(DM dm, boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
691 const std::string x_field, const std::string X_field,
692 const bool ale, const bool field_disp,
693 SmartPetscObj<Vec> &v_energy_ptr);
694
695private:
697};
698
699template <bool D>
700HookeElement::OpCalculateStrain<D>::OpCalculateStrain(
701 const std::string row_field, const std::string col_field,
702 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
703 : VolUserDataOperator(row_field, col_field, OPROW, true),
704 dataAtPts(data_at_pts) {
705 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
706}
707
708template <bool D>
709MoFEMErrorCode HookeElement::OpCalculateStrain<D>::doWork(int row_side,
710 EntityType row_type,
711 EntData &row_data) {
715 // get number of integration points
716 const int nb_integration_pts = getGaussPts().size2();
717 dataAtPts->smallStrainMat->resize(6, nb_integration_pts, false);
718 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
719 auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hMat));
720
721 for (int gg = 0; gg != nb_integration_pts; ++gg) {
722 t_strain(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
723
724 // If displacement field, not field o spatial positons is given
725 if (!D) {
726 t_strain(0, 0) -= 1;
727 t_strain(1, 1) -= 1;
728 t_strain(2, 2) -= 1;
729 }
730
731 ++t_strain;
732 ++t_h;
733 }
735}
736
737template <int S>
738HookeElement::OpAleLhs_dx_dx<S>::OpAleLhs_dx_dx(
739 const std::string row_field, const std::string col_field,
740 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
741 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
742
743template <int S>
744MoFEMErrorCode HookeElement::OpAleLhs_dx_dx<S>::iNtegrate(EntData &row_data,
745 EntData &col_data) {
747
748 // get sub-block (3x3) of local stiffens matrix, here represented by
749 // second order tensor
750 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
752 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
753 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
754 &m(r + 2, c + 2));
755 };
756
761
762 // get element volume
763 double vol = getVolume();
764
765 // get intergrayion weights
766 auto t_w = getFTensor0IntegrationWeight();
767
768 // get derivatives of base functions on rows
769 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
770 const int row_nb_base_fun = row_data.getN().size2();
771
772 // Elastic stiffness tensor (4th rank tensor with minor and major
773 // symmetry)
775 MAT_TO_DDG(dataAtPts->stiffnessMat));
776
777 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
778 auto &det_H = *dataAtPts->detHVec;
779
780 // iterate over integration points
781 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
782
783 // calculate scalar weight times element volume
784 double a = t_w * vol * det_H[gg];
785
786 // iterate over row base functions
787 int rr = 0;
788 for (; rr != nbRows / 3; ++rr) {
789
790 // get sub matrix for the row
791 auto t_m = get_tensor2(K, 3 * rr, 0);
792
793 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
794 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
795
797 // I mix up the indices here so that it behaves like a
798 // Dg. That way I don't have to have a separate wrapper
799 // class Christof_Expr, which simplifies things.
800 t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base_pulled(i));
801
802 // get derivatives of base functions for columns
803 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
804
805 // iterate column base functions
806 for (int cc = 0; cc != nbCols / 3; ++cc) {
807
808 FTensor::Tensor1<double, 3> t_col_diff_base_pulled;
809 t_col_diff_base_pulled(j) = t_col_diff_base(i) * t_invH(i, j);
810
811 // integrate block local stiffens matrix
812 t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base_pulled(k);
813
814 // move to next column base function
815 ++t_col_diff_base;
816
817 // move to next block of local stiffens matrix
818 ++t_m;
819 }
820
821 // move to next row base function
822 ++t_row_diff_base;
823 }
824
825 for (; rr != row_nb_base_fun; ++rr)
826 ++t_row_diff_base;
827
828 // move to next integration weight
829 ++t_w;
830 ++t_D;
831 ++t_invH;
832 }
833
835}
836
837template <int S>
838HookeElement::OpCalculateHomogeneousStiffness<S>::
839 OpCalculateHomogeneousStiffness(
840 const std::string row_field, const std::string col_field,
841 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
842 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
843 : VolUserDataOperator(row_field, col_field, OPROW, true),
844 blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts) {
845 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
846}
847
848template <int S>
849MoFEMErrorCode HookeElement::OpCalculateHomogeneousStiffness<S>::doWork(
850 int row_side, EntityType row_type, EntData &row_data) {
852
853 for (auto &m : (*blockSetsPtr)) {
854 if (m.second.tEts.find(getFEEntityHandle()) != m.second.tEts.end()) {
855
856 dataAtPts->stiffnessMat->resize(36, 1, false);
858 MAT_TO_DDG(dataAtPts->stiffnessMat));
859 const double young = m.second.E;
860 const double poisson = m.second.PoissonRatio;
861
862 // coefficient used in intermediate calculation
863 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
864
869
870 t_D(i, j, k, l) = 0.;
871
872 t_D(0, 0, 0, 0) = 1 - poisson;
873 t_D(1, 1, 1, 1) = 1 - poisson;
874 t_D(2, 2, 2, 2) = 1 - poisson;
875
876 t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * poisson);
877 t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * poisson);
878 t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * poisson);
879
880 t_D(0, 0, 1, 1) = poisson;
881 t_D(1, 1, 0, 0) = poisson;
882 t_D(0, 0, 2, 2) = poisson;
883 t_D(2, 2, 0, 0) = poisson;
884 t_D(1, 1, 2, 2) = poisson;
885 t_D(2, 2, 1, 1) = poisson;
886
887 t_D(i, j, k, l) *= coefficient;
888
889 break;
890 }
891 }
892
894}
895
896template <int S>
897HookeElement::OpCalculateStress<S>::OpCalculateStress(
898 const std::string row_field, const std::string col_field,
899 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
900 : VolUserDataOperator(row_field, col_field, OPROW, true),
901 dataAtPts(data_at_pts) {
902 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
903}
904
905template <int S>
906MoFEMErrorCode HookeElement::OpCalculateStress<S>::doWork(int row_side,
907 EntityType row_type,
908 EntData &row_data) {
910 // get number of integration points
911 const int nb_integration_pts = getGaussPts().size2();
912 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
913 dataAtPts->cauchyStressMat->resize(6, nb_integration_pts, false);
914 auto t_cauchy_stress =
915 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
916
921
922 // elastic stiffness tensor (4th rank tensor with minor and major
923 // symmetry)
925 MAT_TO_DDG(dataAtPts->stiffnessMat));
926 for (int gg = 0; gg != nb_integration_pts; ++gg) {
927 t_cauchy_stress(i, j) = t_D(i, j, k, l) * t_strain(k, l);
928 ++t_strain;
929 ++t_cauchy_stress;
930 ++t_D;
931 }
933}
934
935template <int S>
936HookeElement::OpLhs_dx_dx<S>::OpLhs_dx_dx(
937 const std::string row_field, const std::string col_field,
938 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
939 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
940
941template <int S>
942MoFEMErrorCode HookeElement::OpLhs_dx_dx<S>::iNtegrate(EntData &row_data,
943 EntData &col_data) {
945
946 // get sub-block (3x3) of local stiffens matrix, here represented by
947 // second order tensor
948 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
950 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
951 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
952 &m(r + 2, c + 2));
953 };
954
959
960 // get element volume
961 double vol = getVolume();
962
963 // get intergrayion weights
964 auto t_w = getFTensor0IntegrationWeight();
965
966 // get derivatives of base functions on rows
967 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
968 const int row_nb_base_fun = row_data.getN().size2();
969
970 // Elastic stiffness tensor (4th rank tensor with minor and major
971 // symmetry)
973 MAT_TO_DDG(dataAtPts->stiffnessMat));
974
975 // iterate over integration points
976 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
977
978 // calculate scalar weight times element volume
979 double a = t_w * vol;
980
981 // iterate over row base functions
982 int rr = 0;
983 for (; rr != nbRows / 3; ++rr) {
984
985 // get sub matrix for the row
986 auto t_m = get_tensor2(K, 3 * rr, 0);
987
988 // get derivatives of base functions for columns
989 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
990
992 // I mix up the indices here so that it behaves like a
993 // Dg. That way I don't have to have a separate wrapper
994 // class Christof_Expr, which simplifies things.
995 t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
996
997 // iterate column base functions
998 for (int cc = 0; cc != nbCols / 3; ++cc) {
999
1000 // integrate block local stiffens matrix
1001 t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
1002
1003 // move to next column base function
1004 ++t_col_diff_base;
1005
1006 // move to next block of local stiffens matrix
1007 ++t_m;
1008 }
1009
1010 // move to next row base function
1011 ++t_row_diff_base;
1012 }
1013
1014 for (; rr != row_nb_base_fun; ++rr)
1015 ++t_row_diff_base;
1016
1017 // move to next integration weight
1018 ++t_w;
1019 ++t_D;
1020 }
1021
1023}
1024
1025template <int S>
1026HookeElement::OpAleLhs_dx_dX<S>::OpAleLhs_dx_dX(
1027 const std::string row_field, const std::string col_field,
1028 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1029 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
1030
1031template <int S>
1032MoFEMErrorCode HookeElement::OpAleLhs_dx_dX<S>::iNtegrate(EntData &row_data,
1033 EntData &col_data) {
1035
1036 // get sub-block (3x3) of local stiffens matrix, here represented by
1037 // second order tensor
1038 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1040 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1041 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1042 &m(r + 2, c + 2));
1043 };
1044
1051
1052 // get element volume
1053 double vol = getVolume();
1054
1055 // get intergrayion weights
1056 auto t_w = getFTensor0IntegrationWeight();
1057
1058 // get derivatives of base functions on rows
1059 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1060 const int row_nb_base_fun = row_data.getN().size2();
1061
1062 // Elastic stiffness tensor (4th rank tensor with minor and major
1063 // symmetry)
1065 MAT_TO_DDG(dataAtPts->stiffnessMat));
1066
1067 auto t_cauchy_stress =
1068 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1069 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1070 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1071 auto &det_H = *dataAtPts->detHVec;
1072
1073 // iterate over integration points
1074 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1075
1076 // calculate scalar weight times element volume
1077 double a = t_w * vol * det_H[gg];
1078
1080 t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1081
1082 // iterate over row base functions
1083 int rr = 0;
1084 for (; rr != nbRows / 3; ++rr) {
1085
1086 // get sub matrix for the row
1087 auto t_m = get_tensor2(K, 3 * rr, 0);
1088
1089 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1090 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1091
1092 FTensor::Tensor1<double, 3> t_row_stress;
1093 t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_cauchy_stress(i, j);
1094
1095 FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
1096 t_row_diff_base_pulled_dX(j, k, l) =
1097 -(t_invH(i, k) * t_row_diff_base(i)) * t_invH(l, j);
1098
1099 FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
1100 t_row_dX_stress(i, k, l) =
1101 a * (t_row_diff_base_pulled_dX(j, k, l) * t_cauchy_stress(j, i));
1102
1104 t_row_D(l, j, k) = (a * t_row_diff_base_pulled(i)) * t_D(i, j, k, l);
1105
1106 FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1107 // FIXME: This operator is not implemented, doing operation by hand
1108 // t_row_stress_dX(i, m, n) = t_row_D(i, k, l) * t_F_dX(k, l, m, n);
1109 t_row_stress_dX(i, j, k) = 0;
1110 for (int ii = 0; ii != 3; ++ii)
1111 for (int mm = 0; mm != 3; ++mm)
1112 for (int nn = 0; nn != 3; ++nn) {
1113 auto &v = t_row_stress_dX(ii, mm, nn);
1114 for (int kk = 0; kk != 3; ++kk)
1115 for (int ll = 0; ll != 3; ++ll)
1116 v += t_row_D(ii, kk, ll) * t_F_dX(kk, ll, mm, nn);
1117 }
1118
1119 // get derivatives of base functions for columns
1120 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1121
1122 // iterate column base functions
1123 for (int cc = 0; cc != nbCols / 3; ++cc) {
1124
1125 t_m(i, k) += t_row_stress(i) * (t_invH(j, k) * t_col_diff_base(j));
1126 t_m(i, k) += t_row_dX_stress(i, k, l) * t_col_diff_base(l);
1127 t_m(i, k) += t_row_stress_dX(i, k, l) * t_col_diff_base(l);
1128
1129 // move to next column base function
1130 ++t_col_diff_base;
1131
1132 // move to next block of local stiffens matrix
1133 ++t_m;
1134 }
1135
1136 // move to next row base function
1137 ++t_row_diff_base;
1138 }
1139
1140 for (; rr != row_nb_base_fun; ++rr)
1141 ++t_row_diff_base;
1142
1143 // move to next integration weight
1144 ++t_w;
1145 ++t_D;
1146 ++t_cauchy_stress;
1147 ++t_invH;
1148 ++t_h;
1149 }
1150
1152}
1153
1154template <int S>
1155HookeElement::OpAleLhs_dX_dX<S>::OpAleLhs_dX_dX(
1156 const std::string row_field, const std::string col_field,
1157 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1158 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
1159
1160template <int S>
1161MoFEMErrorCode HookeElement::OpAleLhs_dX_dX<S>::iNtegrate(EntData &row_data,
1162 EntData &col_data) {
1164
1165 // get sub-block (3x3) of local stiffens matrix, here represented by
1166 // second order tensor
1167 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1169 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1170 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1171 &m(r + 2, c + 2));
1172 };
1173
1180
1181 // get element volume
1182 double vol = getVolume();
1183
1184 // get intergrayion weights
1185 auto t_w = getFTensor0IntegrationWeight();
1186
1187 // get derivatives of base functions on rows
1188 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1189 const int row_nb_base_fun = row_data.getN().size2();
1190
1191 // Elastic stiffness tensor (4th rank tensor with minor and major
1192 // symmetry)
1194 MAT_TO_DDG(dataAtPts->stiffnessMat));
1195 auto t_cauchy_stress =
1196 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1197 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
1198 auto t_eshelby_stress =
1199 getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
1200 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1201 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1202 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1203 auto &det_H = *dataAtPts->detHVec;
1204
1205 // iterate over integration points
1206 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1207
1208 // calculate scalar weight times element volume
1209 double a = t_w * vol * det_H[gg];
1210
1212 t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1213
1215 t_D_strain_dX(i, j, m, n) = 0.;
1216 for (int ii = 0; ii != 3; ++ii)
1217 for (int jj = 0; jj != 3; ++jj)
1218 for (int ll = 0; ll != 3; ++ll)
1219 for (int kk = 0; kk != 3; ++kk) {
1220 auto &v = t_D_strain_dX(ii, jj, kk, ll);
1221 for (int mm = 0; mm != 3; ++mm)
1222 for (int nn = 0; nn != 3; ++nn)
1223 v += t_D(ii, jj, mm, nn) * t_F_dX(mm, nn, kk, ll);
1224 }
1225
1226 FTensor::Tensor4<double, 3, 3, 3, 3> t_eshelby_stress_dX;
1227 t_eshelby_stress_dX(i, j, m, n) = t_F(k, i) * t_D_strain_dX(k, j, m, n);
1228
1229 for (int ii = 0; ii != 3; ++ii)
1230 for (int jj = 0; jj != 3; ++jj)
1231 for (int mm = 0; mm != 3; ++mm)
1232 for (int nn = 0; nn != 3; ++nn) {
1233 auto &v = t_eshelby_stress_dX(ii, jj, mm, nn);
1234 for (int kk = 0; kk != 3; ++kk)
1235 v += t_F_dX(kk, ii, mm, nn) * t_cauchy_stress(kk, jj);
1236 }
1237
1238 t_eshelby_stress_dX(i, j, k, l) *= -1;
1239
1241 t_energy_dX(k, l) = t_F_dX(i, j, k, l) * t_cauchy_stress(i, j);
1242 t_energy_dX(k, l) +=
1243 (t_strain(m, n) * t_D(m, n, i, j)) * t_F_dX(i, j, k, l);
1244 t_energy_dX(k, l) /= 2.;
1245
1246 for (int kk = 0; kk != 3; ++kk)
1247 for (int ll = 0; ll != 3; ++ll) {
1248 auto v = t_energy_dX(kk, ll);
1249 for (int ii = 0; ii != 3; ++ii)
1250 t_eshelby_stress_dX(ii, ii, kk, ll) += v;
1251 }
1252
1253 // iterate over row base functions
1254 int rr = 0;
1255 for (; rr != nbRows / 3; ++rr) {
1256
1257 // get sub matrix for the row
1258 auto t_m = get_tensor2(K, 3 * rr, 0);
1259
1260 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1261 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1262
1263 FTensor::Tensor1<double, 3> t_row_stress;
1264 t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1265
1266 FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
1267 t_row_diff_base_pulled_dX(j, k, l) =
1268 -(t_row_diff_base(i) * t_invH(i, k)) * t_invH(l, j);
1269
1270 FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
1271 t_row_dX_stress(i, k, l) =
1272 a * (t_row_diff_base_pulled_dX(j, k, l) * t_eshelby_stress(i, j));
1273
1274 FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1275 t_row_stress_dX(i, m, n) =
1276 a * t_row_diff_base_pulled(j) * t_eshelby_stress_dX(i, j, m, n);
1277
1278 // get derivatives of base functions for columns
1279 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1280
1281 // iterate column base functions
1282 for (int cc = 0; cc != nbCols / 3; ++cc) {
1283
1284 t_m(i, k) += t_row_stress(i) * (t_invH(j, k) * t_col_diff_base(j));
1285 t_m(i, k) += t_row_dX_stress(i, k, l) * t_col_diff_base(l);
1286 t_m(i, k) += t_row_stress_dX(i, k, l) * t_col_diff_base(l);
1287
1288 // move to next column base function
1289 ++t_col_diff_base;
1290
1291 // move to next block of local stiffens matrix
1292 ++t_m;
1293 }
1294
1295 // move to next row base function
1296 ++t_row_diff_base;
1297 }
1298
1299 for (; rr != row_nb_base_fun; ++rr)
1300 ++t_row_diff_base;
1301
1302 // move to next integration weight
1303 ++t_w;
1304 ++t_D;
1305 ++t_cauchy_stress;
1306 ++t_strain;
1307 ++t_eshelby_stress;
1308 ++t_h;
1309 ++t_invH;
1310 ++t_F;
1311 }
1312
1314}
1315
1316template <int S>
1317HookeElement::OpAleLhsPre_dX_dx<S>::OpAleLhsPre_dX_dx(
1318 const std::string row_field, const std::string col_field,
1319 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1320 : VolUserDataOperator(row_field, col_field, OPROW, true),
1321 dataAtPts(data_at_pts) {
1322 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1323}
1324
1325template <int S>
1326MoFEMErrorCode HookeElement::OpAleLhsPre_dX_dx<S>::doWork(int row_side,
1327 EntityType row_type,
1328 EntData &row_data) {
1330
1331 const int nb_integration_pts = row_data.getN().size1();
1332
1333 auto get_eshelby_stress_dx = [this, nb_integration_pts]() {
1335 t_eshelby_stress_dx;
1336 dataAtPts->eshelbyStress_dx->resize(81, nb_integration_pts, false);
1337 int mm = 0;
1338 for (int ii = 0; ii != 3; ++ii)
1339 for (int jj = 0; jj != 3; ++jj)
1340 for (int kk = 0; kk != 3; ++kk)
1341 for (int ll = 0; ll != 3; ++ll)
1342 t_eshelby_stress_dx.ptr(ii, jj, kk, ll) =
1343 &(*dataAtPts->eshelbyStress_dx)(mm++, 0);
1344 return t_eshelby_stress_dx;
1345 };
1346
1347 auto t_eshelby_stress_dx = get_eshelby_stress_dx();
1348
1355
1356 // Elastic stiffness tensor (4th rank tensor with minor and major
1357 // symmetry)
1359 MAT_TO_DDG(dataAtPts->stiffnessMat));
1360 auto t_cauchy_stress =
1361 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1362 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1363 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1364
1365 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1366
1367 t_eshelby_stress_dx(i, j, m, n) =
1368 (t_F(k, i) * t_D(k, j, m, l)) * t_invH(n, l);
1369 for (int ii = 0; ii != 3; ++ii)
1370 for (int jj = 0; jj != 3; ++jj)
1371 for (int mm = 0; mm != 3; ++mm)
1372 for (int nn = 0; nn != 3; ++nn) {
1373 auto &v = t_eshelby_stress_dx(ii, jj, mm, nn);
1374 v += t_invH(nn, ii) * t_cauchy_stress(mm, jj);
1375 }
1376 t_eshelby_stress_dx(i, j, k, l) *= -1;
1377
1379 t_energy_dx(m, n) = t_invH(n, j) * t_cauchy_stress(m, j);
1380
1381 for (int mm = 0; mm != 3; ++mm)
1382 for (int nn = 0; nn != 3; ++nn) {
1383 auto v = t_energy_dx(mm, nn);
1384 for (int ii = 0; ii != 3; ++ii)
1385 t_eshelby_stress_dx(ii, ii, mm, nn) += v;
1386 }
1387
1388 ++t_D;
1389 ++t_invH;
1390 ++t_cauchy_stress;
1391 ++t_eshelby_stress_dx;
1392 ++t_F;
1393 }
1394
1396}
1397
1398template <class ELEMENT>
1399HookeElement::OpPostProcHookeElement<ELEMENT>::OpPostProcHookeElement(
1400 const string row_field, boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1401 map<int, BlockData> &block_sets_ptr, moab::Interface &post_proc_mesh,
1402 std::vector<EntityHandle> &map_gauss_pts, bool is_ale, bool is_field_disp)
1403 : ELEMENT::UserDataOperator(row_field, UserDataOperator::OPROW),
1404 dataAtPts(data_at_pts), blockSetsPtr(block_sets_ptr),
1405 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), isALE(is_ale),
1406 isFieldDisp(is_field_disp) {}
1407
1408template <class ELEMENT>
1409MoFEMErrorCode HookeElement::OpPostProcHookeElement<ELEMENT>::doWork(
1410 int side, EntityType type, EntitiesFieldData::EntData &data) {
1412
1413 if (type != MBVERTEX) {
1415 }
1416
1417 auto tensor_to_tensor = [](const auto &t1, auto &t2) {
1418 t2(0, 0) = t1(0, 0);
1419 t2(1, 1) = t1(1, 1);
1420 t2(2, 2) = t1(2, 2);
1421 t2(0, 1) = t2(1, 0) = t1(1, 0);
1422 t2(0, 2) = t2(2, 0) = t1(2, 0);
1423 t2(1, 2) = t2(2, 1) = t1(2, 1);
1424 };
1425
1426 std::array<double, 9> def_val;
1427 def_val.fill(0);
1428
1429 auto make_tag = [&](auto name, auto size) {
1430 Tag th;
1431 CHKERR postProcMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
1432 MB_TAG_CREAT | MB_TAG_SPARSE,
1433 def_val.data());
1434 return th;
1435 };
1436
1437 auto th_stress = make_tag("STRESS", 9);
1438 auto th_psi = make_tag("ENERGY", 1);
1439
1440 const int nb_integration_pts = mapGaussPts.size();
1441
1446
1447 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1448 auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
1449
1450 dataAtPts->stiffnessMat->resize(36, 1, false);
1452 MAT_TO_DDG(dataAtPts->stiffnessMat));
1453 for (auto &m : (blockSetsPtr)) {
1454 const double young = m.second.E;
1455 const double poisson = m.second.PoissonRatio;
1456
1457 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
1458
1459 t_D(i, j, k, l) = 0.;
1460 t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
1461 t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
1462 0.5 * (1 - 2 * poisson);
1463 t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
1464 t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
1465 t_D(i, j, k, l) *= coefficient;
1466
1467 break; // FIXME: calculates only first block
1468 }
1469
1470 double detH = 0.;
1474 FTensor::Tensor2<double, 3, 3> t_small_strain;
1476 FTensor::Tensor2_symmetric<double, 3> t_small_strain_symm;
1477
1478 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1479
1480 if (isFieldDisp) {
1481 t_h(0, 0) += 1;
1482 t_h(1, 1) += 1;
1483 t_h(2, 2) += 1;
1484 }
1485
1486 if (!isALE) {
1487 t_small_strain_symm(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
1488 } else {
1489 CHKERR determinantTensor3by3(t_H, detH);
1490 CHKERR invertTensor3by3(t_H, detH, t_invH);
1491 t_F(i, j) = t_h(i, k) * t_invH(k, j);
1492 t_small_strain_symm(i, j) = (t_F(i, j) || t_F(j, i)) / 2.;
1493 ++t_H;
1494 }
1495
1496 t_small_strain_symm(0, 0) -= 1;
1497 t_small_strain_symm(1, 1) -= 1;
1498 t_small_strain_symm(2, 2) -= 1;
1499
1500 // symmetric tensors need improvement
1501 t_stress_symm(i, j) = t_D(i, j, k, l) * t_small_strain_symm(k, l);
1502 tensor_to_tensor(t_stress_symm, t_stress);
1503
1504 const double psi = 0.5 * t_stress_symm(i, j) * t_small_strain_symm(i, j);
1505
1506 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
1507 CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
1508 &t_stress(0, 0));
1509
1510 ++t_h;
1511 }
1512
1514}
1515
1516template <int S>
1517HookeElement::OpAnalyticalInternalStrain_dx<S>::OpAnalyticalInternalStrain_dx(
1518 const std::string row_field,
1519 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1520 StrainFunction strain_fun)
1521 : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1522 strainFun(strain_fun) {}
1523
1524template <int S>
1526HookeElement::OpAnalyticalInternalStrain_dx<S>::iNtegrate(EntData &row_data) {
1532
1533 auto get_tensor1 = [](VectorDouble &v, const int r) {
1535 &v(r + 0), &v(r + 1), &v(r + 2));
1536 };
1537
1538 const int nb_integration_pts = getGaussPts().size2();
1539 auto t_coords = getFTensor1CoordsAtGaussPts();
1540
1541 // get element volume
1542 double vol = getVolume();
1543 auto t_w = getFTensor0IntegrationWeight();
1544
1545 nF.resize(nbRows, false);
1546 nF.clear();
1547
1548 // elastic stiffness tensor (4th rank tensor with minor and major
1549 // symmetry)
1551 MAT_TO_DDG(dataAtPts->stiffnessMat));
1552
1553 // get derivatives of base functions on rows
1554 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1555 const int row_nb_base_fun = row_data.getN().size2();
1556
1557 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1558
1559 auto t_fun_strain = strainFun(t_coords);
1561 t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1562
1563 // calculate scalar weight times element volume
1564 double a = t_w * vol;
1565
1566 auto t_nf = get_tensor1(nF, 0);
1567
1568 int rr = 0;
1569 for (; rr != nbRows / 3; ++rr) {
1570 t_nf(i) += a * t_row_diff_base(j) * t_stress(i, j);
1571 ++t_row_diff_base;
1572 ++t_nf;
1573 }
1574
1575 for (; rr != row_nb_base_fun; ++rr)
1576 ++t_row_diff_base;
1577
1578 ++t_w;
1579 ++t_coords;
1580 ++t_D;
1581 }
1582
1584}
1585
1586template <int S>
1587HookeElement::OpAnalyticalInternalAleStrain_dX<S>::
1588 OpAnalyticalInternalAleStrain_dX(
1589 const std::string row_field,
1590 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1591 StrainFunction strain_fun,
1592 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1593 : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1594 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1595
1596template <int S>
1597MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dX<S>::iNtegrate(
1598 EntData &row_data) {
1604
1605 auto get_tensor1 = [](VectorDouble &v, const int r) {
1607 &v(r + 0), &v(r + 1), &v(r + 2));
1608 };
1609
1610 const int nb_integration_pts = getGaussPts().size2();
1611
1612 auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1613 auto t_coords = get_coords();
1614
1615 // get element volume
1616 double vol = getVolume();
1617 auto t_w = getFTensor0IntegrationWeight();
1618
1619 nF.resize(nbRows, false);
1620 nF.clear();
1621
1622 // elastic stiffness tensor (4th rank tensor with minor and major
1623 // symmetry)
1625 MAT_TO_DDG(dataAtPts->stiffnessMat));
1626 auto t_F = getFTensor2FromMat<3, 3>(*(dataAtPts->FMat));
1627 auto &det_H = *dataAtPts->detHVec;
1628 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1629
1630 // get derivatives of base functions on rows
1631 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1632 const int row_nb_base_fun = row_data.getN().size2();
1633
1634 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1635
1636 auto t_fun_strain = strainFun(t_coords);
1638 t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1639 FTensor::Tensor2<double, 3, 3> t_eshelby_stress;
1640 t_eshelby_stress(i, j) = -t_F(k, i) * t_stress(k, j);
1641
1642 // calculate scalar weight times element volume
1643 double a = t_w * vol * det_H[gg];
1644
1645 auto t_nf = get_tensor1(nF, 0);
1646
1647 int rr = 0;
1648 for (; rr != nbRows / 3; ++rr) {
1649 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1650 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1651 t_nf(i) += a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1652 ++t_row_diff_base;
1653 ++t_nf;
1654 }
1655
1656 for (; rr != row_nb_base_fun; ++rr)
1657 ++t_row_diff_base;
1658
1659 ++t_w;
1660 ++t_coords;
1661 ++t_F;
1662 ++t_invH;
1663 ++t_D;
1664 }
1665
1667}
1668
1669template <int S>
1670HookeElement::OpAnalyticalInternalAleStrain_dx<S>::
1671 OpAnalyticalInternalAleStrain_dx(
1672 const std::string row_field,
1673 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1674 StrainFunction strain_fun,
1675 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1676 : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1677 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1678
1679template <int S>
1680MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dx<S>::iNtegrate(
1681 EntData &row_data) {
1687
1688 auto get_tensor1 = [](VectorDouble &v, const int r) {
1690 &v(r + 0), &v(r + 1), &v(r + 2));
1691 };
1692
1693 const int nb_integration_pts = getGaussPts().size2();
1694
1695 auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1696 auto t_coords = get_coords();
1697
1698 // get element volume
1699 double vol = getVolume();
1700 auto t_w = getFTensor0IntegrationWeight();
1701
1702 nF.resize(nbRows, false);
1703 nF.clear();
1704
1705 // elastic stiffness tensor (4th rank tensor with minor and major
1706 // symmetry)
1708 MAT_TO_DDG(dataAtPts->stiffnessMat));
1709 auto &det_H = *dataAtPts->detHVec;
1710 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1711
1712 // get derivatives of base functions on rows
1713 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1714 const int row_nb_base_fun = row_data.getN().size2();
1715
1716 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1717
1718 auto t_fun_strain = strainFun(t_coords);
1720 t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1721
1722 // calculate scalar weight times element volume
1723 double a = t_w * vol * det_H[gg];
1724
1725 auto t_nf = get_tensor1(nF, 0);
1726
1727 int rr = 0;
1728 for (; rr != nbRows / 3; ++rr) {
1729 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1730 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1731 t_nf(i) += a * t_row_diff_base_pulled(j) * t_stress(i, j);
1732 ++t_row_diff_base;
1733 ++t_nf;
1734 }
1735
1736 for (; rr != row_nb_base_fun; ++rr)
1737 ++t_row_diff_base;
1738
1739 ++t_w;
1740 ++t_coords;
1741 ++t_invH;
1742 ++t_D;
1743 }
1744
1746}
1747
1748#endif // __HOOKE_ELEMENT_HPP
static Index< 'n', 3 > n
static Index< 'K', 3 > K
EntitiesFieldData::EntData EntData
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()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
struct ConvectiveMassElement BlockData
FTensor::Index< 'i', SPACE_DIM > i
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:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
UBlasVector< int > VectorInt
Definition: Types.hpp:78
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:1218
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1207
double w(const double g, const double t)
const double D
diffusivity
const double r
rate factor
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
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.
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)