v0.15.0
Loading...
Searching...
No Matches
Classes | Macros | Typedefs | Functions | Variables
HookeElement.hpp File Reference
#include <BasicFiniteElements.hpp>

Go to the source code of this file.

Classes

struct  NonlinearElasticElement
 structure grouping operators and data used for calculation of nonlinear elastic element More...
 
struct  NonlinearElasticElement::BlockData
 data for calculation heat conductivity and heat capacity elements More...
 
struct  ConvectiveMassElement
 structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_elastic_elem More...
 
struct  ConvectiveMassElement::BlockData
 data for calculation inertia forces More...
 
struct  DataAtIntegrationPts
 
struct  OpCalculateStrain< D >
 
struct  OpCalculateStrainAle
 
struct  OpCalculateStress< S >
 
struct  OpCalculateEnergy
 
struct  OpCalculateEshelbyStress
 
struct  OpCalculateHomogeneousStiffness< S >
 
struct  OpCalculateMassMatrix
 Assemble mass matrix for elastic element TODO: CHANGE FORMULA *. More...
 
struct  OpCalculateStiffnessScaledByDensityField
 
struct  OpAssemble
 
struct  OpRhs_dx
 
struct  OpLhs_dx_dx< S >
 
struct  OpAleRhs_dx
 
struct  OpAleLhs_dx_dx< S >
 
struct  OpAleLhs_dx_dX< S >
 
struct  OpAleLhsWithDensity_dx_dX
 
struct  OpAleLhsWithDensity_dX_dX
 
struct  OpAleRhs_dX
 
struct  OpAleLhs_dX_dX< S >
 
struct  OpAleLhsPre_dX_dx< S >
 
struct  OpAleLhs_dX_dx
 
struct  OpAnalyticalInternalStrain_dx< S >
 
struct  OpAnalyticalInternalAleStrain_dX< S >
 
struct  OpAnalyticalInternalAleStrain_dx< S >
 
struct  OpPostProcHookeElement< ELEMENT >
 

Macros

#define MAT_TO_DDG(SM)
 

Typedefs

using MassBlockData = ConvectiveMassElement::BlockData
 
using EntData = EntitiesFieldData::EntData
 
using UserDataOperator = ForcesAndSourcesCore::UserDataOperator
 
using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator
 

Functions

static MoFEMErrorCode setBlocks (MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr)
 
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 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)
 

Variables

struct ConvectiveMassElement BlockData = NonlinearElasticElement::BlockData
 
MatrixDouble invJac
 

Macro Definition Documentation

◆ MAT_TO_DDG

#define MAT_TO_DDG (   SM)
Value:
&(*SM)(0, 0), &(*SM)(1, 0), &(*SM)(2, 0), &(*SM)(3, 0), &(*SM)(4, 0), \
&(*SM)(5, 0), &(*SM)(6, 0), &(*SM)(7, 0), &(*SM)(8, 0), &(*SM)(9, 0), \
&(*SM)(10, 0), &(*SM)(11, 0), &(*SM)(12, 0), &(*SM)(13, 0), \
&(*SM)(14, 0), &(*SM)(15, 0), &(*SM)(16, 0), &(*SM)(17, 0), \
&(*SM)(18, 0), &(*SM)(19, 0), &(*SM)(20, 0), &(*SM)(21, 0), \
&(*SM)(22, 0), &(*SM)(23, 0), &(*SM)(24, 0), &(*SM)(25, 0), \
&(*SM)(26, 0), &(*SM)(27, 0), &(*SM)(28, 0), &(*SM)(29, 0), \
&(*SM)(30, 0), &(*SM)(31, 0), &(*SM)(32, 0), &(*SM)(33, 0), \
&(*SM)(34, 0), &(*SM)(35, 0)
Examples
mofem/users_modules/basic_finite_elements/src/HookeElement.hpp, and mofem/users_modules/basic_finite_elements/src/impl/HookeElement.cpp.

Definition at line 142 of file HookeElement.hpp.

