153
155 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
156
158
159 protected:
160 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
161 };
162
164
166 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
168
170
171 protected:
172 boost::shared_ptr<DataAtIntegrationPts>
dataAtPts;
174 };
175
177
179 const std::string row_field, const std::string col_field,
180 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
181
183
184 protected:
185 boost::shared_ptr<DataAtIntegrationPts>
dataAtPts;
186 };
187
188 template <int S = 0>
190
192 const std::string row_field, const std::string col_field,
193 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
194 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
195
197
198 protected:
199 boost::shared_ptr<map<int, BlockData>>
201
202 boost::shared_ptr<DataAtIntegrationPts>
dataAtPts;
203 };
204
205
206
207
208
209
210
212 : public VolumeElementForcesAndSourcesCore::UserDataOperator {
213
218
219 boost::shared_ptr<DataAtIntegrationPts>
commonData;
220
222 const std::string col_field,
BlockData &data,
224 boost::shared_ptr<DataAtIntegrationPts> &common_data,
225 bool symm = true)
227 row_field, col_field, OPROWCOL, symm),
229
230 PetscErrorCode
doWork(
int row_side,
int col_side, EntityType row_type,
231 EntityType col_type,
235
238 &
m(3 * r + 0, 3 *
c + 0), &
m(3 * r + 0, 3 *
c + 1),
239 &
m(3 * r + 0, 3 *
c + 2), &
m(3 * r + 1, 3 *
c + 0),
240 &
m(3 * r + 1, 3 *
c + 1), &
m(3 * r + 1, 3 *
c + 2),
241 &
m(3 * r + 2, 3 *
c + 0), &
m(3 * r + 2, 3 *
c + 1),
242 &
m(3 * r + 2, 3 *
c + 2));
243 };
244
245 const int row_nb_dofs = row_data.
getIndices().size();
246 if (!row_nb_dofs)
248 const int col_nb_dofs = col_data.
getIndices().size();
249 if (!col_nb_dofs)
251 if (
dAta.tEts.find(getFEEntityHandle()) ==
dAta.tEts.end()) {
253 }
256 }
257
258 const bool diagonal_block =
259 (row_type == col_type) && (row_side == col_side);
260
261
262 locK.resize(row_nb_dofs, col_nb_dofs,
false);
264
265 const int row_nb_gauss_pts = row_data.
getN().size1();
266 const int row_nb_base_functions = row_data.
getN().size2();
267
272
274
275
276 auto t_w = getFTensor0IntegrationWeight();
277
278
279 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
280
282
283
284 double w = getVolume() * t_w;
285
286 for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
288 for (int col_bb = 0; col_bb != col_nb_dofs / 3; col_bb++) {
289 auto t_assemble = get_tensor2(
locK, row_bb, col_bb);
290 t_assemble(
i,
j) += density * t_row_base_func * t_col_base_func *
w;
291
292 ++t_col_base_func;
293 }
294
295 ++t_row_base_func;
296 }
297
298 ++t_w;
299 }
300
302 ADD_VALUES);
303
304
305 if (row_type != col_type || row_side != col_side) {
306 translocK.resize(col_nb_dofs, row_nb_dofs,
false);
308
310 ADD_VALUES);
311 }
312
314 }
315 };
316
318 protected:
319 boost::shared_ptr<map<int, BlockData>>
321
322 boost::shared_ptr<DataAtIntegrationPts>
dataAtPts;
323
327
328 public:
330 const std::string row_field, const std::string col_field,
331 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
332 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
333 boost::shared_ptr<VectorDouble> rho_at_gauss_pts, const double rho_n,
334 const double rho_0);
335
337 };
338
340
341 OpAssemble(
const std::string row_field,
const std::string col_field,
342 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
343 const char type, bool symm = false);
344
345
346
347
348
349
350
351
352
353
354
356 EntityType col_type,
EntData &row_data,
358
360
361 protected:
362
366
367 boost::shared_ptr<DataAtIntegrationPts>
dataAtPts;
368
371
376
378
380
381
382
383
384
385
386
387
389
390
391
392
393
394
395
396
398 };
399
401
402 OpRhs_dx(
const std::string row_field,
const std::string col_field,
403 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
404
405 protected:
407 };
408
410
411 OpLhs_dx_dx(
const std::string row_field,
const std::string col_field,
412 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
413
414 protected:
415
416
417
418
419
420
422 };
423
425
426 OpAleRhs_dx(
const std::string row_field,
const std::string col_field,
427 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
428
429 protected:
431 };
432
434
435 OpAleLhs_dx_dx(
const std::string row_field,
const std::string col_field,
436 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
437
438 protected:
439
440
441
442
443
444
446 };
447
449
450 OpAleLhs_dx_dX(
const std::string row_field,
const std::string col_field,
451 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
452
453 protected:
454
455
456
457
458
459
461 };
462
464
469
471 const std::string row_field, const std::string col_field,
472 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
473 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
474 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
475 const double rho_n, const double rho_0);
476
477 protected:
478
479
480
481
482
483
485 };
486
488
493
495 const std::string row_field, const std::string col_field,
496 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
497 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
498 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
499 const double rho_n, const double rho_0);
500
501 protected:
502
503
504
505
506
507
509 };
510
512
513 OpAleRhs_dX(
const std::string row_field,
const std::string col_field,
514 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
515
516 protected:
518 };
519
521
522 OpAleLhs_dX_dX(
const std::string row_field,
const std::string col_field,
523 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
524
525 protected:
526
527
528
529
530
531
533 };
534
536
538 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
539
541
542 private:
543 boost::shared_ptr<DataAtIntegrationPts>
dataAtPts;
544 };
545
547
548 OpAleLhs_dX_dx(
const std::string row_field,
const std::string col_field,
549 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
551
552 protected:
553
554
555
556
557
558
560 };
561
563
564 typedef boost::function<
565
567
569
570 )
571
572 >
574
576 const std::string row_field,
577 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
579
580 protected:
583 };
584
586
587 typedef boost::function<
588
590
592
593 )
594
595 >
597
599 const std::string row_field,
600 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
602 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
603
604 protected:
608 };
609
611
612 typedef boost::function<
613
615
617
618 )
619
620 >
622
624 const std::string row_field,
625 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
627 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
628
629 protected:
633 };
634
635 template <class ELEMENT>
637 boost::shared_ptr<DataAtIntegrationPts>
dataAtPts;
638 map<int, BlockData>
644
646 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
647 map<int, BlockData> &block_sets_ptr,
648 moab::Interface &post_proc_mesh,
649 std::vector<EntityHandle> &map_gauss_pts,
650 bool is_ale = false, bool is_field_disp = true);
651
654 };
655
658 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr);
659
662 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
663 const std::string element_name, const std::string x_field,
664 const std::string X_field, const bool ale);
665
667 setOperators(boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr,
668 boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr,
669 boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
670 const std::string x_field, const std::string X_field,
671 const bool ale, const bool field_disp,
672 const EntityType type = MBTET,
673 boost::shared_ptr<DataAtIntegrationPts> data_at_pts = nullptr);
674
676 calculateEnergy(DM dm, boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
677 const std::string x_field, const std::string X_field,
678 const bool ale, const bool field_disp,
680
681private:
683};
684
685template <bool D>
686HookeElement::OpCalculateStrain<D>::OpCalculateStrain(
687 const std::string row_field, const std::string col_field,
688 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
690 dataAtPts(data_at_pts) {
691 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
692}
693
694template <bool D>
695MoFEMErrorCode HookeElement::OpCalculateStrain<D>::doWork(
int row_side,
696 EntityType row_type,
701
702 const int nb_integration_pts = getGaussPts().size2();
703 dataAtPts->smallStrainMat->resize(6, nb_integration_pts, false);
704 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
705 auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hMat));
706
707 for (int gg = 0; gg != nb_integration_pts; ++gg) {
708 t_strain(
i,
j) = (t_h(
i,
j) || t_h(
j,
i)) / 2.;
709
710
712 t_strain(0, 0) -= 1;
713 t_strain(1, 1) -= 1;
714 t_strain(2, 2) -= 1;
715 }
716
717 ++t_strain;
718 ++t_h;
719 }
721}
722
723template <int S>
724HookeElement::OpAleLhs_dx_dx<S>::OpAleLhs_dx_dx(
725 const std::string row_field, const std::string col_field,
726 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
727 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
728
729template <int S>
733
734
735
738 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
739 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
741 };
742
747
748
749 double vol = getVolume();
750
751
752 auto t_w = getFTensor0IntegrationWeight();
753
754
756 const int row_nb_base_fun = row_data.
getN().size2();
757
758
759
762
763 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
764 auto &det_H = *dataAtPts->detHVec;
765
766
767 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
768
769
770 double a = t_w * vol * det_H[gg];
771
772
773 int rr = 0;
774 for (; rr != nbRows / 3; ++rr) {
775
776
777 auto t_m = get_tensor2(
K, 3 * rr, 0);
778
780 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
781
783
784
785
786 t_rowD(
l,
j,
k) = t_D(
i,
j,
k,
l) * (
a * t_row_diff_base_pulled(
i));
787
788
790
791
792 for (int cc = 0; cc != nbCols / 3; ++cc) {
793
795 t_col_diff_base_pulled(
j) = t_col_diff_base(
i) * t_invH(
i,
j);
796
797
798 t_m(
i,
j) += t_rowD(
i,
j,
k) * t_col_diff_base_pulled(
k);
799
800
801 ++t_col_diff_base;
802
803
804 ++t_m;
805 }
806
807
808 ++t_row_diff_base;
809 }
810
811 for (; rr != row_nb_base_fun; ++rr)
812 ++t_row_diff_base;
813
814
815 ++t_w;
816 ++t_D;
817 ++t_invH;
818 }
819
821}
822
823template <int S>
824HookeElement::OpCalculateHomogeneousStiffness<S>::
825 OpCalculateHomogeneousStiffness(
826 const std::string row_field, const std::string col_field,
827 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
828 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
830 blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts) {
831 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
832}
833
834template <int S>
835MoFEMErrorCode HookeElement::OpCalculateHomogeneousStiffness<S>::doWork(
836 int row_side, EntityType row_type,
EntData &row_data) {
838
839 for (
auto &
m : (*blockSetsPtr)) {
840 if (
m.second.tEts.find(getFEEntityHandle()) !=
m.second.tEts.end()) {
841
842 dataAtPts->stiffnessMat->resize(36, 1, false);
845 const double young =
m.second.E;
846 const double poisson =
m.second.PoissonRatio;
847
848
849 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
850
855
856 t_D(
i,
j,
k,
l) = 0.;
857
858 t_D(0, 0, 0, 0) = 1 - poisson;
859 t_D(1, 1, 1, 1) = 1 - poisson;
860 t_D(2, 2, 2, 2) = 1 - poisson;
861
862 t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * poisson);
863 t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * poisson);
864 t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * poisson);
865
866 t_D(0, 0, 1, 1) = poisson;
867 t_D(1, 1, 0, 0) = poisson;
868 t_D(0, 0, 2, 2) = poisson;
869 t_D(2, 2, 0, 0) = poisson;
870 t_D(1, 1, 2, 2) = poisson;
871 t_D(2, 2, 1, 1) = poisson;
872
873 t_D(
i,
j,
k,
l) *= coefficient;
874
875 break;
876 }
877 }
878
880}
881
882template <int S>
883HookeElement::OpCalculateStress<S>::OpCalculateStress(
884 const std::string row_field, const std::string col_field,
885 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
887 dataAtPts(data_at_pts) {
888 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
889}
890
891template <int S>
892MoFEMErrorCode HookeElement::OpCalculateStress<S>::doWork(
int row_side,
893 EntityType row_type,
896
897 const int nb_integration_pts = getGaussPts().size2();
898 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
899 dataAtPts->cauchyStressMat->resize(6, nb_integration_pts, false);
900 auto t_cauchy_stress =
901 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
902
907
908
909
912 for (int gg = 0; gg != nb_integration_pts; ++gg) {
913 t_cauchy_stress(
i,
j) = t_D(
i,
j,
k,
l) * t_strain(
k,
l);
914 ++t_strain;
915 ++t_cauchy_stress;
916 ++t_D;
917 }
919}
920
921template <int S>
922HookeElement::OpLhs_dx_dx<S>::OpLhs_dx_dx(
923 const std::string row_field, const std::string col_field,
924 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
925 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
926
927template <int S>
931
932
933
936 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
937 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
939 };
940
945
946
947 double vol = getVolume();
948
949
950 auto t_w = getFTensor0IntegrationWeight();
951
952
954 const int row_nb_base_fun = row_data.
getN().size2();
955
956
957
960
961
962 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
963
964
965 double a = t_w * vol;
966
967
968 int rr = 0;
969 for (; rr != nbRows / 3; ++rr) {
970
971
972 auto t_m = get_tensor2(
K, 3 * rr, 0);
973
974
976
978
979
980
981 t_rowD(
l,
j,
k) = t_D(
i,
j,
k,
l) * (
a * t_row_diff_base(
i));
982
983
984 for (int cc = 0; cc != nbCols / 3; ++cc) {
985
986
987 t_m(
i,
j) += t_rowD(
i,
j,
k) * t_col_diff_base(
k);
988
989
990 ++t_col_diff_base;
991
992
993 ++t_m;
994 }
995
996
997 ++t_row_diff_base;
998 }
999
1000 for (; rr != row_nb_base_fun; ++rr)
1001 ++t_row_diff_base;
1002
1003
1004 ++t_w;
1005 ++t_D;
1006 }
1007
1009}
1010
1011template <int S>
1012HookeElement::OpAleLhs_dx_dX<S>::OpAleLhs_dx_dX(
1013 const std::string row_field, const std::string col_field,
1014 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1015 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
1016
1017template <int S>
1021
1022
1023
1026 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
1027 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
1029 };
1030
1037
1038
1039 double vol = getVolume();
1040
1041
1042 auto t_w = getFTensor0IntegrationWeight();
1043
1044
1046 const int row_nb_base_fun = row_data.
getN().size2();
1047
1048
1049
1052
1053 auto t_cauchy_stress =
1054 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1055 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1056 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1057 auto &det_H = *dataAtPts->detHVec;
1058
1059
1060 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1061
1062
1063 double a = t_w * vol * det_H[gg];
1064
1066 t_F_dX(
i,
j,
k,
l) = -(t_h(
i,
m) * t_invH(
m,
k)) * t_invH(
l,
j);
1067
1068
1069 int rr = 0;
1070 for (; rr != nbRows / 3; ++rr) {
1071
1072
1073 auto t_m = get_tensor2(
K, 3 * rr, 0);
1074
1076 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1077
1079 t_row_stress(
i) =
a * t_row_diff_base_pulled(
j) * t_cauchy_stress(
i,
j);
1080
1082 t_row_diff_base_pulled_dX(
j,
k,
l) =
1083 -(t_invH(
i,
k) * t_row_diff_base(
i)) * t_invH(
l,
j);
1084
1086 t_row_dX_stress(
i,
k,
l) =
1087 a * (t_row_diff_base_pulled_dX(
j,
k,
l) * t_cauchy_stress(
j,
i));
1088
1090 t_row_D(
l,
j,
k) = (
a * t_row_diff_base_pulled(
i)) * t_D(
i,
j,
k,
l);
1091
1093
1094
1095 t_row_stress_dX(
i,
j,
k) = 0;
1096 for (int ii = 0; ii != 3; ++ii)
1097 for (int mm = 0; mm != 3; ++mm)
1098 for (int nn = 0; nn != 3; ++nn) {
1099 auto &
v = t_row_stress_dX(ii, mm, nn);
1100 for (int kk = 0; kk != 3; ++kk)
1101 for (int ll = 0; ll != 3; ++ll)
1102 v += t_row_D(ii, kk, ll) * t_F_dX(kk, ll, mm, nn);
1103 }
1104
1105
1107
1108
1109 for (int cc = 0; cc != nbCols / 3; ++cc) {
1110
1111 t_m(
i,
k) += t_row_stress(
i) * (t_invH(
j,
k) * t_col_diff_base(
j));
1112 t_m(
i,
k) += t_row_dX_stress(
i,
k,
l) * t_col_diff_base(
l);
1113 t_m(
i,
k) += t_row_stress_dX(
i,
k,
l) * t_col_diff_base(
l);
1114
1115
1116 ++t_col_diff_base;
1117
1118
1119 ++t_m;
1120 }
1121
1122
1123 ++t_row_diff_base;
1124 }
1125
1126 for (; rr != row_nb_base_fun; ++rr)
1127 ++t_row_diff_base;
1128
1129
1130 ++t_w;
1131 ++t_D;
1132 ++t_cauchy_stress;
1133 ++t_invH;
1134 ++t_h;
1135 }
1136
1138}
1139
1140template <int S>
1141HookeElement::OpAleLhs_dX_dX<S>::OpAleLhs_dX_dX(
1142 const std::string row_field, const std::string col_field,
1143 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1144 :
OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
1145
1146template <int S>
1150
1151
1152
1155 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
1156 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
1158 };
1159
1166
1167
1168 double vol = getVolume();
1169
1170
1171 auto t_w = getFTensor0IntegrationWeight();
1172
1173
1175 const int row_nb_base_fun = row_data.
getN().size2();
1176
1177
1178
1181 auto t_cauchy_stress =
1182 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1183 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
1184 auto t_eshelby_stress =
1185 getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
1186 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1187 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1188 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1189 auto &det_H = *dataAtPts->detHVec;
1190
1191
1192 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1193
1194
1195 double a = t_w * vol * det_H[gg];
1196
1198 t_F_dX(
i,
j,
k,
l) = -(t_h(
i,
m) * t_invH(
m,
k)) * t_invH(
l,
j);
1199
1201 t_D_strain_dX(
i,
j,
m,
n) = 0.;
1202 for (int ii = 0; ii != 3; ++ii)
1203 for (int jj = 0; jj != 3; ++jj)
1204 for (int ll = 0; ll != 3; ++ll)
1205 for (int kk = 0; kk != 3; ++kk) {
1206 auto &
v = t_D_strain_dX(ii, jj, kk, ll);
1207 for (int mm = 0; mm != 3; ++mm)
1208 for (int nn = 0; nn != 3; ++nn)
1209 v += t_D(ii, jj, mm, nn) * t_F_dX(mm, nn, kk, ll);
1210 }
1211
1213 t_eshelby_stress_dX(
i,
j,
m,
n) = t_F(
k,
i) * t_D_strain_dX(
k,
j,
m,
n);
1214
1215 for (int ii = 0; ii != 3; ++ii)
1216 for (int jj = 0; jj != 3; ++jj)
1217 for (int mm = 0; mm != 3; ++mm)
1218 for (int nn = 0; nn != 3; ++nn) {
1219 auto &
v = t_eshelby_stress_dX(ii, jj, mm, nn);
1220 for (int kk = 0; kk != 3; ++kk)
1221 v += t_F_dX(kk, ii, mm, nn) * t_cauchy_stress(kk, jj);
1222 }
1223
1224 t_eshelby_stress_dX(
i,
j,
k,
l) *= -1;
1225
1227 t_energy_dX(
k,
l) = t_F_dX(
i,
j,
k,
l) * t_cauchy_stress(
i,
j);
1228 t_energy_dX(
k,
l) +=
1229 (t_strain(
m,
n) * t_D(
m,
n,
i,
j)) * t_F_dX(
i,
j,
k,
l);
1230 t_energy_dX(
k,
l) /= 2.;
1231
1232 for (int kk = 0; kk != 3; ++kk)
1233 for (int ll = 0; ll != 3; ++ll) {
1234 auto v = t_energy_dX(kk, ll);
1235 for (int ii = 0; ii != 3; ++ii)
1236 t_eshelby_stress_dX(ii, ii, kk, ll) +=
v;
1237 }
1238
1239
1240 int rr = 0;
1241 for (; rr != nbRows / 3; ++rr) {
1242
1243
1244 auto t_m = get_tensor2(
K, 3 * rr, 0);
1245
1247 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1248
1250 t_row_stress(
i) =
a * t_row_diff_base_pulled(
j) * t_eshelby_stress(
i,
j);
1251
1253 t_row_diff_base_pulled_dX(
j,
k,
l) =
1254 -(t_row_diff_base(
i) * t_invH(
i,
k)) * t_invH(
l,
j);
1255
1257 t_row_dX_stress(
i,
k,
l) =
1258 a * (t_row_diff_base_pulled_dX(
j,
k,
l) * t_eshelby_stress(
i,
j));
1259
1261 t_row_stress_dX(
i,
m,
n) =
1262 a * t_row_diff_base_pulled(
j) * t_eshelby_stress_dX(
i,
j,
m,
n);
1263
1264
1266
1267
1268 for (int cc = 0; cc != nbCols / 3; ++cc) {
1269
1270 t_m(
i,
k) += t_row_stress(
i) * (t_invH(
j,
k) * t_col_diff_base(
j));
1271 t_m(
i,
k) += t_row_dX_stress(
i,
k,
l) * t_col_diff_base(
l);
1272 t_m(
i,
k) += t_row_stress_dX(
i,
k,
l) * t_col_diff_base(
l);
1273
1274
1275 ++t_col_diff_base;
1276
1277
1278 ++t_m;
1279 }
1280
1281
1282 ++t_row_diff_base;
1283 }
1284
1285 for (; rr != row_nb_base_fun; ++rr)
1286 ++t_row_diff_base;
1287
1288
1289 ++t_w;
1290 ++t_D;
1291 ++t_cauchy_stress;
1292 ++t_strain;
1293 ++t_eshelby_stress;
1294 ++t_h;
1295 ++t_invH;
1296 ++t_F;
1297 }
1298
1300}
1301
1302template <int S>
1303HookeElement::OpAleLhsPre_dX_dx<S>::OpAleLhsPre_dX_dx(
1304 const std::string row_field, const std::string col_field,
1305 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1307 dataAtPts(data_at_pts) {
1308 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1309}
1310
1311template <int S>
1312MoFEMErrorCode HookeElement::OpAleLhsPre_dX_dx<S>::doWork(
int row_side,
1313 EntityType row_type,
1316
1317 const int nb_integration_pts = row_data.
getN().size1();
1318
1319 auto get_eshelby_stress_dx = [this, nb_integration_pts]() {
1321 t_eshelby_stress_dx;
1322 dataAtPts->eshelbyStress_dx->resize(81, nb_integration_pts, false);
1323 int mm = 0;
1324 for (int ii = 0; ii != 3; ++ii)
1325 for (int jj = 0; jj != 3; ++jj)
1326 for (int kk = 0; kk != 3; ++kk)
1327 for (int ll = 0; ll != 3; ++ll)
1328 t_eshelby_stress_dx.ptr(ii, jj, kk, ll) =
1329 &(*dataAtPts->eshelbyStress_dx)(mm++, 0);
1330 return t_eshelby_stress_dx;
1331 };
1332
1333 auto t_eshelby_stress_dx = get_eshelby_stress_dx();
1334
1341
1342
1343
1346 auto t_cauchy_stress =
1347 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1348 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1349 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1350
1351 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1352
1353 t_eshelby_stress_dx(
i,
j,
m,
n) =
1354 (t_F(
k,
i) * t_D(
k,
j,
m,
l)) * t_invH(
n,
l);
1355 for (int ii = 0; ii != 3; ++ii)
1356 for (int jj = 0; jj != 3; ++jj)
1357 for (int mm = 0; mm != 3; ++mm)
1358 for (int nn = 0; nn != 3; ++nn) {
1359 auto &
v = t_eshelby_stress_dx(ii, jj, mm, nn);
1360 v += t_invH(nn, ii) * t_cauchy_stress(mm, jj);
1361 }
1362 t_eshelby_stress_dx(
i,
j,
k,
l) *= -1;
1363
1365 t_energy_dx(
m,
n) = t_invH(
n,
j) * t_cauchy_stress(
m,
j);
1366
1367 for (int mm = 0; mm != 3; ++mm)
1368 for (int nn = 0; nn != 3; ++nn) {
1369 auto v = t_energy_dx(mm, nn);
1370 for (int ii = 0; ii != 3; ++ii)
1371 t_eshelby_stress_dx(ii, ii, mm, nn) +=
v;
1372 }
1373
1374 ++t_D;
1375 ++t_invH;
1376 ++t_cauchy_stress;
1377 ++t_eshelby_stress_dx;
1378 ++t_F;
1379 }
1380
1382}
1383
1384template <class ELEMENT>
1385HookeElement::OpPostProcHookeElement<ELEMENT>::OpPostProcHookeElement(
1386 const string row_field, boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1387 map<int, BlockData> &block_sets_ptr, moab::Interface &post_proc_mesh,
1388 std::vector<EntityHandle> &map_gauss_pts, bool is_ale, bool is_field_disp)
1390 dataAtPts(data_at_pts), blockSetsPtr(block_sets_ptr),
1391 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), isALE(is_ale),
1392 isFieldDisp(is_field_disp) {}
1393
1394template <class ELEMENT>
1395MoFEMErrorCode HookeElement::OpPostProcHookeElement<ELEMENT>::doWork(
1398
1399 if (type != MBVERTEX) {
1401 }
1402
1403 auto tensor_to_tensor = [](const auto &t1, auto &t2) {
1404 t2(0, 0) = t1(0, 0);
1405 t2(1, 1) = t1(1, 1);
1406 t2(2, 2) = t1(2, 2);
1407 t2(0, 1) = t2(1, 0) = t1(1, 0);
1408 t2(0, 2) = t2(2, 0) = t1(2, 0);
1409 t2(1, 2) = t2(2, 1) = t1(2, 1);
1410 };
1411
1412 std::array<double, 9> def_val;
1413 def_val.fill(0);
1414
1415 auto make_tag = [&](auto name, auto size) {
1417 CHKERR postProcMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE,
th,
1418 MB_TAG_CREAT | MB_TAG_SPARSE,
1419 def_val.data());
1421 };
1422
1423 auto th_stress = make_tag("STRESS", 9);
1424 auto th_psi = make_tag("ENERGY", 1);
1425
1426 const int nb_integration_pts = mapGaussPts.size();
1427
1432
1433 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1434 auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
1435
1436 dataAtPts->stiffnessMat->resize(36, 1, false);
1439
1443 if (type == MBTRI || type == MBQUAD) {
1445 auto &m_field = this->getPtrFE()->mField;
1446 CHKERR m_field.get_moab().get_adjacencies(&ent, 1, 3,
false, ents,
1447 moab::Interface::UNION);
1448#ifndef NDEBUG
1449 if (ents.empty())
1451 "Could not find a 3D element adjacent to a given face element");
1452#endif
1453 ent_3d = ents.front();
1454 }
1455
1456 bool found_block = false;
1457 int block_id = -1;
1458 for (
auto &
m : (blockSetsPtr)) {
1459 if (
m.second.tEts.find(ent_3d) !=
m.second.tEts.end()) {
1460 const double young =
m.second.E;
1461 const double poisson =
m.second.PoissonRatio;
1462 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
1463 block_id =
m.second.iD;
1464
1465 t_D(
i,
j,
k,
l) = 0.;
1466 t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
1467 t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
1468 0.5 * (1 - 2 * poisson);
1469 t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
1470 t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
1471 t_D(
i,
j,
k,
l) *= coefficient;
1472
1473 found_block = true;
1474 break;
1475 }
1476 }
1477 if (!found_block)
1479 "Element not found in any of material blocksets");
1480
1481 int def_val_int = 0;
1483 CHKERR postProcMesh.tag_get_handle(
"MAT_ELASTIC", 1, MB_TYPE_INTEGER, tag_mat,
1484 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
1485 double detH = 0.;
1492
1493 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1494
1495 if (isFieldDisp) {
1496 t_h(0, 0) += 1;
1497 t_h(1, 1) += 1;
1498 t_h(2, 2) += 1;
1499 }
1500
1501 if (!isALE) {
1502 t_small_strain_symm(
i,
j) = (t_h(
i,
j) || t_h(
j,
i)) / 2.;
1503 } else {
1506 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
1507 t_small_strain_symm(
i,
j) = (t_F(
i,
j) || t_F(
j,
i)) / 2.;
1508 ++t_H;
1509 }
1510
1511 t_small_strain_symm(0, 0) -= 1;
1512 t_small_strain_symm(1, 1) -= 1;
1513 t_small_strain_symm(2, 2) -= 1;
1514
1515
1516 t_stress_symm(
i,
j) = t_D(
i,
j,
k,
l) * t_small_strain_symm(
k,
l);
1517 tensor_to_tensor(t_stress_symm, t_stress);
1518
1519 const double psi = 0.5 * t_stress_symm(
i,
j) * t_small_strain_symm(
i,
j);
1520
1521 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
1522 CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
1523 &t_stress(0, 0));
1524 CHKERR postProcMesh.tag_set_data(tag_mat, &mapGaussPts[gg], 1,
1525 &block_id);
1526
1527 ++t_h;
1528 }
1529
1531}
1532
1533template <int S>
1534HookeElement::OpAnalyticalInternalStrain_dx<S>::OpAnalyticalInternalStrain_dx(
1535 const std::string row_field,
1536 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1537 StrainFunction strain_fun)
1538 :
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1539 strainFun(strain_fun) {}
1540
1541template <int S>
1543HookeElement::OpAnalyticalInternalStrain_dx<S>::iNtegrate(
EntData &row_data) {
1549
1552 &
v(r + 0), &
v(r + 1), &
v(r + 2));
1553 };
1554
1555 const int nb_integration_pts = getGaussPts().size2();
1556 auto t_coords = getFTensor1CoordsAtGaussPts();
1557
1558
1559 double vol = getVolume();
1560 auto t_w = getFTensor0IntegrationWeight();
1561
1562 nF.resize(nbRows, false);
1563 nF.clear();
1564
1565
1566
1569
1570
1572 const int row_nb_base_fun = row_data.
getN().size2();
1573
1574 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1575
1576 auto t_fun_strain = strainFun(t_coords);
1578 t_stress(
i,
j) = -t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
1579
1580
1581 double a = t_w * vol;
1582
1583 auto t_nf = get_tensor1(nF, 0);
1584
1585 int rr = 0;
1586 for (; rr != nbRows / 3; ++rr) {
1587 t_nf(
i) +=
a * t_row_diff_base(
j) * t_stress(
i,
j);
1588 ++t_row_diff_base;
1589 ++t_nf;
1590 }
1591
1592 for (; rr != row_nb_base_fun; ++rr)
1593 ++t_row_diff_base;
1594
1595 ++t_w;
1596 ++t_coords;
1597 ++t_D;
1598 }
1599
1601}
1602
1603template <int S>
1604HookeElement::OpAnalyticalInternalAleStrain_dX<S>::
1605 OpAnalyticalInternalAleStrain_dX(
1606 const std::string row_field,
1607 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1608 StrainFunction strain_fun,
1609 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1610 :
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1611 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1612
1613template <int S>
1614MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dX<S>::iNtegrate(
1621
1624 &
v(r + 0), &
v(r + 1), &
v(r + 2));
1625 };
1626
1627 const int nb_integration_pts = getGaussPts().size2();
1628
1629 auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1630 auto t_coords = get_coords();
1631
1632
1633 double vol = getVolume();
1634 auto t_w = getFTensor0IntegrationWeight();
1635
1636 nF.resize(nbRows, false);
1637 nF.clear();
1638
1639
1640
1643 auto t_F = getFTensor2FromMat<3, 3>(*(dataAtPts->FMat));
1644 auto &det_H = *dataAtPts->detHVec;
1645 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1646
1647
1649 const int row_nb_base_fun = row_data.
getN().size2();
1650
1651 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1652
1653 auto t_fun_strain = strainFun(t_coords);
1655 t_stress(
i,
j) = -t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
1657 t_eshelby_stress(
i,
j) = -t_F(
k,
i) * t_stress(
k,
j);
1658
1659
1660 double a = t_w * vol * det_H[gg];
1661
1662 auto t_nf = get_tensor1(nF, 0);
1663
1664 int rr = 0;
1665 for (; rr != nbRows / 3; ++rr) {
1667 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1668 t_nf(
i) +=
a * t_row_diff_base_pulled(
j) * t_eshelby_stress(
i,
j);
1669 ++t_row_diff_base;
1670 ++t_nf;
1671 }
1672
1673 for (; rr != row_nb_base_fun; ++rr)
1674 ++t_row_diff_base;
1675
1676 ++t_w;
1677 ++t_coords;
1678 ++t_F;
1679 ++t_invH;
1680 ++t_D;
1681 }
1682
1684}
1685
1686template <int S>
1687HookeElement::OpAnalyticalInternalAleStrain_dx<S>::
1688 OpAnalyticalInternalAleStrain_dx(
1689 const std::string row_field,
1690 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1691 StrainFunction strain_fun,
1692 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1693 :
OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1694 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1695
1696template <int S>
1697MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dx<S>::iNtegrate(
1704
1707 &
v(r + 0), &
v(r + 1), &
v(r + 2));
1708 };
1709
1710 const int nb_integration_pts = getGaussPts().size2();
1711
1712 auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1713 auto t_coords = get_coords();
1714
1715
1716 double vol = getVolume();
1717 auto t_w = getFTensor0IntegrationWeight();
1718
1719 nF.resize(nbRows, false);
1720 nF.clear();
1721
1722
1723
1726 auto &det_H = *dataAtPts->detHVec;
1727 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1728
1729
1731 const int row_nb_base_fun = row_data.
getN().size2();
1732
1733 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1734
1735 auto t_fun_strain = strainFun(t_coords);
1737 t_stress(
i,
j) = -t_D(
i,
j,
k,
l) * t_fun_strain(
k,
l);
1738
1739
1740 double a = t_w * vol * det_H[gg];
1741
1742 auto t_nf = get_tensor1(nF, 0);
1743
1744 int rr = 0;
1745 for (; rr != nbRows / 3; ++rr) {
1747 t_row_diff_base_pulled(
i) = t_row_diff_base(
j) * t_invH(
j,
i);
1748 t_nf(
i) +=
a * t_row_diff_base_pulled(
j) * t_stress(
i,
j);
1749 ++t_row_diff_base;
1750 ++t_nf;
1751 }
1752
1753 for (; rr != row_nb_base_fun; ++rr)
1754 ++t_row_diff_base;
1755
1756 ++t_w;
1757 ++t_coords;
1758 ++t_invH;
1759 ++t_D;
1760 }
1761
1763}
1764
1765#endif
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)
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr)
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
#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
#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)
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.
UBlasVector< int > VectorInt
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.
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
Range tEts
elements in block set
double rho0
reference density
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 degrees of freedom on entity.
@ OPROWCOL
operator doWork is executed on FE rows &columns
intrusive_ptr for managing petsc objects
Volume finite element base.
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
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for spatial momentum.
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for material momentum.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for material momentum.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for spatial momentum.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate B^T D B operator.
MoFEMErrorCode iNtegrate(EntData &row_data)
MoFEMErrorCode iNtegrate(EntData &row_data)
MoFEMErrorCode iNtegrate(EntData &row_data)
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
MoFEMErrorCode iNtegrate(EntData &row_data)
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
MoFEMErrorCode iNtegrate(EntData &row_data)
int nbCols
number if dof on column
virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
bool isDiag
true if this block is on diagonal
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
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
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
data at integration 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< map< int, BlockData > > blockSetsPtr
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Assemble mass matrix for elastic element TODO: CHANGE FORMULA *.
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
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
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 B^T D B operator.
std::vector< EntityHandle > & mapGaussPts
map< int, BlockData > & blockSetsPtr
moab::Interface & postProcMesh
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode iNtegrate(EntData &row_data)