v0.14.0
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>
19 using namespace MoFEM;
20 #include <HookeElement.hpp>
21 
22 HookeElement::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 
30 MoFEMErrorCode 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 
72 HookeElement::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 
81 MoFEMErrorCode 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 
130 HookeElement::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 
138 MoFEMErrorCode 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 
169 HookeElement::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 
176 MoFEMErrorCode 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 
221 MoFEMErrorCode 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 
247 MoFEMErrorCode HookeElement::OpAssemble::iNtegrate(EntData &row_data,
248  EntData &col_data) {
250  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
252 };
253 
254 MoFEMErrorCode HookeElement::OpAssemble::iNtegrate(EntData &row_data) {
256  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
258 };
259 
260 MoFEMErrorCode 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 
315 MoFEMErrorCode 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 
342 HookeElement::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 
347 MoFEMErrorCode 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 
395 HookeElement::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 
400 MoFEMErrorCode 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 
451 HookeElement::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 
456 MoFEMErrorCode 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 
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 
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 
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 
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 = createVectorMPI(m_field_ptr->get_comm(), PETSC_DECIDE, 1);
699 
700  boost::shared_ptr<DataAtIntegrationPts> data_at_pts(
701  new DataAtIntegrationPts());
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 
714  using FatPrismElementForcesAndSourcesCore::
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,
724  EntityType type) {
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 
787 MoFEMErrorCode 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 
885 HookeElement::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 
896 HookeElement::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 
992 HookeElement::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 
1003 HookeElement::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 
1094 HookeElement::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 
1108 MoFEMErrorCode 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 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
setBlocks
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
OpLhs_dx_dx
Definition: HookeElement.hpp:410
rho
double rho
Definition: plastic.cpp:140
MoFEM::addHOOpsVol
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:674
OpAleLhsPre_dX_dx
Definition: HookeElement.hpp:536
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
MoFEM::Mat_Elastic
Elastic material data structure.
Definition: MaterialBlocks.hpp:139
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM.hpp
calculateEnergy
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)
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::OpCalculateInvJacForFatPrism
Calculate inverse of jacobian for face element.
Definition: UserDataOperators.hpp:3145
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::invertTensor3by3
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:1588
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
OpCalculateStrainAle
Definition: HookeElement.hpp:130
MoFEM::CoreInterface::add_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
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
OpAleRhs_dx
Definition: HookeElement.hpp:425
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MoFEM::OpSetInvJacH1ForFatPrism
Transform local reference derivatives of shape functions to global derivatives.
Definition: UserDataOperators.hpp:3172
OpCalculateHomogeneousStiffness
Definition: HookeElement.hpp:190
a
constexpr double a
Definition: approx_sphere.cpp:30
double
convert.type
type
Definition: convert.py:64
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:172
OpAleRhs_dX
Definition: HookeElement.hpp:512
OpCalculateEshelbyStress
Definition: HookeElement.hpp:177
OpCalculateEnergy
Definition: HookeElement.hpp:164
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
HookeElement.hpp
MoFEM::VectorDofs
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Definition: EntitiesFieldData.hpp:23
MoFEM::Mat_Elastic::data
_data_ data
Definition: MaterialBlocks.hpp:155
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
PrismFE
Definition: prism_elements_from_surface.cpp:62
MoFEM::EntitiesFieldData::EntData::getFieldDofs
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
Definition: EntitiesFieldData.hpp:1269
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
MoFEM::CoreInterface::add_ents_to_finite_element_by_dim
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
MoFEM::FatPrismElementForcesAndSourcesCore
FatPrism finite element.
Definition: FatPrismElementForcesAndSourcesCore.hpp:31
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
FTensor::Index< 'i', 3 >
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
convert.n
n
Definition: convert.py:82
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
FTensor::Tensor0
Definition: Tensor0.hpp:16
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1318
OpCalculateStrain
Definition: HookeElement.hpp:119
MoFEM::DMoFEMGetInterfacePtr
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMoFEM.cpp:414
MAT_TO_DDG
#define MAT_TO_DDG(SM)
Definition: HookeElement.hpp:142
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
FTensor::Ddg
Definition: Ddg_value.hpp:7
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
OpAleLhs_dx_dX
Definition: HookeElement.hpp:449
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
addElasticElement
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)
Definition: HookeElement.cpp:533
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
OpCalculateStress
Definition: HookeElement.hpp:153
setOperators
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)
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
OpAssemble
Definition: HookeElement.hpp:340
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
OpAleLhs_dx_dx
Definition: HookeElement.hpp:434
MoFEM::SmartPetscObj< Vec >
OpAleLhs_dX_dX
Definition: HookeElement.hpp:521
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
OpRhs_dx
Definition: HookeElement.hpp:401
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
F
@ F
Definition: free_surface.cpp:394
OpAleLhs_dX_dx
Definition: HookeElement.hpp:547