152 : public VolUserDataOperator {
153
154 OpCalculateStress(const std::string row_field, const std::string col_field,
155 boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
156
157 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
158
159 protected:
160 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
161 };
162
163 struct OpCalculateEnergy : public VolUserDataOperator {
164
165 OpCalculateEnergy(const std::string row_field, const std::string col_field,
166 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
168
169 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
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
182 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
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
196 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
197
198 protected:
199 boost::shared_ptr<map<int, BlockData>>
200 blockSetsPtr; ///< Structure keeping data about problem, like
201 ///< material parameters
202 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
203 };
204
205 /** * @brief Assemble mass matrix for elastic element TODO: CHANGE FORMULA *
206 * \f[
207 * {\bf{M}} = \int\limits_\Omega
208 * \f]
209 *
210 */
212 : public VolumeElementForcesAndSourcesCore::UserDataOperator {
213
218
219 boost::shared_ptr<DataAtIntegrationPts> commonData;
220
221 OpCalculateMassMatrix(const std::string row_field,
222 const std::string col_field, BlockData &data,
223 MassBlockData &mass_data,
224 boost::shared_ptr<DataAtIntegrationPts> &common_data,
225 bool symm = true)
227 row_field, col_field, OPROWCOL, symm),
228 commonData(common_data), dAta(data), massData(mass_data) {}
229
230 PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
231 EntityType col_type,
233 EntitiesFieldData::EntData &col_data) {
235
236 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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 }
254 if (massData.tEts.find(getFEEntityHandle()) == massData.tEts.end()) {
256 }
257
258 const bool diagonal_block =
259 (row_type == col_type) && (row_side == col_side);
260 // get number of integration points
261 // Set size can clear local tangent matrix
262 locK.resize(row_nb_dofs, col_nb_dofs, false);
263 locK.clear();
264
265 const int row_nb_gauss_pts = row_data.getN().size1();
266 const int row_nb_base_functions = row_data.getN().size2();
267
268 FTensor::Index<'i', 3> i;
269 FTensor::Index<'j', 3> j;
270 FTensor::Index<'k', 3> k;
271 FTensor::Index<'l', 3> l;
272
273 double density = massData.rho0;
274
275 // get integration weights
276 auto t_w = getFTensor0IntegrationWeight();
277
278 // integrate local matrix for entity block
279 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
280
281 auto t_row_base_func = row_data.getFTensor0N(gg, 0);
282
283 // Get volume and integration weight
284 double w = getVolume() * t_w;
285
286 for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
287 auto t_col_base_func = col_data.getFTensor0N(gg, 0);
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 // Next base function for column
292 ++t_col_base_func;
293 }
294 // Next base function for row
295 ++t_row_base_func;
296 }
297 // Next integration point for getting weight
298 ++t_w;
299 }
300
301 CHKERR MatSetValues(getKSPB(), row_data, col_data, &locK(0, 0),
302 ADD_VALUES);
303
304 // is symmetric
305 if (row_type != col_type || row_side != col_side) {
306 translocK.resize(col_nb_dofs, row_nb_dofs, false);
307 noalias(translocK) = trans(locK);
308
309 CHKERR MatSetValues(getKSPB(), col_data, row_data, &translocK(0, 0),
310 ADD_VALUES);
311 }
312
314 }
315 };
316
318 protected:
319 boost::shared_ptr<map<int, BlockData>>
320 blockSetsPtr; ///< Structure keeping data about problem, like
321 ///< material parameters
322 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
323
324 boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
325 const double rhoN; ///< exponent n in E(p) = E * (p / p_0)^n
326 const double rHo0; ///< p_0 reference density in E(p) = E * (p / p_0)^n
327 // // where p is density, E - youngs modulus
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
336 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
337 };
338
339 struct OpAssemble : public VolUserDataOperator {
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 * \brief Do calculations for give operator
347 * @param row_side row side number (local number) of entity on element
348 * @param col_side column side number (local number) of entity on element
349 * @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
350 * @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
351 * @param row_data data for row
352 * @param col_data data for column
353 * @return error code
354 */
355 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
356 EntityType col_type, EntData &row_data,
357 EntData &col_data);
358
359 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
360
361 protected:
362 // Finite element stiffness sub-matrix K_ij
366
367 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
368
371
372 int nbRows; ///< number of dofs on rows
373 int nbCols; ///< number if dof on column
374 int nbIntegrationPts; ///< number of integration points
375 bool isDiag; ///< true if this block is on diagonal
376
377 virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
378
379 virtual MoFEMErrorCode iNtegrate(EntData &row_data);
380
381 /**
382 * \brief Assemble local entity block matrix
383 * @param row_data row data (consist base functions on row entity)
384 * @param col_data column data (consist base functions on column
385 * entity)
386 * @return error code
387 */
388 MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data);
389
390 /**
391 * \brief Assemble local entity right-hand vector
392 * @param row_data row data (consist base functions on row entity)
393 * @param col_data column data (consist base functions on column
394 * entity)
395 * @return error code
396 */
398 };
399
400 struct OpRhs_dx : public OpAssemble {
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
409 template <int S = 0> struct OpLhs_dx_dx : public OpAssemble {
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 * \brief Integrate B^T D B operator
417 * @param row_data row data (consist base functions on row entity)
418 * @param col_data column data (consist base functions on column entity)
419 * @return error code
420 */
421 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
422 };
423
424 struct OpAleRhs_dx : public OpAssemble {
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
433 template <int S = 0> struct OpAleLhs_dx_dx : public OpAssemble {
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 * \brief Integrate B^T D B operator
441 * @param row_data row data (consist base functions on row entity)
442 * @param col_data column data (consist base functions on column entity)
443 * @return error code
444 */
445 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
446 };
447
448 template <int S = 0> struct OpAleLhs_dx_dX : public OpAssemble {
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 * \brief Integrate tangent stiffness for spatial momentum
456 * @param row_data row data (consist base functions on row entity)
457 * @param col_data column data (consist base functions on column entity)
458 * @return error code
459 */
460 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
461 };
462
463 struct OpAleLhsWithDensity_dx_dX : public OpAssemble {
464
465 boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
466 boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
467 const double rhoN;
468 const double rHo0;
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 * \brief Integrate tangent stiffness for spatial momentum
480 * @param row_data row data (consist base functions on row entity)
481 * @param col_data column data (consist base functions on column entity)
482 * @return error code
483 */
484 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
485 };
486
487 struct OpAleLhsWithDensity_dX_dX : public OpAssemble {
488
489 boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
490 boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
491 const double rhoN;
492 const double rHo0;
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 * \brief Integrate tangent stiffness for material momentum
504 * @param row_data row data (consist base functions on row entity)
505 * @param col_data column data (consist base functions on column entity)
506 * @return error code
507 */
508 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
509 };
510
511 struct OpAleRhs_dX : public OpAssemble {
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
520 template <int S = 0> struct OpAleLhs_dX_dX : public OpAssemble {
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 * \brief Integrate tangent stiffness for material momentum
528 * @param row_data row data (consist base functions on row entity)
529 * @param col_data column data (consist base functions on column entity)
530 * @return error code
531 */
532 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
533 };
534
535 template <int S = 0> struct OpAleLhsPre_dX_dx : public VolUserDataOperator {
536
537 OpAleLhsPre_dX_dx(const std::string row_field, const std::string col_field,
538 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
539
540 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
541
542 private:
543 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
544 };
545
546 struct OpAleLhs_dX_dx : public OpAssemble {
547
548 OpAleLhs_dX_dx(const std::string row_field, const std::string col_field,
549 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
550 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
551
552 protected:
553 /**
554 * \brief Integrate tangent stiffness for material momentum
555 * @param row_data row data (consist base functions on row entity)
556 * @param col_data column data (consist base functions on column entity)
557 * @return error code
558 */
559 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
560 };
561
562 template <int S> struct OpAnalyticalInternalStrain_dx : public OpAssemble {
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,
578 StrainFunction strain_fun);
579
580 protected:
583 };
584
585 template <int S> struct OpAnalyticalInternalAleStrain_dX : public OpAssemble {
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,
601 StrainFunction strain_fun,
602 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
603
604 protected:
607 boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
608 };
609
610 template <int S> struct OpAnalyticalInternalAleStrain_dx : public OpAssemble {
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,
626 StrainFunction strain_fun,
627 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
628
629 protected:
632 boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
633 };
634
635 template <class ELEMENT>
636 struct OpPostProcHookeElement : public ELEMENT::UserDataOperator {
637 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
638 map<int, BlockData>
639 &blockSetsPtr; // FIXME: (works only with the first block)
640 moab::Interface &postProcMesh;
641 std::vector<EntityHandle> &mapGaussPts;
642 bool isALE;
643 bool isFieldDisp;
644
645 OpPostProcHookeElement(const string row_field,
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
652 MoFEMErrorCode doWork(int side, EntityType type,
654 };
655
656 static MoFEMErrorCode
658 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr);
659
660 static MoFEMErrorCode
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
666 static MoFEMErrorCode
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
675 static MoFEMErrorCode
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,
679 SmartPetscObj<Vec> &v_energy_ptr);
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)
689 : VolUserDataOperator(row_field, col_field, OPROW, true),
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,
697 EntData &row_data) {
699 FTensor::Index<'i', 3> i;
700 FTensor::Index<'j', 3> j;
701 // get number of integration points
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 // If displacement field, not field o spatial positons is given
711 if (!D) {
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>
730MoFEMErrorCode HookeElement::OpAleLhs_dx_dx<S>::iNtegrate(EntData &row_data,
731 EntData &col_data) {
733
734 // get sub-block (3x3) of local stiffens matrix, here represented by
735 // second order tensor
736 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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),
740 &m(r + 2, c + 2));
741 };
742
743 FTensor::Index<'i', 3> i;
744 FTensor::Index<'j', 3> j;
745 FTensor::Index<'k', 3> k;
746 FTensor::Index<'l', 3> l;
747
748 // get element volume
749 double vol = getVolume();
750
751 // get intergrayion weights
752 auto t_w = getFTensor0IntegrationWeight();
753
754 // get derivatives of base functions on rows
755 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
756 const int row_nb_base_fun = row_data.getN().size2();
757
758 // Elastic stiffness tensor (4th rank tensor with minor and major
759 // symmetry)
761 MAT_TO_DDG(dataAtPts->stiffnessMat));
762
763 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
764 auto &det_H = *dataAtPts->detHVec;
765
766 // iterate over integration points
767 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
768
769 // calculate scalar weight times element volume
770 double a = t_w * vol * det_H[gg];
771
772 // iterate over row base functions
773 int rr = 0;
774 for (; rr != nbRows / 3; ++rr) {
775
776 // get sub matrix for the row
777 auto t_m = get_tensor2(K, 3 * rr, 0);
778
779 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
780 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
781
783 // I mix up the indices here so that it behaves like a
784 // Dg. That way I don't have to have a separate wrapper
785 // class Christof_Expr, which simplifies things.
786 t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base_pulled(i));
787
788 // get derivatives of base functions for columns
789 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
790
791 // iterate column base functions
792 for (int cc = 0; cc != nbCols / 3; ++cc) {
793
794 FTensor::Tensor1<double, 3> t_col_diff_base_pulled;
795 t_col_diff_base_pulled(j) = t_col_diff_base(i) * t_invH(i, j);
796
797 // integrate block local stiffens matrix
798 t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base_pulled(k);
799
800 // move to next column base function
801 ++t_col_diff_base;
802
803 // move to next block of local stiffens matrix
804 ++t_m;
805 }
806
807 // move to next row base function
808 ++t_row_diff_base;
809 }
810
811 for (; rr != row_nb_base_fun; ++rr)
812 ++t_row_diff_base;
813
814 // move to next integration weight
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)
829 : VolUserDataOperator(row_field, col_field, OPROW, true),
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);
844 MAT_TO_DDG(dataAtPts->stiffnessMat));
845 const double young = m.second.E;
846 const double poisson = m.second.PoissonRatio;
847
848 // coefficient used in intermediate calculation
849 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
850
851 FTensor::Index<'i', 3> i;
852 FTensor::Index<'j', 3> j;
853 FTensor::Index<'k', 3> k;
854 FTensor::Index<'l', 3> l;
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)
886 : VolUserDataOperator(row_field, col_field, OPROW, true),
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,
894 EntData &row_data) {
896 // get number of integration points
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
903 FTensor::Index<'i', 3> i;
904 FTensor::Index<'j', 3> j;
905 FTensor::Index<'k', 3> k;
906 FTensor::Index<'l', 3> l;
907
908 // elastic stiffness tensor (4th rank tensor with minor and major
909 // symmetry)
911 MAT_TO_DDG(dataAtPts->stiffnessMat));
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>
928MoFEMErrorCode HookeElement::OpLhs_dx_dx<S>::iNtegrate(EntData &row_data,
929 EntData &col_data) {
931
932 // get sub-block (3x3) of local stiffens matrix, here represented by
933 // second order tensor
934 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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),
938 &m(r + 2, c + 2));
939 };
940
941 FTensor::Index<'i', 3> i;
942 FTensor::Index<'j', 3> j;
943 FTensor::Index<'k', 3> k;
944 FTensor::Index<'l', 3> l;
945
946 // get element volume
947 double vol = getVolume();
948
949 // get intergrayion weights
950 auto t_w = getFTensor0IntegrationWeight();
951
952 // get derivatives of base functions on rows
953 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
954 const int row_nb_base_fun = row_data.getN().size2();
955
956 // Elastic stiffness tensor (4th rank tensor with minor and major
957 // symmetry)
959 MAT_TO_DDG(dataAtPts->stiffnessMat));
960
961 // iterate over integration points
962 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
963
964 // calculate scalar weight times element volume
965 double a = t_w * vol;
966
967 // iterate over row base functions
968 int rr = 0;
969 for (; rr != nbRows / 3; ++rr) {
970
971 // get sub matrix for the row
972 auto t_m = get_tensor2(K, 3 * rr, 0);
973
974 // get derivatives of base functions for columns
975 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
976
978 // I mix up the indices here so that it behaves like a
979 // Dg. That way I don't have to have a separate wrapper
980 // class Christof_Expr, which simplifies things.
981 t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
982
983 // iterate column base functions
984 for (int cc = 0; cc != nbCols / 3; ++cc) {
985
986 // integrate block local stiffens matrix
987 t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
988
989 // move to next column base function
990 ++t_col_diff_base;
991
992 // move to next block of local stiffens matrix
993 ++t_m;
994 }
995
996 // move to next row base function
997 ++t_row_diff_base;
998 }
999
1000 for (; rr != row_nb_base_fun; ++rr)
1001 ++t_row_diff_base;
1002
1003 // move to next integration weight
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>
1018MoFEMErrorCode HookeElement::OpAleLhs_dx_dX<S>::iNtegrate(EntData &row_data,
1019 EntData &col_data) {
1021
1022 // get sub-block (3x3) of local stiffens matrix, here represented by
1023 // second order tensor
1024 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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),
1028 &m(r + 2, c + 2));
1029 };
1030
1031 FTensor::Index<'i', 3> i;
1032 FTensor::Index<'j', 3> j;
1033 FTensor::Index<'k', 3> k;
1034 FTensor::Index<'l', 3> l;
1035 FTensor::Index<'m', 3> m;
1036 FTensor::Index<'n', 3> n;
1037
1038 // get element volume
1039 double vol = getVolume();
1040
1041 // get intergrayion weights
1042 auto t_w = getFTensor0IntegrationWeight();
1043
1044 // get derivatives of base functions on rows
1045 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1046 const int row_nb_base_fun = row_data.getN().size2();
1047
1048 // Elastic stiffness tensor (4th rank tensor with minor and major
1049 // symmetry)
1051 MAT_TO_DDG(dataAtPts->stiffnessMat));
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 // iterate over integration points
1060 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1061
1062 // calculate scalar weight times element volume
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 // iterate over row base functions
1069 int rr = 0;
1070 for (; rr != nbRows / 3; ++rr) {
1071
1072 // get sub matrix for the row
1073 auto t_m = get_tensor2(K, 3 * rr, 0);
1074
1075 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1076 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1077
1078 FTensor::Tensor1<double, 3> t_row_stress;
1079 t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_cauchy_stress(i, j);
1080
1081 FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
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
1085 FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
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
1092 FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1093 // FIXME: This operator is not implemented, doing operation by hand
1094 // t_row_stress_dX(i, m, n) = t_row_D(i, k, l) * t_F_dX(k, l, m, n);
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 // get derivatives of base functions for columns
1106 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1107
1108 // iterate column base functions
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 // move to next column base function
1116 ++t_col_diff_base;
1117
1118 // move to next block of local stiffens matrix
1119 ++t_m;
1120 }
1121
1122 // move to next row base function
1123 ++t_row_diff_base;
1124 }
1125
1126 for (; rr != row_nb_base_fun; ++rr)
1127 ++t_row_diff_base;
1128
1129 // move to next integration weight
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>
1147MoFEMErrorCode HookeElement::OpAleLhs_dX_dX<S>::iNtegrate(EntData &row_data,
1148 EntData &col_data) {
1150
1151 // get sub-block (3x3) of local stiffens matrix, here represented by
1152 // second order tensor
1153 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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),
1157 &m(r + 2, c + 2));
1158 };
1159
1160 FTensor::Index<'i', 3> i;
1161 FTensor::Index<'j', 3> j;
1162 FTensor::Index<'k', 3> k;
1163 FTensor::Index<'l', 3> l;
1164 FTensor::Index<'m', 3> m;
1165 FTensor::Index<'n', 3> n;
1166
1167 // get element volume
1168 double vol = getVolume();
1169
1170 // get intergrayion weights
1171 auto t_w = getFTensor0IntegrationWeight();
1172
1173 // get derivatives of base functions on rows
1174 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1175 const int row_nb_base_fun = row_data.getN().size2();
1176
1177 // Elastic stiffness tensor (4th rank tensor with minor and major
1178 // symmetry)
1180 MAT_TO_DDG(dataAtPts->stiffnessMat));
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 // iterate over integration points
1192 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1193
1194 // calculate scalar weight times element volume
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
1212 FTensor::Tensor4<double, 3, 3, 3, 3> t_eshelby_stress_dX;
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 // iterate over row base functions
1240 int rr = 0;
1241 for (; rr != nbRows / 3; ++rr) {
1242
1243 // get sub matrix for the row
1244 auto t_m = get_tensor2(K, 3 * rr, 0);
1245
1246 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1247 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1248
1249 FTensor::Tensor1<double, 3> t_row_stress;
1250 t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1251
1252 FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
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
1256 FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
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
1260 FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
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 // get derivatives of base functions for columns
1265 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1266
1267 // iterate column base functions
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 // move to next column base function
1275 ++t_col_diff_base;
1276
1277 // move to next block of local stiffens matrix
1278 ++t_m;
1279 }
1280
1281 // move to next row base function
1282 ++t_row_diff_base;
1283 }
1284
1285 for (; rr != row_nb_base_fun; ++rr)
1286 ++t_row_diff_base;
1287
1288 // move to next integration weight
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)
1306 : VolUserDataOperator(row_field, col_field, OPROW, true),
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,
1314 EntData &row_data) {
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
1335 FTensor::Index<'i', 3> i;
1336 FTensor::Index<'j', 3> j;
1337 FTensor::Index<'k', 3> k;
1338 FTensor::Index<'l', 3> l;
1339 FTensor::Index<'m', 3> m;
1340 FTensor::Index<'n', 3> n;
1341
1342 // Elastic stiffness tensor (4th rank tensor with minor and major
1343 // symmetry)
1345 MAT_TO_DDG(dataAtPts->stiffnessMat));
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)
1389 : ELEMENT::UserDataOperator(row_field, UserDataOperator::OPROW),
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(
1396 int side, EntityType type, EntitiesFieldData::EntData &data) {
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) {
1416 Tag th;
1417 CHKERR postProcMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
1418 MB_TAG_CREAT | MB_TAG_SPARSE,
1419 def_val.data());
1420 return th;
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
1428 FTensor::Index<'i', 3> i;
1429 FTensor::Index<'j', 3> j;
1430 FTensor::Index<'k', 3> k;
1431 FTensor::Index<'l', 3> l;
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);
1438 MAT_TO_DDG(dataAtPts->stiffnessMat));
1439
1440 EntityHandle ent = this->getFEEntityHandle();
1441 auto type = type_from_handle(ent);
1442 EntityHandle ent_3d = ent;
1443 if (type == MBTRI || type == MBQUAD) {
1444 Range ents;
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())
1450 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
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)
1478 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1479 "Element not found in any of material blocksets");
1480
1481 int def_val_int = 0;
1482 Tag tag_mat;
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.;
1489 FTensor::Tensor2<double, 3, 3> t_small_strain;
1491 FTensor::Tensor2_symmetric<double, 3> t_small_strain_symm;
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 {
1504 CHKERR determinantTensor3by3(t_H, detH);
1505 CHKERR invertTensor3by3(t_H, detH, t_invH);
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 // symmetric tensors need improvement
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) {
1544 FTensor::Index<'i', 3> i;
1545 FTensor::Index<'j', 3> j;
1546 FTensor::Index<'k', 3> k;
1547 FTensor::Index<'l', 3> l;
1549
1550 auto get_tensor1 = [](VectorDouble &v, const int r) {
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 // get element volume
1559 double vol = getVolume();
1560 auto t_w = getFTensor0IntegrationWeight();
1561
1562 nF.resize(nbRows, false);
1563 nF.clear();
1564
1565 // elastic stiffness tensor (4th rank tensor with minor and major
1566 // symmetry)
1568 MAT_TO_DDG(dataAtPts->stiffnessMat));
1569
1570 // get derivatives of base functions on rows
1571 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
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 // calculate scalar weight times element volume
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(
1615 EntData &row_data) {
1616 FTensor::Index<'i', 3> i;
1617 FTensor::Index<'j', 3> j;
1618 FTensor::Index<'k', 3> k;
1619 FTensor::Index<'l', 3> l;
1621
1622 auto get_tensor1 = [](VectorDouble &v, const int r) {
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 // get element volume
1633 double vol = getVolume();
1634 auto t_w = getFTensor0IntegrationWeight();
1635
1636 nF.resize(nbRows, false);
1637 nF.clear();
1638
1639 // elastic stiffness tensor (4th rank tensor with minor and major
1640 // symmetry)
1642 MAT_TO_DDG(dataAtPts->stiffnessMat));
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 // get derivatives of base functions on rows
1648 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
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);
1656 FTensor::Tensor2<double, 3, 3> t_eshelby_stress;
1657 t_eshelby_stress(i, j) = -t_F(k, i) * t_stress(k, j);
1658
1659 // calculate scalar weight times element volume
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) {
1666 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
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(
1698 EntData &row_data) {
1699 FTensor::Index<'i', 3> i;
1700 FTensor::Index<'j', 3> j;
1701 FTensor::Index<'k', 3> k;
1702 FTensor::Index<'l', 3> l;
1704
1705 auto get_tensor1 = [](VectorDouble &v, const int r) {
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 // get element volume
1716 double vol = getVolume();
1717 auto t_w = getFTensor0IntegrationWeight();
1718
1719 nF.resize(nbRows, false);
1720 nF.clear();
1721
1722 // elastic stiffness tensor (4th rank tensor with minor and major
1723 // symmetry)
1725 MAT_TO_DDG(dataAtPts->stiffnessMat));
1726 auto &det_H = *dataAtPts->detHVec;
1727 auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1728
1729 // get derivatives of base functions on rows
1730 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
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 // calculate scalar weight times element volume
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) {
1746 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
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 // __HOOKE_ELEMENT_HPP
static MoFEMErrorCode addElasticElement(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr, const std::string element_name, const std::string x_field, const std::string X_field, const bool ale)
static MoFEMErrorCode setOperators(boost::shared_ptr< ForcesAndSourcesCore > fe_lhs_ptr, boost::shared_ptr< ForcesAndSourcesCore > fe_rhs_ptr, boost::shared_ptr< map< int, BlockData > > block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, const EntityType type=MBTET, boost::shared_ptr< DataAtIntegrationPts > data_at_pts=nullptr)
#define MAT_TO_DDG(SM)
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr)
MatrixDouble invJac
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData > > block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, SmartPetscObj< Vec > &v_energy_ptr)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
double D
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasVector< int > VectorInt
Definition Types.hpp:67
auto type_from_handle(const EntityHandle h)
get type from entity handle
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
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.
int r
Definition sdf.py:8
FTensor::Index< 'm', 3 > m
data for calculation inertia forces
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
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
VectorDouble nF
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Assemble local entity block matrix.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Do calculations for give operator.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
int nbIntegrationPts
number of integration points
int nbRows
number of dofs on rows
MatrixDouble transK
VectorInt colIndices
VectorInt rowIndices
MatrixDouble K
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)
MassBlockData & massData
boost::shared_ptr< DataAtIntegrationPts > commonData
const double rHo0
p_0 reference density in E(p) = E * (p / p_0)^n
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
const double rhoN
exponent n in E(p) = E * (p / p_0)^n
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)

Typedef Documentation

◆ EntData

Definition at line 74 of file HookeElement.hpp.

◆ MassBlockData

Definition at line 72 of file HookeElement.hpp.

◆ UserDataOperator

using UserDataOperator = ForcesAndSourcesCore::UserDataOperator

Definition at line 75 of file HookeElement.hpp.

◆ VolUserDataOperator

Definition at line 76 of file HookeElement.hpp.

Function Documentation

◆ 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 
)
static

◆ 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 
)
static

◆ setBlocks()

static MoFEMErrorCode setBlocks ( MoFEM::Interface m_field,
boost::shared_ptr< map< int, BlockData > > &  block_sets_ptr 
)
static

◆ 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 
)
static

Variable Documentation

◆ invJac

MatrixDouble invJac
private