v0.13.1
Loading...
Searching...
No Matches
HookeElement.cpp
Go to the documentation of this file.
1/** \file HookeElement.cpp
2 * \example HookeElement.cpp
3 * \brief Operators and data structures for linear elastic analysis
4 *
5 * See as well header file HookeElement.hpp
6 *
7 * Implemention of operators for Hooke material. Implementation is extended to
8 * the case when the mesh is moving as results of topological changes, also the
9 * calculation of material forces and associated tangent matrices are added to
10 * implementation.
11 *
12 * In other words spatial deformation is small but topological changes large.
13
14 */
15
16
17
18#include <MoFEM.hpp>
19using namespace MoFEM;
20#include <HookeElement.hpp>
21
22HookeElement::OpCalculateStrainAle::OpCalculateStrainAle(
23 const std::string row_field, const std::string col_field,
24 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
25 : VolUserDataOperator(row_field, col_field, OPROW, false),
26 dataAtPts(data_at_pts) {
27 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
28}
29
30MoFEMErrorCode HookeElement::OpCalculateStrainAle::doWork(int row_side,
31 EntityType row_type,
32 EntData &row_data) {
37 // get number of integration points
38 const int nb_integration_pts = getGaussPts().size2();
39 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
40 auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
41
42 dataAtPts->detHVec->resize(nb_integration_pts, false);
43 dataAtPts->invHMat->resize(9, nb_integration_pts, false);
44 dataAtPts->FMat->resize(9, nb_integration_pts, false);
45 dataAtPts->smallStrainMat->resize(6, nb_integration_pts, false);
46
47 auto t_detH = getFTensor0FromVec(*dataAtPts->detHVec);
48 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
49 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
50 auto t_strain = getFTensor2SymmetricFromMat<3>(*dataAtPts->smallStrainMat);
51
52 for (int gg = 0; gg != nb_integration_pts; ++gg) {
53 CHKERR determinantTensor3by3(t_H, t_detH);
54 CHKERR invertTensor3by3(t_H, t_detH, t_invH);
55 t_F(i, j) = t_h(i, k) * t_invH(k, j);
56 t_strain(i, j) = (t_F(i, j) || t_F(j, i)) / 2.;
57
58 t_strain(0, 0) -= 1;
59 t_strain(1, 1) -= 1;
60 t_strain(2, 2) -= 1;
61
62 ++t_strain;
63 ++t_h;
64 ++t_H;
65 ++t_detH;
66 ++t_invH;
67 ++t_F;
68 }
70}
71
72HookeElement::OpCalculateEnergy::OpCalculateEnergy(
73 const std::string row_field, const std::string col_field,
74 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
75 SmartPetscObj<Vec> ghost_vec)
76 : VolUserDataOperator(row_field, col_field, OPROW, true),
77 dataAtPts(data_at_pts), ghostVec(ghost_vec, true) {
78 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
79}
80
81MoFEMErrorCode HookeElement::OpCalculateEnergy::doWork(int row_side,
82 EntityType row_type,
83 EntData &row_data) {
85
86 // get number of integration points
87 const int nb_integration_pts = getGaussPts().size2();
88 auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
89 auto t_cauchy_stress =
90 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
91 dataAtPts->energyVec->resize(nb_integration_pts, false);
93 &*(dataAtPts->energyVec->data().begin()));
94
97
98 for (int gg = 0; gg != nb_integration_pts; ++gg) {
99 t_energy = (t_strain(i, j) * t_cauchy_stress(i, j)) / 2.;
100 ++t_strain;
101 ++t_cauchy_stress;
102 ++t_energy;
103 }
104
105 if (ghostVec.get()) {
106 // get element volume
107 double vol = getVolume();
108 // get intergrayion weights
109 auto t_w = getFTensor0IntegrationWeight();
110 auto &det_H = *dataAtPts->detHVec;
112 &*(dataAtPts->energyVec->data().begin()));
113 double energy = 0;
114 for (int gg = 0; gg != nb_integration_pts; ++gg) {
115 // calculate scalar weight times element volume
116 double a = t_w * vol;
117 if (det_H.size()) {
118 a *= det_H[gg];
119 }
120 energy += a * t_energy;
121 ++t_energy;
122 ++t_w;
123 }
124 CHKERR VecSetValue(ghostVec, 0, energy, ADD_VALUES);
125 }
126
128}
129
130HookeElement::OpCalculateEshelbyStress::OpCalculateEshelbyStress(
131 const std::string row_field, const std::string col_field,
132 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
133 : VolUserDataOperator(row_field, col_field, OPROW, true),
134 dataAtPts(data_at_pts) {
135 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
136}
137
138MoFEMErrorCode HookeElement::OpCalculateEshelbyStress::doWork(
139 int row_side, EntityType row_type, EntData &row_data) {
141 // get number of integration points
142 const int nb_integration_pts = getGaussPts().size2();
144 &*(dataAtPts->energyVec->data().begin()));
145 auto t_cauchy_stress =
146 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
147 auto t_F = getFTensor2FromMat<3, 3>(*(dataAtPts->FMat));
148 dataAtPts->eshelbyStressMat->resize(9, nb_integration_pts, false);
149 auto t_eshelby_stress =
150 getFTensor2FromMat<3, 3>(*(dataAtPts->eshelbyStressMat));
151
155
156 for (int gg = 0; gg != nb_integration_pts; ++gg) {
157 t_eshelby_stress(i, j) = -t_F(k, i) * t_cauchy_stress(k, j);
158 t_eshelby_stress(0, 0) += t_energy;
159 t_eshelby_stress(1, 1) += t_energy;
160 t_eshelby_stress(2, 2) += t_energy;
161 ++t_cauchy_stress;
162 ++t_energy;
163 ++t_eshelby_stress;
164 ++t_F;
165 }
167}
168
169HookeElement::OpAssemble::OpAssemble(
170 const std::string row_field, const std::string col_field,
171 boost::shared_ptr<DataAtIntegrationPts> data_at_pts, const char type,
172 bool symm)
173 : VolUserDataOperator(row_field, col_field, type, symm),
174 dataAtPts(data_at_pts) {}
175
176MoFEMErrorCode HookeElement::OpAssemble::doWork(int row_side, int col_side,
177 EntityType row_type,
178 EntityType col_type,
179 EntData &row_data,
180 EntData &col_data) {
181
183
184 // get number of dofs on row
185 nbRows = row_data.getIndices().size();
186 // if no dofs on row, exit that work, nothing to do here
187 if (!nbRows)
189
190 // get number of dofs on column
191 nbCols = col_data.getIndices().size();
192 // if no dofs on Columbia, exit nothing to do here
193 if (!nbCols)
195
196 // K_ij matrix will have 3 times the number of degrees of freedom of the
197 // i-th entity set (nbRows)
198 // and 3 times the number of degrees of freedom of the j-th entity set
199 // (nbCols)
200 K.resize(nbRows, nbCols, false);
201 K.clear();
202
203 // get number of integration points
204 nbIntegrationPts = getGaussPts().size2();
205 // check if entity block is on matrix diagonal
206 if (row_side == col_side && row_type == col_type) {
207 isDiag = true;
208 } else {
209 isDiag = false;
210 }
211
212 // integrate local matrix for entity block
213 CHKERR iNtegrate(row_data, col_data);
214
215 // assemble local matrix
216 CHKERR aSsemble(row_data, col_data);
217
219}
220
221MoFEMErrorCode HookeElement::OpAssemble::doWork(int row_side,
222 EntityType row_type,
223 EntData &row_data) {
225
226 // get number of dofs on row
227 nbRows = row_data.getIndices().size();
228 // if no dofs on row, exit that work, nothing to do here
229 if (!nbRows)
231
232 nF.resize(nbRows, false);
233 nF.clear();
234
235 // get number of integration points
236 nbIntegrationPts = getGaussPts().size2();
237
238 // integrate local matrix for entity block
239 CHKERR iNtegrate(row_data);
240
241 // assemble local matrix
242 CHKERR aSsemble(row_data);
243
245}
246
247MoFEMErrorCode HookeElement::OpAssemble::iNtegrate(EntData &row_data,
248 EntData &col_data) {
250 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
252};
253
254MoFEMErrorCode HookeElement::OpAssemble::iNtegrate(EntData &row_data) {
256 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
258};
259
260MoFEMErrorCode HookeElement::OpAssemble::aSsemble(EntData &row_data,
261 EntData &col_data) {
263
264 // get pointer to first global index on row
265 const int *row_indices = &*row_data.getIndices().data().begin();
266 // get pointer to first global index on column
267 const int *col_indices = &*col_data.getIndices().data().begin();
268
269 auto &data = *dataAtPts;
270 if (!data.forcesOnlyOnEntitiesRow.empty()) {
271 rowIndices.resize(nbRows, false);
272 noalias(rowIndices) = row_data.getIndices();
273 row_indices = &rowIndices[0];
274 VectorDofs &dofs = row_data.getFieldDofs();
275 VectorDofs::iterator dit = dofs.begin();
276 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
277 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
278 data.forcesOnlyOnEntitiesRow.end()) {
279 rowIndices[ii] = -1;
280 }
281 }
282 }
283
284 if (!data.forcesOnlyOnEntitiesCol.empty()) {
285 colIndices.resize(nbCols, false);
286 noalias(colIndices) = col_data.getIndices();
287 col_indices = &colIndices[0];
288 VectorDofs &dofs = col_data.getFieldDofs();
289 VectorDofs::iterator dit = dofs.begin();
290 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
291 if (data.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
292 data.forcesOnlyOnEntitiesCol.end()) {
293 colIndices[ii] = -1;
294 }
295 }
296 }
297
298 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
299 : getFEMethod()->snes_B;
300 // assemble local matrix
301 CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
302 &*K.data().begin(), ADD_VALUES);
303
304 if (!isDiag && sYmm) {
305 // if not diagonal term and since global matrix is symmetric assemble
306 // transpose term.
307 transK.resize(K.size2(), K.size1(), false);
308 noalias(transK) = trans(K);
309 CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
310 &*transK.data().begin(), ADD_VALUES);
311 }
313}
314
315MoFEMErrorCode HookeElement::OpAssemble::aSsemble(EntData &row_data) {
317
318 // get pointer to first global index on row
319 const int *row_indices = &*row_data.getIndices().data().begin();
320
321 auto &data = *dataAtPts;
322 if (!data.forcesOnlyOnEntitiesRow.empty()) {
323 rowIndices.resize(nbRows, false);
324 noalias(rowIndices) = row_data.getIndices();
325 row_indices = &rowIndices[0];
326 VectorDofs &dofs = row_data.getFieldDofs();
327 VectorDofs::iterator dit = dofs.begin();
328 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
329 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
330 data.forcesOnlyOnEntitiesRow.end()) {
331 rowIndices[ii] = -1;
332 }
333 }
334 }
335
336 Vec F = getFEMethod()->snes_f;
337 // assemble local matrix
338 CHKERR VecSetValues(F, nbRows, row_indices, &*nF.data().begin(), ADD_VALUES);
340}
341
342HookeElement::OpRhs_dx::OpRhs_dx(
343 const std::string row_field, const std::string col_field,
344 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
345 : OpAssemble(row_field, col_field, data_at_pts, OPROW) {}
346
347MoFEMErrorCode HookeElement::OpRhs_dx::iNtegrate(EntData &row_data) {
349
350 auto get_tensor1 = [](VectorDouble &v, const int r) {
352 &v(r + 0), &v(r + 1), &v(r + 2));
353 };
354
359
360 // get element volume
361 double vol = getVolume();
362 // get intergrayion weights
363 auto t_w = getFTensor0IntegrationWeight();
364
365 // get derivatives of base functions on rows
366 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
367 const int row_nb_base_fun = row_data.getN().size2();
368 auto t_cauchy_stress =
369 getFTensor2SymmetricFromMat<3>(*dataAtPts->cauchyStressMat);
370
371 // iterate over integration points
372 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
373
374 // calculate scalar weight times element volume
375 double a = t_w * vol;
376 auto t_nf = get_tensor1(nF, 0);
377
378 int rr = 0;
379 for (; rr != nbRows / 3; ++rr) {
380 t_nf(i) += a * t_row_diff_base(j) * t_cauchy_stress(i, j);
381 ++t_row_diff_base;
382 ++t_nf;
383 }
384
385 for (; rr != row_nb_base_fun; ++rr)
386 ++t_row_diff_base;
387
388 ++t_w;
389 ++t_cauchy_stress;
390 }
391
393}
394
395HookeElement::OpAleRhs_dx::OpAleRhs_dx(
396 const std::string row_field, const std::string col_field,
397 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
398 : OpAssemble(row_field, col_field, data_at_pts, OPROW) {}
399
400MoFEMErrorCode HookeElement::OpAleRhs_dx::iNtegrate(EntData &row_data) {
402
403 auto get_tensor1 = [](VectorDouble &v, const int r) {
405 &v(r + 0), &v(r + 1), &v(r + 2));
406 };
407
410
411 // get element volume
412 double vol = getVolume();
413 // get intergrayion weights
414 auto t_w = getFTensor0IntegrationWeight();
415
416 // get derivatives of base functions on rows
417 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
418 const int row_nb_base_fun = row_data.getN().size2();
419 auto t_cauchy_stress =
420 getFTensor2SymmetricFromMat<3>(*dataAtPts->cauchyStressMat);
421 auto &det_H = *dataAtPts->detHVec;
422 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
423
424 // iterate over integration points
425 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
426
427 // calculate scalar weight times element volume
428 double a = t_w * vol * det_H[gg];
429 auto t_nf = get_tensor1(nF, 0);
430
431 int rr = 0;
432 for (; rr != nbRows / 3; ++rr) {
433 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
434 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
435 t_nf(i) += a * t_row_diff_base_pulled(j) * t_cauchy_stress(i, j);
436 ++t_row_diff_base;
437 ++t_nf;
438 }
439
440 for (; rr != row_nb_base_fun; ++rr)
441 ++t_row_diff_base;
442
443 ++t_w;
444 ++t_cauchy_stress;
445 ++t_invH;
446 }
447
449}
450
451HookeElement::OpAleRhs_dX::OpAleRhs_dX(
452 const std::string row_field, const std::string col_field,
453 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
454 : OpAssemble(row_field, col_field, data_at_pts, OPROW) {}
455
456MoFEMErrorCode HookeElement::OpAleRhs_dX::iNtegrate(EntData &row_data) {
458
459 auto get_tensor1 = [](VectorDouble &v, const int r) {
461 &v(r + 0), &v(r + 1), &v(r + 2));
462 };
463
466
467 // get element volume
468 double vol = getVolume();
469 // get intergrayion weights
470 auto t_w = getFTensor0IntegrationWeight();
471
472 // get derivatives of base functions on rows
473 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
474 const int row_nb_base_fun = row_data.getN().size2();
475 auto t_eshelby_stress =
476 getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
477 auto &det_H = *dataAtPts->detHVec;
478 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
479
480 // iterate over integration points
481 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
482
483 // calculate scalar weight times element volume
484 double a = t_w * vol * det_H[gg];
485 auto t_nf = get_tensor1(nF, 0);
486
487 int rr = 0;
488 for (; rr != nbRows / 3; ++rr) {
489 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
490 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
491 t_nf(i) += a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
492 ++t_row_diff_base;
493 ++t_nf;
494 }
495
496 for (; rr != row_nb_base_fun; ++rr)
497 ++t_row_diff_base;
498
499 ++t_w;
500 ++t_eshelby_stress;
501 ++t_invH;
502 }
503
505}
506
507MoFEMErrorCode HookeElement::setBlocks(
508 MoFEM::Interface &m_field,
509 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr) {
511
512 if (!block_sets_ptr)
513 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
514 "Pointer to block of sets is null");
515
517 m_field, BLOCKSET | MAT_ELASTICSET, it)) {
518 Mat_Elastic mydata;
519 CHKERR it->getAttributeDataStructure(mydata);
520 int id = it->getMeshsetId();
521 auto &block_data = (*block_sets_ptr)[id];
522 EntityHandle meshset = it->getMeshset();
523 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 3,
524 block_data.tEts, true);
525 block_data.iD = id;
526 block_data.E = mydata.data.Young;
527 block_data.PoissonRatio = mydata.data.Poisson;
528 }
529
531}
532
533MoFEMErrorCode HookeElement::addElasticElement(
534 MoFEM::Interface &m_field,
535 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
536 const std::string element_name, const std::string x_field,
537 const std::string X_field, const bool ale) {
539
540 if (!block_sets_ptr)
541 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
542 "Pointer to block of sets is null");
543
544 CHKERR m_field.add_finite_element(element_name, MF_ZERO);
545 CHKERR m_field.modify_finite_element_add_field_row(element_name, x_field);
546 CHKERR m_field.modify_finite_element_add_field_col(element_name, x_field);
547 CHKERR m_field.modify_finite_element_add_field_data(element_name, x_field);
548 if (m_field.check_field(X_field)) {
549 if (ale) {
550 CHKERR m_field.modify_finite_element_add_field_row(element_name, X_field);
551 CHKERR m_field.modify_finite_element_add_field_col(element_name, X_field);
552 }
553 CHKERR m_field.modify_finite_element_add_field_data(element_name, X_field);
554 }
555
556 for (auto &m : (*block_sets_ptr)) {
557 CHKERR m_field.add_ents_to_finite_element_by_dim(m.second.tEts, 3,
558 element_name);
559 }
560
562}
563
564MoFEMErrorCode HookeElement::setOperators(
565 boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr,
566 boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr,
567 boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
568 const std::string x_field, const std::string X_field, const bool ale,
569 const bool field_disp, const EntityType type,
570 boost::shared_ptr<DataAtIntegrationPts> data_at_pts) {
572
573 if (!block_sets_ptr)
574 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
575 "Pointer to block of sets is null");
576
577 if (!data_at_pts)
578 data_at_pts = boost::make_shared<DataAtIntegrationPts>();
579
580 if (fe_lhs_ptr) {
581 if (ale == PETSC_FALSE) {
582 if (type == MBPRISM) {
583 boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
584 fe_lhs_ptr->getOpPtrVector().push_back(
585 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
586 fe_lhs_ptr->getOpPtrVector().push_back(
587 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
588 }
589 fe_lhs_ptr->getOpPtrVector().push_back(
591 x_field, x_field, block_sets_ptr, data_at_pts));
592 fe_lhs_ptr->getOpPtrVector().push_back(
593 new OpLhs_dx_dx<0>(x_field, x_field, data_at_pts));
594 } else {
595 if (type == MBPRISM) {
596 boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
597 fe_lhs_ptr->getOpPtrVector().push_back(
598 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
599 fe_lhs_ptr->getOpPtrVector().push_back(
600 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
601 }
602 fe_lhs_ptr->getOpPtrVector().push_back(
603 new OpCalculateVectorFieldGradient<3, 3>(X_field, data_at_pts->HMat));
604 fe_lhs_ptr->getOpPtrVector().push_back(
606 x_field, x_field, block_sets_ptr, data_at_pts));
607 fe_lhs_ptr->getOpPtrVector().push_back(
608 new OpCalculateVectorFieldGradient<3, 3>(x_field, data_at_pts->hMat));
609 fe_lhs_ptr->getOpPtrVector().push_back(
610 new OpCalculateStrainAle(x_field, x_field, data_at_pts));
611 fe_lhs_ptr->getOpPtrVector().push_back(
612 new OpCalculateStress<0>(x_field, x_field, data_at_pts));
613 fe_lhs_ptr->getOpPtrVector().push_back(
614 new OpAleLhs_dx_dx<0>(x_field, x_field, data_at_pts));
615 fe_lhs_ptr->getOpPtrVector().push_back(
616 new OpAleLhs_dx_dX<0>(x_field, X_field, data_at_pts));
617 fe_lhs_ptr->getOpPtrVector().push_back(
618 new OpCalculateEnergy(X_field, X_field, data_at_pts));
619 fe_lhs_ptr->getOpPtrVector().push_back(
620 new OpCalculateEshelbyStress(X_field, X_field, data_at_pts));
621 fe_lhs_ptr->getOpPtrVector().push_back(
622 new OpAleLhs_dX_dX<0>(X_field, X_field, data_at_pts));
623 fe_lhs_ptr->getOpPtrVector().push_back(
624 new OpAleLhsPre_dX_dx<0>(X_field, x_field, data_at_pts));
625 fe_lhs_ptr->getOpPtrVector().push_back(
626 new OpAleLhs_dX_dx(X_field, x_field, data_at_pts));
627 }
628 }
629
630 if (fe_rhs_ptr) {
631
632 if (ale == PETSC_FALSE) {
633 if (type == MBPRISM) {
634 boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
635 fe_rhs_ptr->getOpPtrVector().push_back(
636 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
637 fe_rhs_ptr->getOpPtrVector().push_back(
638 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
639 }
640 fe_rhs_ptr->getOpPtrVector().push_back(
641 new OpCalculateVectorFieldGradient<3, 3>(x_field, data_at_pts->hMat));
642 fe_rhs_ptr->getOpPtrVector().push_back(
643 new OpCalculateHomogeneousStiffness<0>(x_field, x_field,
644 block_sets_ptr, data_at_pts));
645 if (field_disp) {
646 fe_rhs_ptr->getOpPtrVector().push_back(
647 new OpCalculateStrain<true>(x_field, x_field, data_at_pts));
648 } else {
649 fe_rhs_ptr->getOpPtrVector().push_back(
650 new OpCalculateStrain<false>(x_field, x_field, data_at_pts));
651 }
652 fe_rhs_ptr->getOpPtrVector().push_back(
653 new OpCalculateStress<0>(x_field, x_field, data_at_pts));
654 fe_rhs_ptr->getOpPtrVector().push_back(
655 new OpRhs_dx(x_field, x_field, data_at_pts));
656 } else {
657 if (type == MBPRISM) {
658 boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
659 fe_rhs_ptr->getOpPtrVector().push_back(
660 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
661 fe_rhs_ptr->getOpPtrVector().push_back(
662 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
663 }
664 fe_rhs_ptr->getOpPtrVector().push_back(
665 new OpCalculateVectorFieldGradient<3, 3>(X_field, data_at_pts->HMat));
666 fe_rhs_ptr->getOpPtrVector().push_back(
667 new OpCalculateHomogeneousStiffness<0>(x_field, x_field,
668 block_sets_ptr, data_at_pts));
669 fe_rhs_ptr->getOpPtrVector().push_back(
670 new OpCalculateVectorFieldGradient<3, 3>(x_field, data_at_pts->hMat));
671 fe_rhs_ptr->getOpPtrVector().push_back(
672 new OpCalculateStrainAle(x_field, x_field, data_at_pts));
673 fe_rhs_ptr->getOpPtrVector().push_back(
674 new OpCalculateStress<0>(x_field, x_field, data_at_pts));
675 fe_rhs_ptr->getOpPtrVector().push_back(
676 new OpAleRhs_dx(x_field, x_field, data_at_pts));
677 fe_rhs_ptr->getOpPtrVector().push_back(
678 new OpCalculateEnergy(X_field, X_field, data_at_pts));
679 fe_rhs_ptr->getOpPtrVector().push_back(
680 new OpCalculateEshelbyStress(X_field, X_field, data_at_pts));
681 fe_rhs_ptr->getOpPtrVector().push_back(
682 new OpAleRhs_dX(X_field, X_field, data_at_pts));
683 }
684 }
685
687}
688
689MoFEMErrorCode HookeElement::calculateEnergy(
690 DM dm, boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
691 const std::string x_field, const std::string X_field, const bool ale,
692 const bool field_disp, SmartPetscObj<Vec> &v_energy) {
694
695 MoFEM::Interface *m_field_ptr;
696 CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
697
698 v_energy = createSmartVectorMPI(m_field_ptr->get_comm(), PETSC_DECIDE, 1);
699
700 boost::shared_ptr<DataAtIntegrationPts> data_at_pts(
702
703 auto fe_ptr =
704 boost::make_shared<VolumeElementForcesAndSourcesCore>(*m_field_ptr);
705 fe_ptr->getRuleHook = [](const double, const double, const double o) {
706 return 2 * o;
707 };
708
709 if (m_field_ptr->check_field("MESH_NODE_POSITIONS")) CHKERR addHOOpsVol(
710 "MESH_NODE_POSITIONS", *fe_ptr, true, false, false, false);
711
713
715 FatPrismElementForcesAndSourcesCore;
716 int getRuleTrianglesOnly(int order) { return 2 * order; }
717 int getRuleThroughThickness(int order) { return 2 * order; }
718 };
719
720 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_ptr(
721 new PrismFE(*m_field_ptr));
722
723 auto push_ops = [&](boost::shared_ptr<ForcesAndSourcesCore> fe_ptr,
726 boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
727 if (ale == PETSC_FALSE) {
728 if (type == MBPRISM) {
729 fe_ptr->getOpPtrVector().push_back(
730 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
731 fe_ptr->getOpPtrVector().push_back(
732 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
733 }
734 fe_ptr->getOpPtrVector().push_back(
735 new OpCalculateVectorFieldGradient<3, 3>(x_field, data_at_pts->hMat));
736 fe_ptr->getOpPtrVector().push_back(new OpCalculateHomogeneousStiffness<0>(
737 x_field, x_field, block_sets_ptr, data_at_pts));
738 if (field_disp) {
739 fe_ptr->getOpPtrVector().push_back(
740 new OpCalculateStrain<true>(x_field, x_field, data_at_pts));
741 } else {
742 fe_ptr->getOpPtrVector().push_back(
743 new OpCalculateStrain<false>(x_field, x_field, data_at_pts));
744 }
745 fe_ptr->getOpPtrVector().push_back(
746 new OpCalculateStress<0>(x_field, x_field, data_at_pts));
747 fe_ptr->getOpPtrVector().push_back(
748 new OpCalculateEnergy(X_field, X_field, data_at_pts, v_energy));
749 } else {
750 if (type == MBPRISM) {
751 fe_ptr->getOpPtrVector().push_back(
752 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
753 fe_ptr->getOpPtrVector().push_back(
754 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
755 }
756 fe_ptr->getOpPtrVector().push_back(
757 new OpCalculateVectorFieldGradient<3, 3>(X_field, data_at_pts->HMat));
758 fe_ptr->getOpPtrVector().push_back(new OpCalculateHomogeneousStiffness<0>(
759 x_field, x_field, block_sets_ptr, data_at_pts));
760 fe_ptr->getOpPtrVector().push_back(
761 new OpCalculateVectorFieldGradient<3, 3>(x_field, data_at_pts->hMat));
762 fe_ptr->getOpPtrVector().push_back(
763 new OpCalculateStrainAle(x_field, x_field, data_at_pts));
764 fe_ptr->getOpPtrVector().push_back(
765 new OpCalculateStress<0>(x_field, x_field, data_at_pts));
766 fe_ptr->getOpPtrVector().push_back(
767 new OpCalculateEnergy(X_field, X_field, data_at_pts, v_energy));
768 }
770 };
771
772 CHKERR push_ops(fe_ptr, MBTET);
773 CHKERR push_ops(prism_fe_ptr, MBPRISM);
774
775 CHKERR VecZeroEntries(v_energy);
776
777 fe_ptr->snes_ctx = SnesMethod::CTX_SNESNONE;
778 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_ptr);
779 CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_ptr);
780
781 CHKERR VecAssemblyBegin(v_energy);
782 CHKERR VecAssemblyEnd(v_energy);
783
785}
786
787MoFEMErrorCode HookeElement::OpAleLhs_dX_dx::iNtegrate(EntData &row_data,
788 EntData &col_data) {
790
791 // get sub-block (3x3) of local stiffens matrix, here represented by
792 // second order tensor
793 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
795 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
796 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
797 &m(r + 2, c + 2));
798 };
799
806
807 // get element volume
808 double vol = getVolume();
809
810 // get intergrayion weights
811 auto t_w = getFTensor0IntegrationWeight();
812
813 // get derivatives of base functions on rows
814 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
815 const int row_nb_base_fun = row_data.getN().size2();
816
817 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
818 auto &det_H = *dataAtPts->detHVec;
819
820 auto get_eshelby_stress_dx = [this]() {
822 t_eshelby_stress_dx;
823 int mm = 0;
824 for (int ii = 0; ii != 3; ++ii)
825 for (int jj = 0; jj != 3; ++jj)
826 for (int kk = 0; kk != 3; ++kk)
827 for (int ll = 0; ll != 3; ++ll)
828 t_eshelby_stress_dx.ptr(ii, jj, kk, ll) =
829 &(*dataAtPts->eshelbyStress_dx)(mm++, 0);
830 return t_eshelby_stress_dx;
831 };
832
833 auto t_eshelby_stress_dx = get_eshelby_stress_dx();
834
835 // iterate over integration points
836 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
837
838 // calculate scalar weight times element volume
839 double a = t_w * vol * det_H[gg];
840
841 // iterate over row base functions
842 int rr = 0;
843 for (; rr != nbRows / 3; ++rr) {
844
845 // get sub matrix for the row
846 auto t_m = get_tensor2(K, 3 * rr, 0);
847
848 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
849 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
850
851 FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dx;
852 t_row_stress_dx(i, k, l) =
853 a * t_row_diff_base_pulled(j) * t_eshelby_stress_dx(i, j, k, l);
854
855 // get derivatives of base functions for columns
856 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
857
858 // iterate column base functions
859 for (int cc = 0; cc != nbCols / 3; ++cc) {
860
861 t_m(i, k) += t_row_stress_dx(i, k, l) * t_col_diff_base(l);
862
863 // move to next column base function
864 ++t_col_diff_base;
865
866 // move to next block of local stiffens matrix
867 ++t_m;
868 }
869
870 // move to next row base function
871 ++t_row_diff_base;
872 }
873
874 for (; rr != row_nb_base_fun; ++rr)
875 ++t_row_diff_base;
876
877 ++t_w;
878 ++t_invH;
879 ++t_eshelby_stress_dx;
880 }
881
883}
884
885HookeElement::OpAleLhsWithDensity_dX_dX::OpAleLhsWithDensity_dX_dX(
886 const std::string row_field, const std::string col_field,
887 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
888 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
889 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts, const double rho_n,
890 const double rho_0)
891 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false),
892 rhoAtGaussPtsPtr(rho_at_gauss_pts),
893 rhoGradAtGaussPtsPtr(rho_grad_at_gauss_pts), rhoN(rho_n), rHo0(rho_0) {}
894
896HookeElement::OpAleLhsWithDensity_dX_dX::iNtegrate(EntData &row_data,
897 EntData &col_data) {
899
900 // get sub-block (3x3) of local stiffens matrix, here represented by
901 // second order tensor
902 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
904 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
905 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
906 &m(r + 2, c + 2));
907 };
908
913
914 // get element volume
915 double vol = getVolume();
916
917 // get intergrayion weights
918 auto t_w = getFTensor0IntegrationWeight();
919
920 // get derivatives of base functions on rows
921 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
922 const int row_nb_base_fun = row_data.getN().size2();
923
924 // Elastic stiffness tensor (4th rank tensor with minor and major
925 // symmetry)
926 auto rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
927 auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
928
929 auto t_eshelby_stress =
930 getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
931 // auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
932 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
933 auto &det_H = *dataAtPts->detHVec;
934
935 // iterate over integration points
936 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
937
938 // calculate scalar weight times element volume
939 double a = t_w * vol * det_H[gg];
940
941 const double stress_dho_coef = (rhoN / rho);
942 // (rhoN / rHo0) * pow(rho / rHo0, rhoN - 1.) * (1. / pow(rho / rHo0,
943 // rhoN));
944
945 // iterate over row base functions
946 int rr = 0;
947 for (; rr != nbRows / 3; ++rr) {
948
949 // get sub matrix for the row
950 auto t_m = get_tensor2(K, 3 * rr, 0);
951
952 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
953 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
954
955 FTensor::Tensor1<double, 3> t_row_stress;
956 t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
957
958 // get derivatives of base functions for columns
959 // auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
960 auto t_col_base = col_data.getFTensor0N(gg, 0);
961
962 // iterate column base functions
963 for (int cc = 0; cc != nbCols / 3; ++cc) {
964
965 t_m(i, k) +=
966 t_row_stress(i) * stress_dho_coef * t_grad_rho(k) * t_col_base;
967 // move to next column base function
968 ++t_col_base;
969
970 // move to next block of local stiffens matrix
971 ++t_m;
972 }
973
974 // move to next row base function
975 ++t_row_diff_base;
976 }
977
978 for (; rr != row_nb_base_fun; ++rr)
979 ++t_row_diff_base;
980
981 // move to next integration weight
982 ++t_w;
983 ++t_eshelby_stress;
984 ++t_invH;
985 ++rho;
986 ++t_grad_rho;
987 }
988
990}
991
992HookeElement::OpAleLhsWithDensity_dx_dX::OpAleLhsWithDensity_dx_dX(
993 const std::string row_field, const std::string col_field,
994 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
995 boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
996 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts, const double rho_n,
997 const double rho_0)
998 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false),
999 rhoAtGaussPtsPtr(rho_at_gauss_pts),
1000 rhoGradAtGaussPtsPtr(rho_grad_at_gauss_pts), rhoN(rho_n), rHo0(rho_0) {}
1001
1003HookeElement::OpAleLhsWithDensity_dx_dX::iNtegrate(EntData &row_data,
1004 EntData &col_data) {
1006
1007 // get sub-block (3x3) of local stiffens matrix, here represented by
1008 // second order tensor
1009 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1011 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1012 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1013 &m(r + 2, c + 2));
1014 };
1015
1020
1021 // get element volume
1022 double vol = getVolume();
1023
1024 // get integration weights
1025 auto t_w = getFTensor0IntegrationWeight();
1026
1027 // get derivatives of base functions on rows
1028 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1029 const int row_nb_base_fun = row_data.getN().size2();
1030
1031 auto rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
1032 auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
1033 auto t_cauchy_stress =
1034 getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1035 // auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1036 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1037 auto &det_H = *dataAtPts->detHVec;
1038
1039 // iterate over integration points
1040 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1041
1042 // calculate scalar weight times element volume
1043 double a = t_w * vol * det_H[gg];
1044
1045 const double stress_dho_coef = (rhoN / rho);
1046 // (rhoN / rHo0) * pow(rho / rHo0, rhoN - 1.) * (1. / pow(rho / rHo0,
1047 // rhoN)); iterate over row base functions
1048 int rr = 0;
1049 for (; rr != nbRows / 3; ++rr) {
1050
1051 // get sub matrix for the row
1052 auto t_m = get_tensor2(K, 3 * rr, 0);
1053
1054 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1055 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1056
1057 FTensor::Tensor1<double, 3> t_row_stress;
1058 t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_cauchy_stress(i, j);
1059
1060 // get derivatives of base functions for columns
1061 auto t_col_base = col_data.getFTensor0N(gg, 0);
1062 // iterate column base functions
1063 for (int cc = 0; cc != nbCols / 3; ++cc) {
1064
1065 t_m(i, k) +=
1066 t_row_stress(i) * stress_dho_coef * t_grad_rho(k) * t_col_base;
1067
1068 ++t_col_base;
1069
1070 // move to next block of local stiffens matrix
1071 ++t_m;
1072 }
1073
1074 // move to next row base function
1075 ++t_row_diff_base;
1076 }
1077
1078 for (; rr != row_nb_base_fun; ++rr)
1079 ++t_row_diff_base;
1080
1081 // move to next integration weight
1082 ++t_w;
1083 // ++t_D;
1084 ++t_cauchy_stress;
1085 ++t_invH;
1086 // ++t_h;
1087 ++rho;
1088 ++t_grad_rho;
1089 }
1090
1092}
1093
1094HookeElement::OpCalculateStiffnessScaledByDensityField::
1095 OpCalculateStiffnessScaledByDensityField(
1096 const std::string row_field, const std::string col_field,
1097 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
1098 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1099 boost::shared_ptr<VectorDouble> rho_at_gauss_pts, const double rho_n,
1100 const double rho_0)
1101
1102 : VolUserDataOperator(row_field, col_field, OPROW, true),
1103 blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts),
1104 rhoAtGaussPtsPtr(rho_at_gauss_pts), rhoN(rho_n), rHo0(rho_0) {
1105 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1106}
1107
1108MoFEMErrorCode HookeElement::OpCalculateStiffnessScaledByDensityField::doWork(
1109 int row_side, EntityType row_type, EntData &row_data) {
1111
1112 if (!rhoAtGaussPtsPtr)
1113 SETERRQ(PETSC_COMM_SELF, 1, "Calculate density with MWLS first.");
1114
1115 for (auto &m : (*blockSetsPtr)) {
1116
1117 if (m.second.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1118 m.second.tEts.end()) {
1119 continue;
1120 }
1121
1122 const int nb_integration_pts = getGaussPts().size2();
1123 dataAtPts->stiffnessMat->resize(36, nb_integration_pts, false);
1124
1126 MAT_TO_DDG(dataAtPts->stiffnessMat));
1127 const double young = m.second.E;
1128 const double poisson = m.second.PoissonRatio;
1129
1130 auto rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
1131
1132 // coefficient used in intermediate calculation
1133 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
1134
1139
1140 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1141
1142 t_D(i, j, k, l) = 0.;
1143
1144 t_D(0, 0, 0, 0) = 1 - poisson;
1145 t_D(1, 1, 1, 1) = 1 - poisson;
1146 t_D(2, 2, 2, 2) = 1 - poisson;
1147
1148 t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * poisson);
1149 t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * poisson);
1150 t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * poisson);
1151
1152 t_D(0, 0, 1, 1) = poisson;
1153 t_D(1, 1, 0, 0) = poisson;
1154 t_D(0, 0, 2, 2) = poisson;
1155 t_D(2, 2, 0, 0) = poisson;
1156 t_D(1, 1, 2, 2) = poisson;
1157 t_D(2, 2, 1, 1) = poisson;
1158 // here the coefficient is modified to take density into account for
1159 // porous materials: E(p) = E * (p / p_0)^n
1160 t_D(i, j, k, l) *= coefficient * pow(rho / rHo0, rhoN);
1161
1162 ++t_D;
1163 ++rho;
1164 }
1165 }
1166
1168}
static Index< 'o', 3 > o
#define MAT_TO_DDG(SM)
constexpr double a
@ MF_ZERO
Definition: definitions.h:98
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:533
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMMoFEM.cpp:361
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual bool check_field(const std::string &name) const =0
check if field is in database
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
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)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
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:1323
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1312
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
const double r
rate factor
double rho
Definition: plastic.cpp:98
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
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 VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Elastic material data structure.
Calculate inverse of jacobian for face element.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Transform local reference derivatives of shape functions to global derivatives.
intrusive_ptr for managing petsc objects