v0.16.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,
167 SmartPetscObj<Vec> ghost_vec = SmartPetscObj<Vec>());
168
169 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
170
171 protected:
172 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
173 SmartPetscObj<Vec> ghostVec;
174 };
175
176 struct OpCalculateEshelbyStress : public VolUserDataOperator {
177
178 OpCalculateEshelbyStress(
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>
189 struct OpCalculateHomogeneousStiffness : public VolUserDataOperator {
190
191 OpCalculateHomogeneousStiffness(
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 */
211 struct OpCalculateMassMatrix
212 : public VolumeElementForcesAndSourcesCore::UserDataOperator {
213
214 MatrixDouble locK;
215 MatrixDouble translocK;
216 BlockData &dAta;
217 MassBlockData &massData;
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)
226 : VolumeElementForcesAndSourcesCore::UserDataOperator(
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,
232 EntitiesFieldData::EntData &row_data,
233 EntitiesFieldData::EntData &col_data) {
234 MoFEMFunctionBegin;
235
236 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
237 return FTensor::Tensor2<double *, 3, 3>(
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)
247 MoFEMFunctionReturnHot(0);
248 const int col_nb_dofs = col_data.getIndices().size();
249 if (!col_nb_dofs)
250 MoFEMFunctionReturnHot(0);
251 if (dAta.tEts.find(getFEEntityHandle()) == dAta.tEts.end()) {
252 MoFEMFunctionReturnHot(0);
253 }
254 if (massData.tEts.find(getFEEntityHandle()) == massData.tEts.end()) {
255 MoFEMFunctionReturnHot(0);
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
313 MoFEMFunctionReturn(0);
314 }
315 };
316
317 struct OpCalculateStiffnessScaledByDensityField : public VolUserDataOperator {
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:
329 OpCalculateStiffnessScaledByDensityField(
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
363 MatrixDouble K;
364 MatrixDouble transK;
365 VectorDouble nF;
366
367 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
368
369 VectorInt rowIndices;
370 VectorInt colIndices;
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 */
397 MoFEMErrorCode aSsemble(EntData &row_data);
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:
406 MoFEMErrorCode iNtegrate(EntData &row_data);
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:
430 MoFEMErrorCode iNtegrate(EntData &row_data);
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
470 OpAleLhsWithDensity_dx_dX(
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
494 OpAleLhsWithDensity_dX_dX(
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:
517 MoFEMErrorCode iNtegrate(EntData &row_data);
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
566 FTensor::Tensor2_symmetric<double, 3>(
567
568 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> &t_coords
569
570 )
571
572 >
573 StrainFunction;
574
575 OpAnalyticalInternalStrain_dx(
576 const std::string row_field,
577 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
578 StrainFunction strain_fun);
579
580 protected:
581 MoFEMErrorCode iNtegrate(EntData &row_data);
582 StrainFunction strainFun;
583 };
584
585 template <int S> struct OpAnalyticalInternalAleStrain_dX : public OpAssemble {
586
587 typedef boost::function<
588
589 FTensor::Tensor2_symmetric<double, 3>(
590
591 FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> &t_coords
592
593 )
594
595 >
596 StrainFunction;
597
598 OpAnalyticalInternalAleStrain_dX(
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:
605 MoFEMErrorCode iNtegrate(EntData &row_data);
606 StrainFunction strainFun;
607 boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
608 };
609
610 template <int S> struct OpAnalyticalInternalAleStrain_dx : public OpAssemble {
611
612 typedef boost::function<
613
614 FTensor::Tensor2_symmetric<double, 3>(
615
616 FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> &t_coords
617
618 )
619
620 >
621 StrainFunction;
622
623 OpAnalyticalInternalAleStrain_dx(
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:
630 MoFEMErrorCode iNtegrate(EntData &row_data);
631 StrainFunction strainFun;
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,
653 EntitiesFieldData::EntData &data);
654 };
655
656 static MoFEMErrorCode
657 setBlocks(MoFEM::Interface &m_field,
658 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr);
659
660 static MoFEMErrorCode
661 addElasticElement(MoFEM::Interface &m_field,
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:
682 MatrixDouble invJac;
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) {
698 MoFEMFunctionBegin;
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, -1, CoeffsByGauss>(
705 *(dataAtPts->smallStrainMat));
706 auto t_h =
707 getFTensor2FromMat<3, 3>(*(dataAtPts->hMat));
708
709 for (int gg = 0; gg != nb_integration_pts; ++gg) {
710 t_strain(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
711
712 // If displacement field, not field o spatial positons is given
713 if (!D) {
714 t_strain(0, 0) -= 1;
715 t_strain(1, 1) -= 1;
716 t_strain(2, 2) -= 1;
717 }
718
719 ++t_strain;
720 ++t_h;
721 }
722 MoFEMFunctionReturn(0);
723}
724
725template <int S>
726HookeElement::OpAleLhs_dx_dx<S>::OpAleLhs_dx_dx(
727 const std::string row_field, const std::string col_field,
728 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
729 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
730
731template <int S>
732MoFEMErrorCode HookeElement::OpAleLhs_dx_dx<S>::iNtegrate(EntData &row_data,
733 EntData &col_data) {
734 MoFEMFunctionBegin;
735
736 // get sub-block (3x3) of local stiffens matrix, here represented by
737 // second order tensor
738 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
739 return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
740 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
741 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
742 &m(r + 2, c + 2));
743 };
744
745 FTensor::Index<'i', 3> i;
746 FTensor::Index<'j', 3> j;
747 FTensor::Index<'k', 3> k;
748 FTensor::Index<'l', 3> l;
749
750 // get element volume
751 double vol = getVolume();
752
753 // get intergrayion weights
754 auto t_w = getFTensor0IntegrationWeight();
755
756 // get derivatives of base functions on rows
757 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
758 const int row_nb_base_fun = row_data.getN().size2();
759
760 // Elastic stiffness tensor (4th rank tensor with minor and major
761 // symmetry)
762 FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
763 MAT_TO_DDG(dataAtPts->stiffnessMat));
764
765 auto t_invH =
766 getFTensor2FromMat<3, 3, -1, CoeffsByGauss>(*dataAtPts->invHMat);
767 auto &det_H = *dataAtPts->detHVec;
768
769 // iterate over integration points
770 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
771
772 // calculate scalar weight times element volume
773 double a = t_w * vol * det_H[gg];
774
775 // iterate over row base functions
776 int rr = 0;
777 for (; rr != nbRows / 3; ++rr) {
778
779 // get sub matrix for the row
780 auto t_m = get_tensor2(K, 3 * rr, 0);
781
782 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
783 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
784
785 FTensor::Christof<double, 3, 3> t_rowD;
786 // I mix up the indices here so that it behaves like a
787 // Dg. That way I don't have to have a separate wrapper
788 // class Christof_Expr, which simplifies things.
789 t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base_pulled(i));
790
791 // get derivatives of base functions for columns
792 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
793
794 // iterate column base functions
795 for (int cc = 0; cc != nbCols / 3; ++cc) {
796
797 FTensor::Tensor1<double, 3> t_col_diff_base_pulled;
798 t_col_diff_base_pulled(j) = t_col_diff_base(i) * t_invH(i, j);
799
800 // integrate block local stiffens matrix
801 t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base_pulled(k);
802
803 // move to next column base function
804 ++t_col_diff_base;
805
806 // move to next block of local stiffens matrix
807 ++t_m;
808 }
809
810 // move to next row base function
811 ++t_row_diff_base;
812 }
813
814 for (; rr != row_nb_base_fun; ++rr)
815 ++t_row_diff_base;
816
817 // move to next integration weight
818 ++t_w;
819 ++t_D;
820 ++t_invH;
821 }
822
823 MoFEMFunctionReturn(0);
824}
825
826template <int S>
827HookeElement::OpCalculateHomogeneousStiffness<S>::
828 OpCalculateHomogeneousStiffness(
829 const std::string row_field, const std::string col_field,
830 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
831 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
832 : VolUserDataOperator(row_field, col_field, OPROW, true),
833 blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts) {
834 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
835}
836
837template <int S>
838MoFEMErrorCode HookeElement::OpCalculateHomogeneousStiffness<S>::doWork(
839 int row_side, EntityType row_type, EntData &row_data) {
840 MoFEMFunctionBegin;
841
842 for (auto &m : (*blockSetsPtr)) {
843 if (m.second.tEts.find(getFEEntityHandle()) != m.second.tEts.end()) {
844
845 dataAtPts->stiffnessMat->resize(36, 1, false);
846 FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
847 MAT_TO_DDG(dataAtPts->stiffnessMat));
848 const double young = m.second.E;
849 const double poisson = m.second.PoissonRatio;
850
851 // coefficient used in intermediate calculation
852 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
853
854 FTensor::Index<'i', 3> i;
855 FTensor::Index<'j', 3> j;
856 FTensor::Index<'k', 3> k;
857 FTensor::Index<'l', 3> l;
858
859 t_D(i, j, k, l) = 0.;
860
861 t_D(0, 0, 0, 0) = 1 - poisson;
862 t_D(1, 1, 1, 1) = 1 - poisson;
863 t_D(2, 2, 2, 2) = 1 - poisson;
864
865 t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * poisson);
866 t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * poisson);
867 t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * poisson);
868
869 t_D(0, 0, 1, 1) = poisson;
870 t_D(1, 1, 0, 0) = poisson;
871 t_D(0, 0, 2, 2) = poisson;
872 t_D(2, 2, 0, 0) = poisson;
873 t_D(1, 1, 2, 2) = poisson;
874 t_D(2, 2, 1, 1) = poisson;
875
876 t_D(i, j, k, l) *= coefficient;
877
878 break;
879 }
880 }
881
882 MoFEMFunctionReturn(0);
883}
884
885template <int S>
886HookeElement::OpCalculateStress<S>::OpCalculateStress(
887 const std::string row_field, const std::string col_field,
888 boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
889 : VolUserDataOperator(row_field, col_field, OPROW, true),
890 dataAtPts(data_at_pts) {
891 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
892}
893
894template <int S>
895MoFEMErrorCode HookeElement::OpCalculateStress<S>::doWork(int row_side,
896 EntityType row_type,
897 EntData &row_data) {
898 MoFEMFunctionBegin;
899 // get number of integration points
900 const int nb_integration_pts = getGaussPts().size2();
901 auto t_strain = getFTensor2SymmetricFromMat<3, -1, CoeffsByGauss>(
902 *(dataAtPts->smallStrainMat));
903 dataAtPts->cauchyStressMat->resize(6, nb_integration_pts, false);
904 auto t_cauchy_stress =
905 getFTensor2SymmetricFromMat<3, -1, CoeffsByGauss>(
906 *(dataAtPts->cauchyStressMat));
907
908 FTensor::Index<'i', 3> i;
909 FTensor::Index<'j', 3> j;
910 FTensor::Index<'k', 3> k;
911 FTensor::Index<'l', 3> l;
912
913 // elastic stiffness tensor (4th rank tensor with minor and major
914 // symmetry)
915 FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
916 MAT_TO_DDG(dataAtPts->stiffnessMat));
917 for (int gg = 0; gg != nb_integration_pts; ++gg) {
918 t_cauchy_stress(i, j) = t_D(i, j, k, l) * t_strain(k, l);
919 ++t_strain;
920 ++t_cauchy_stress;
921 ++t_D;
922 }
923 MoFEMFunctionReturn(0);
924}
925
926template <int S>
927HookeElement::OpLhs_dx_dx<S>::OpLhs_dx_dx(
928 const std::string row_field, const std::string col_field,
929 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
930 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
931
932template <int S>
933MoFEMErrorCode HookeElement::OpLhs_dx_dx<S>::iNtegrate(EntData &row_data,
934 EntData &col_data) {
935 MoFEMFunctionBegin;
936
937 // get sub-block (3x3) of local stiffens matrix, here represented by
938 // second order tensor
939 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
940 return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
941 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
942 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
943 &m(r + 2, c + 2));
944 };
945
946 FTensor::Index<'i', 3> i;
947 FTensor::Index<'j', 3> j;
948 FTensor::Index<'k', 3> k;
949 FTensor::Index<'l', 3> l;
950
951 // get element volume
952 double vol = getVolume();
953
954 // get intergrayion weights
955 auto t_w = getFTensor0IntegrationWeight();
956
957 // get derivatives of base functions on rows
958 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
959 const int row_nb_base_fun = row_data.getN().size2();
960
961 // Elastic stiffness tensor (4th rank tensor with minor and major
962 // symmetry)
963 FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
964 MAT_TO_DDG(dataAtPts->stiffnessMat));
965
966 // iterate over integration points
967 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
968
969 // calculate scalar weight times element volume
970 double a = t_w * vol;
971
972 // iterate over row base functions
973 int rr = 0;
974 for (; rr != nbRows / 3; ++rr) {
975
976 // get sub matrix for the row
977 auto t_m = get_tensor2(K, 3 * rr, 0);
978
979 // get derivatives of base functions for columns
980 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
981
982 FTensor::Christof<double, 3, 3> t_rowD;
983 // I mix up the indices here so that it behaves like a
984 // Dg. That way I don't have to have a separate wrapper
985 // class Christof_Expr, which simplifies things.
986 t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
987
988 // iterate column base functions
989 for (int cc = 0; cc != nbCols / 3; ++cc) {
990
991 // integrate block local stiffens matrix
992 t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
993
994 // move to next column base function
995 ++t_col_diff_base;
996
997 // move to next block of local stiffens matrix
998 ++t_m;
999 }
1000
1001 // move to next row base function
1002 ++t_row_diff_base;
1003 }
1004
1005 for (; rr != row_nb_base_fun; ++rr)
1006 ++t_row_diff_base;
1007
1008 // move to next integration weight
1009 ++t_w;
1010 ++t_D;
1011 }
1012
1013 MoFEMFunctionReturn(0);
1014}
1015
1016template <int S>
1017HookeElement::OpAleLhs_dx_dX<S>::OpAleLhs_dx_dX(
1018 const std::string row_field, const std::string col_field,
1019 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1020 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
1021
1022template <int S>
1023MoFEMErrorCode HookeElement::OpAleLhs_dx_dX<S>::iNtegrate(EntData &row_data,
1024 EntData &col_data) {
1025 MoFEMFunctionBegin;
1026
1027 // get sub-block (3x3) of local stiffens matrix, here represented by
1028 // second order tensor
1029 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1030 return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
1031 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1032 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1033 &m(r + 2, c + 2));
1034 };
1035
1036 FTensor::Index<'i', 3> i;
1037 FTensor::Index<'j', 3> j;
1038 FTensor::Index<'k', 3> k;
1039 FTensor::Index<'l', 3> l;
1040 FTensor::Index<'m', 3> m;
1041 FTensor::Index<'n', 3> n;
1042
1043 // get element volume
1044 double vol = getVolume();
1045
1046 // get intergrayion weights
1047 auto t_w = getFTensor0IntegrationWeight();
1048
1049 // get derivatives of base functions on rows
1050 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1051 const int row_nb_base_fun = row_data.getN().size2();
1052
1053 // Elastic stiffness tensor (4th rank tensor with minor and major
1054 // symmetry)
1055 FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
1056 MAT_TO_DDG(dataAtPts->stiffnessMat));
1057
1058 auto t_cauchy_stress =
1059 getFTensor2SymmetricFromMat<3, -1, CoeffsByGauss>(
1060 *(dataAtPts->cauchyStressMat));
1061 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1062 auto t_invH =
1063 getFTensor2FromMat<3, 3, -1, CoeffsByGauss>(*dataAtPts->invHMat);
1064 auto &det_H = *dataAtPts->detHVec;
1065
1066 // iterate over integration points
1067 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1068
1069 // calculate scalar weight times element volume
1070 double a = t_w * vol * det_H[gg];
1071
1072 FTensor::Tensor4<double, 3, 3, 3, 3> t_F_dX;
1073 t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1074
1075 // iterate over row base functions
1076 int rr = 0;
1077 for (; rr != nbRows / 3; ++rr) {
1078
1079 // get sub matrix for the row
1080 auto t_m = get_tensor2(K, 3 * rr, 0);
1081
1082 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1083 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1084
1085 FTensor::Tensor1<double, 3> t_row_stress;
1086 t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_cauchy_stress(i, j);
1087
1088 FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
1089 t_row_diff_base_pulled_dX(j, k, l) =
1090 -(t_invH(i, k) * t_row_diff_base(i)) * t_invH(l, j);
1091
1092 FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
1093 t_row_dX_stress(i, k, l) =
1094 a * (t_row_diff_base_pulled_dX(j, k, l) * t_cauchy_stress(j, i));
1095
1096 FTensor::Christof<double, 3, 3> t_row_D;
1097 t_row_D(l, j, k) = (a * t_row_diff_base_pulled(i)) * t_D(i, j, k, l);
1098
1099 FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1100 // FIXME: This operator is not implemented, doing operation by hand
1101 // t_row_stress_dX(i, m, n) = t_row_D(i, k, l) * t_F_dX(k, l, m, n);
1102 t_row_stress_dX(i, j, k) = 0;
1103 for (int ii = 0; ii != 3; ++ii)
1104 for (int mm = 0; mm != 3; ++mm)
1105 for (int nn = 0; nn != 3; ++nn) {
1106 auto &v = t_row_stress_dX(ii, mm, nn);
1107 for (int kk = 0; kk != 3; ++kk)
1108 for (int ll = 0; ll != 3; ++ll)
1109 v += t_row_D(ii, kk, ll) * t_F_dX(kk, ll, mm, nn);
1110 }
1111
1112 // get derivatives of base functions for columns
1113 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1114
1115 // iterate column base functions
1116 for (int cc = 0; cc != nbCols / 3; ++cc) {
1117
1118 t_m(i, k) += t_row_stress(i) * (t_invH(j, k) * t_col_diff_base(j));
1119 t_m(i, k) += t_row_dX_stress(i, k, l) * t_col_diff_base(l);
1120 t_m(i, k) += t_row_stress_dX(i, k, l) * t_col_diff_base(l);
1121
1122 // move to next column base function
1123 ++t_col_diff_base;
1124
1125 // move to next block of local stiffens matrix
1126 ++t_m;
1127 }
1128
1129 // move to next row base function
1130 ++t_row_diff_base;
1131 }
1132
1133 for (; rr != row_nb_base_fun; ++rr)
1134 ++t_row_diff_base;
1135
1136 // move to next integration weight
1137 ++t_w;
1138 ++t_D;
1139 ++t_cauchy_stress;
1140 ++t_invH;
1141 ++t_h;
1142 }
1143
1144 MoFEMFunctionReturn(0);
1145}
1146
1147template <int S>
1148HookeElement::OpAleLhs_dX_dX<S>::OpAleLhs_dX_dX(
1149 const std::string row_field, const std::string col_field,
1150 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1151 : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
1152
1153template <int S>
1154MoFEMErrorCode HookeElement::OpAleLhs_dX_dX<S>::iNtegrate(EntData &row_data,
1155 EntData &col_data) {
1156 MoFEMFunctionBegin;
1157
1158 // get sub-block (3x3) of local stiffens matrix, here represented by
1159 // second order tensor
1160 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1161 return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
1162 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1163 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1164 &m(r + 2, c + 2));
1165 };
1166
1167 FTensor::Index<'i', 3> i;
1168 FTensor::Index<'j', 3> j;
1169 FTensor::Index<'k', 3> k;
1170 FTensor::Index<'l', 3> l;
1171 FTensor::Index<'m', 3> m;
1172 FTensor::Index<'n', 3> n;
1173
1174 // get element volume
1175 double vol = getVolume();
1176
1177 // get intergrayion weights
1178 auto t_w = getFTensor0IntegrationWeight();
1179
1180 // get derivatives of base functions on rows
1181 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1182 const int row_nb_base_fun = row_data.getN().size2();
1183
1184 // Elastic stiffness tensor (4th rank tensor with minor and major
1185 // symmetry)
1186 FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
1187 MAT_TO_DDG(dataAtPts->stiffnessMat));
1188 auto t_cauchy_stress =
1189 getFTensor2SymmetricFromMat<3, -1, CoeffsByGauss>(
1190 *(dataAtPts->cauchyStressMat));
1191 auto t_strain = getFTensor2SymmetricFromMat<3, -1, CoeffsByGauss>(
1192 *(dataAtPts->smallStrainMat));
1193 auto t_eshelby_stress =
1194 getFTensor2FromMat<3, 3, -1, CoeffsByGauss>(
1195 *dataAtPts->eshelbyStressMat);
1196 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1197 auto t_invH =
1198 getFTensor2FromMat<3, 3, -1, CoeffsByGauss>(*dataAtPts->invHMat);
1199 auto t_F = getFTensor2FromMat<3, 3, -1, CoeffsByGauss>(*dataAtPts->FMat);
1200 auto &det_H = *dataAtPts->detHVec;
1201
1202 // iterate over integration points
1203 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1204
1205 // calculate scalar weight times element volume
1206 double a = t_w * vol * det_H[gg];
1207
1208 FTensor::Tensor4<double, 3, 3, 3, 3> t_F_dX;
1209 t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1210
1211 FTensor::Tensor4<double, 3, 3, 3, 3> t_D_strain_dX;
1212 t_D_strain_dX(i, j, m, n) = 0.;
1213 for (int ii = 0; ii != 3; ++ii)
1214 for (int jj = 0; jj != 3; ++jj)
1215 for (int ll = 0; ll != 3; ++ll)
1216 for (int kk = 0; kk != 3; ++kk) {
1217 auto &v = t_D_strain_dX(ii, jj, kk, ll);
1218 for (int mm = 0; mm != 3; ++mm)
1219 for (int nn = 0; nn != 3; ++nn)
1220 v += t_D(ii, jj, mm, nn) * t_F_dX(mm, nn, kk, ll);
1221 }
1222
1223 FTensor::Tensor4<double, 3, 3, 3, 3> t_eshelby_stress_dX;
1224 t_eshelby_stress_dX(i, j, m, n) = t_F(k, i) * t_D_strain_dX(k, j, m, n);
1225
1226 for (int ii = 0; ii != 3; ++ii)
1227 for (int jj = 0; jj != 3; ++jj)
1228 for (int mm = 0; mm != 3; ++mm)
1229 for (int nn = 0; nn != 3; ++nn) {
1230 auto &v = t_eshelby_stress_dX(ii, jj, mm, nn);
1231 for (int kk = 0; kk != 3; ++kk)
1232 v += t_F_dX(kk, ii, mm, nn) * t_cauchy_stress(kk, jj);
1233 }
1234
1235 t_eshelby_stress_dX(i, j, k, l) *= -1;
1236
1237 FTensor::Tensor2<double, 3, 3> t_energy_dX;
1238 t_energy_dX(k, l) = t_F_dX(i, j, k, l) * t_cauchy_stress(i, j);
1239 t_energy_dX(k, l) +=
1240 (t_strain(m, n) * t_D(m, n, i, j)) * t_F_dX(i, j, k, l);
1241 t_energy_dX(k, l) /= 2.;
1242
1243 for (int kk = 0; kk != 3; ++kk)
1244 for (int ll = 0; ll != 3; ++ll) {
1245 auto v = t_energy_dX(kk, ll);
1246 for (int ii = 0; ii != 3; ++ii)
1247 t_eshelby_stress_dX(ii, ii, kk, ll) += v;
1248 }
1249
1250 // iterate over row base functions
1251 int rr = 0;
1252 for (; rr != nbRows / 3; ++rr) {
1253
1254 // get sub matrix for the row
1255 auto t_m = get_tensor2(K, 3 * rr, 0);
1256
1257 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1258 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1259
1260 FTensor::Tensor1<double, 3> t_row_stress;
1261 t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1262
1263 FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
1264 t_row_diff_base_pulled_dX(j, k, l) =
1265 -(t_row_diff_base(i) * t_invH(i, k)) * t_invH(l, j);
1266
1267 FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
1268 t_row_dX_stress(i, k, l) =
1269 a * (t_row_diff_base_pulled_dX(j, k, l) * t_eshelby_stress(i, j));
1270
1271 FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1272 t_row_stress_dX(i, m, n) =
1273 a * t_row_diff_base_pulled(j) * t_eshelby_stress_dX(i, j, m, n);
1274
1275 // get derivatives of base functions for columns
1276 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1277
1278 // iterate column base functions
1279 for (int cc = 0; cc != nbCols / 3; ++cc) {
1280
1281 t_m(i, k) += t_row_stress(i) * (t_invH(j, k) * t_col_diff_base(j));
1282 t_m(i, k) += t_row_dX_stress(i, k, l) * t_col_diff_base(l);
1283 t_m(i, k) += t_row_stress_dX(i, k, l) * t_col_diff_base(l);
1284
1285 // move to next column base function
1286 ++t_col_diff_base;
1287
1288 // move to next block of local stiffens matrix
1289 ++t_m;
1290 }
1291
1292 // move to next row base function
1293 ++t_row_diff_base;
1294 }
1295
1296 for (; rr != row_nb_base_fun; ++rr)
1297 ++t_row_diff_base;
1298
1299 // move to next integration weight
1300 ++t_w;
1301 ++t_D;
1302 ++t_cauchy_stress;
1303 ++t_strain;
1304 ++t_eshelby_stress;
1305 ++t_h;
1306 ++t_invH;
1307 ++t_F;
1308 }
1309
1310 MoFEMFunctionReturn(0);
1311}
1312
1313template <int S>
1314HookeElement::OpAleLhsPre_dX_dx<S>::OpAleLhsPre_dX_dx(
1315 const std::string row_field, const std::string col_field,
1316 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1317 : VolUserDataOperator(row_field, col_field, OPROW, true),
1318 dataAtPts(data_at_pts) {
1319 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1320}
1321
1322template <int S>
1323MoFEMErrorCode HookeElement::OpAleLhsPre_dX_dx<S>::doWork(int row_side,
1324 EntityType row_type,
1325 EntData &row_data) {
1326 MoFEMFunctionBegin;
1327
1328 const int nb_integration_pts = row_data.getN().size1();
1329
1330 auto get_eshelby_stress_dx = [this, nb_integration_pts]() {
1331 FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>
1332 t_eshelby_stress_dx;
1333 dataAtPts->eshelbyStress_dx->resize(81, nb_integration_pts, false);
1334 int mm = 0;
1335 for (int ii = 0; ii != 3; ++ii)
1336 for (int jj = 0; jj != 3; ++jj)
1337 for (int kk = 0; kk != 3; ++kk)
1338 for (int ll = 0; ll != 3; ++ll)
1339 t_eshelby_stress_dx.ptr(ii, jj, kk, ll) =
1340 &(*dataAtPts->eshelbyStress_dx)(mm++, 0);
1341 return t_eshelby_stress_dx;
1342 };
1343
1344 auto t_eshelby_stress_dx = get_eshelby_stress_dx();
1345
1346 FTensor::Index<'i', 3> i;
1347 FTensor::Index<'j', 3> j;
1348 FTensor::Index<'k', 3> k;
1349 FTensor::Index<'l', 3> l;
1350 FTensor::Index<'m', 3> m;
1351 FTensor::Index<'n', 3> n;
1352
1353 // Elastic stiffness tensor (4th rank tensor with minor and major
1354 // symmetry)
1355 FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
1356 MAT_TO_DDG(dataAtPts->stiffnessMat));
1357 auto t_cauchy_stress =
1358 getFTensor2SymmetricFromMat<3, -1, CoeffsByGauss>(
1359 *(dataAtPts->cauchyStressMat));
1360 auto t_invH =
1361 getFTensor2FromMat<3, 3, -1, CoeffsByGauss>(*dataAtPts->invHMat);
1362 auto t_F = getFTensor2FromMat<3, 3, -1, CoeffsByGauss>(*dataAtPts->FMat);
1363
1364 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1365
1366 t_eshelby_stress_dx(i, j, m, n) =
1367 (t_F(k, i) * t_D(k, j, m, l)) * t_invH(n, l);
1368 for (int ii = 0; ii != 3; ++ii)
1369 for (int jj = 0; jj != 3; ++jj)
1370 for (int mm = 0; mm != 3; ++mm)
1371 for (int nn = 0; nn != 3; ++nn) {
1372 auto &v = t_eshelby_stress_dx(ii, jj, mm, nn);
1373 v += t_invH(nn, ii) * t_cauchy_stress(mm, jj);
1374 }
1375 t_eshelby_stress_dx(i, j, k, l) *= -1;
1376
1377 FTensor::Tensor2<double, 3, 3> t_energy_dx;
1378 t_energy_dx(m, n) = t_invH(n, j) * t_cauchy_stress(m, j);
1379
1380 for (int mm = 0; mm != 3; ++mm)
1381 for (int nn = 0; nn != 3; ++nn) {
1382 auto v = t_energy_dx(mm, nn);
1383 for (int ii = 0; ii != 3; ++ii)
1384 t_eshelby_stress_dx(ii, ii, mm, nn) += v;
1385 }
1386
1387 ++t_D;
1388 ++t_invH;
1389 ++t_cauchy_stress;
1390 ++t_eshelby_stress_dx;
1391 ++t_F;
1392 }
1393
1394 MoFEMFunctionReturn(0);
1395}
1396
1397template <class ELEMENT>
1398HookeElement::OpPostProcHookeElement<ELEMENT>::OpPostProcHookeElement(
1399 const string row_field, boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1400 map<int, BlockData> &block_sets_ptr, moab::Interface &post_proc_mesh,
1401 std::vector<EntityHandle> &map_gauss_pts, bool is_ale, bool is_field_disp)
1402 : ELEMENT::UserDataOperator(row_field, UserDataOperator::OPROW),
1403 dataAtPts(data_at_pts), blockSetsPtr(block_sets_ptr),
1404 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), isALE(is_ale),
1405 isFieldDisp(is_field_disp) {}
1406
1407template <class ELEMENT>
1408MoFEMErrorCode HookeElement::OpPostProcHookeElement<ELEMENT>::doWork(
1409 int side, EntityType type, EntitiesFieldData::EntData &data) {
1410 MoFEMFunctionBegin;
1411
1412 if (type != MBVERTEX) {
1413 MoFEMFunctionReturnHot(0);
1414 }
1415
1416 auto tensor_to_tensor = [](const auto &t1, auto &t2) {
1417 t2(0, 0) = t1(0, 0);
1418 t2(1, 1) = t1(1, 1);
1419 t2(2, 2) = t1(2, 2);
1420 t2(0, 1) = t2(1, 0) = t1(1, 0);
1421 t2(0, 2) = t2(2, 0) = t1(2, 0);
1422 t2(1, 2) = t2(2, 1) = t1(2, 1);
1423 };
1424
1425 std::array<double, 9> def_val;
1426 def_val.fill(0);
1427
1428 auto make_tag = [&](auto name, auto size) {
1429 Tag th;
1430 CHKERR postProcMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
1431 MB_TAG_CREAT | MB_TAG_SPARSE,
1432 def_val.data());
1433 return th;
1434 };
1435
1436 auto th_stress = make_tag("STRESS", 9);
1437 auto th_psi = make_tag("ENERGY", 1);
1438
1439 const int nb_integration_pts = mapGaussPts.size();
1440
1441 FTensor::Index<'i', 3> i;
1442 FTensor::Index<'j', 3> j;
1443 FTensor::Index<'k', 3> k;
1444 FTensor::Index<'l', 3> l;
1445
1446 auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1447 auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
1448
1449 dataAtPts->stiffnessMat->resize(36, 1, false);
1450 FTensor::Ddg<FTensor::PackPtr<double *, 1>, 3, 3> t_D(
1451 MAT_TO_DDG(dataAtPts->stiffnessMat));
1452
1453 EntityHandle ent = this->getFEEntityHandle();
1454 auto type = type_from_handle(ent);
1455 EntityHandle ent_3d = ent;
1456 if (type == MBTRI || type == MBQUAD) {
1457 Range ents;
1458 auto &m_field = this->getPtrFE()->mField;
1459 CHKERR m_field.get_moab().get_adjacencies(&ent, 1, 3, false, ents,
1460 moab::Interface::UNION);
1461#ifndef NDEBUG
1462 if (ents.empty())
1463 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1464 "Could not find a 3D element adjacent to a given face element");
1465#endif
1466 ent_3d = ents.front();
1467 }
1468
1469 bool found_block = false;
1470 int block_id = -1;
1471 for (auto &m : (blockSetsPtr)) {
1472 if (m.second.tEts.find(ent_3d) != m.second.tEts.end()) {
1473 const double young = m.second.E;
1474 const double poisson = m.second.PoissonRatio;
1475 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
1476 block_id = m.second.iD;
1477
1478 t_D(i, j, k, l) = 0.;
1479 t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
1480 t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
1481 0.5 * (1 - 2 * poisson);
1482 t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
1483 t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
1484 t_D(i, j, k, l) *= coefficient;
1485
1486 found_block = true;
1487 break;
1488 }
1489 }
1490 if (!found_block)
1491 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1492 "Element not found in any of material blocksets");
1493
1494 int def_val_int = 0;
1495 Tag tag_mat;
1496 CHKERR postProcMesh.tag_get_handle("MAT_ELASTIC", 1, MB_TYPE_INTEGER, tag_mat,
1497 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
1498 double detH = 0.;
1499 FTensor::Tensor2<double, 3, 3> t_invH;
1500 FTensor::Tensor2<double, 3, 3> t_F;
1501 FTensor::Tensor2<double, 3, 3> t_stress;
1502 FTensor::Tensor2<double, 3, 3> t_small_strain;
1503 FTensor::Tensor2_symmetric<double, 3> t_stress_symm;
1504 FTensor::Tensor2_symmetric<double, 3> t_small_strain_symm;
1505
1506 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1507
1508 if (isFieldDisp) {
1509 t_h(0, 0) += 1;
1510 t_h(1, 1) += 1;
1511 t_h(2, 2) += 1;
1512 }
1513
1514 if (!isALE) {
1515 t_small_strain_symm(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
1516 } else {
1517 CHKERR determinantTensor3by3(t_H, detH);
1518 CHKERR invertTensor3by3(t_H, detH, t_invH);
1519 t_F(i, j) = t_h(i, k) * t_invH(k, j);
1520 t_small_strain_symm(i, j) = (t_F(i, j) || t_F(j, i)) / 2.;
1521 ++t_H;
1522 }
1523
1524 t_small_strain_symm(0, 0) -= 1;
1525 t_small_strain_symm(1, 1) -= 1;
1526 t_small_strain_symm(2, 2) -= 1;
1527
1528 // symmetric tensors need improvement
1529 t_stress_symm(i, j) = t_D(i, j, k, l) * t_small_strain_symm(k, l);
1530 tensor_to_tensor(t_stress_symm, t_stress);
1531
1532 const double psi = 0.5 * t_stress_symm(i, j) * t_small_strain_symm(i, j);
1533
1534 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
1535 CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
1536 &t_stress(0, 0));
1537 CHKERR postProcMesh.tag_set_data(tag_mat, &mapGaussPts[gg], 1,
1538 &block_id);
1539
1540 ++t_h;
1541 }
1542
1543 MoFEMFunctionReturn(0);
1544}
1545
1546template <int S>
1547HookeElement::OpAnalyticalInternalStrain_dx<S>::OpAnalyticalInternalStrain_dx(
1548 const std::string row_field,
1549 boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1550 StrainFunction strain_fun)
1551 : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1552 strainFun(strain_fun) {}
1553
1554template <int S>
1555MoFEMErrorCode
1556HookeElement::OpAnalyticalInternalStrain_dx<S>::iNtegrate(EntData &row_data) {
1557 FTensor::Index<'i', 3> i;
1558 FTensor::Index<'j', 3> j;
1559 FTensor::Index<'k', 3> k;
1560 FTensor::Index<'l', 3> l;
1561 MoFEMFunctionBegin;
1562
1563 auto get_tensor1 = [](VectorDouble &v, const int r) {
1564 return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
1565 &v(r + 0), &v(r + 1), &v(r + 2));
1566 };
1567
1568 const int nb_integration_pts = getGaussPts().size2();
1569 auto t_coords = getFTensor1CoordsAtGaussPts();
1570
1571 // get element volume
1572 double vol = getVolume();
1573 auto t_w = getFTensor0IntegrationWeight();
1574
1575 nF.resize(nbRows, false);
1576 nF.clear();
1577
1578 // elastic stiffness tensor (4th rank tensor with minor and major
1579 // symmetry)
1580 FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
1581 MAT_TO_DDG(dataAtPts->stiffnessMat));
1582
1583 // get derivatives of base functions on rows
1584 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1585 const int row_nb_base_fun = row_data.getN().size2();
1586
1587 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1588
1589 auto t_fun_strain = strainFun(t_coords);
1590 FTensor::Tensor2_symmetric<double, 3> t_stress;
1591 t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1592
1593 // calculate scalar weight times element volume
1594 double a = t_w * vol;
1595
1596 auto t_nf = get_tensor1(nF, 0);
1597
1598 int rr = 0;
1599 for (; rr != nbRows / 3; ++rr) {
1600 t_nf(i) += a * t_row_diff_base(j) * t_stress(i, j);
1601 ++t_row_diff_base;
1602 ++t_nf;
1603 }
1604
1605 for (; rr != row_nb_base_fun; ++rr)
1606 ++t_row_diff_base;
1607
1608 ++t_w;
1609 ++t_coords;
1610 ++t_D;
1611 }
1612
1613 MoFEMFunctionReturn(0);
1614}
1615
1616template <int S>
1617HookeElement::OpAnalyticalInternalAleStrain_dX<S>::
1618 OpAnalyticalInternalAleStrain_dX(
1619 const std::string row_field,
1620 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1621 StrainFunction strain_fun,
1622 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1623 : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1624 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1625
1626template <int S>
1627MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dX<S>::iNtegrate(
1628 EntData &row_data) {
1629 FTensor::Index<'i', 3> i;
1630 FTensor::Index<'j', 3> j;
1631 FTensor::Index<'k', 3> k;
1632 FTensor::Index<'l', 3> l;
1633 MoFEMFunctionBegin;
1634
1635 auto get_tensor1 = [](VectorDouble &v, const int r) {
1636 return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
1637 &v(r + 0), &v(r + 1), &v(r + 2));
1638 };
1639
1640 const int nb_integration_pts = getGaussPts().size2();
1641
1642 auto get_coords = [&]() {
1643 return getFTensor1FromMat<3>(*matPosAtPtsPtr);
1644 };
1645 auto t_coords = get_coords();
1646
1647 // get element volume
1648 double vol = getVolume();
1649 auto t_w = getFTensor0IntegrationWeight();
1650
1651 nF.resize(nbRows, false);
1652 nF.clear();
1653
1654 // elastic stiffness tensor (4th rank tensor with minor and major
1655 // symmetry)
1656 FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
1657 MAT_TO_DDG(dataAtPts->stiffnessMat));
1658 auto t_F = getFTensor2FromMat<3, 3, -1, CoeffsByGauss>(
1659 *(dataAtPts->FMat));
1660 auto &det_H = *dataAtPts->detHVec;
1661 auto t_invH =
1662 getFTensor2FromMat<3, 3, -1, CoeffsByGauss>(*dataAtPts->invHMat);
1663
1664 // get derivatives of base functions on rows
1665 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1666 const int row_nb_base_fun = row_data.getN().size2();
1667
1668 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1669
1670 auto t_fun_strain = strainFun(t_coords);
1671 FTensor::Tensor2_symmetric<double, 3> t_stress;
1672 t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1673 FTensor::Tensor2<double, 3, 3> t_eshelby_stress;
1674 t_eshelby_stress(i, j) = -t_F(k, i) * t_stress(k, j);
1675
1676 // calculate scalar weight times element volume
1677 double a = t_w * vol * det_H[gg];
1678
1679 auto t_nf = get_tensor1(nF, 0);
1680
1681 int rr = 0;
1682 for (; rr != nbRows / 3; ++rr) {
1683 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1684 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1685 t_nf(i) += a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1686 ++t_row_diff_base;
1687 ++t_nf;
1688 }
1689
1690 for (; rr != row_nb_base_fun; ++rr)
1691 ++t_row_diff_base;
1692
1693 ++t_w;
1694 ++t_coords;
1695 ++t_F;
1696 ++t_invH;
1697 ++t_D;
1698 }
1699
1700 MoFEMFunctionReturn(0);
1701}
1702
1703template <int S>
1704HookeElement::OpAnalyticalInternalAleStrain_dx<S>::
1705 OpAnalyticalInternalAleStrain_dx(
1706 const std::string row_field,
1707 boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1708 StrainFunction strain_fun,
1709 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1710 : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1711 strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1712
1713template <int S>
1714MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dx<S>::iNtegrate(
1715 EntData &row_data) {
1716 FTensor::Index<'i', 3> i;
1717 FTensor::Index<'j', 3> j;
1718 FTensor::Index<'k', 3> k;
1719 FTensor::Index<'l', 3> l;
1720 MoFEMFunctionBegin;
1721
1722 auto get_tensor1 = [](VectorDouble &v, const int r) {
1723 return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
1724 &v(r + 0), &v(r + 1), &v(r + 2));
1725 };
1726
1727 const int nb_integration_pts = getGaussPts().size2();
1728
1729 auto get_coords = [&]() {
1730 return getFTensor1FromMat<3>(*matPosAtPtsPtr);
1731 };
1732 auto t_coords = get_coords();
1733
1734 // get element volume
1735 double vol = getVolume();
1736 auto t_w = getFTensor0IntegrationWeight();
1737
1738 nF.resize(nbRows, false);
1739 nF.clear();
1740
1741 // elastic stiffness tensor (4th rank tensor with minor and major
1742 // symmetry)
1743 FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
1744 MAT_TO_DDG(dataAtPts->stiffnessMat));
1745 auto &det_H = *dataAtPts->detHVec;
1746 auto t_invH =
1747 getFTensor2FromMat<3, 3, -1, CoeffsByGauss>(*dataAtPts->invHMat);
1748
1749 // get derivatives of base functions on rows
1750 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1751 const int row_nb_base_fun = row_data.getN().size2();
1752
1753 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1754
1755 auto t_fun_strain = strainFun(t_coords);
1756 FTensor::Tensor2_symmetric<double, 3> t_stress;
1757 t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1758
1759 // calculate scalar weight times element volume
1760 double a = t_w * vol * det_H[gg];
1761
1762 auto t_nf = get_tensor1(nF, 0);
1763
1764 int rr = 0;
1765 for (; rr != nbRows / 3; ++rr) {
1766 FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1767 t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1768 t_nf(i) += a * t_row_diff_base_pulled(j) * t_stress(i, j);
1769 ++t_row_diff_base;
1770 ++t_nf;
1771 }
1772
1773 for (; rr != row_nb_base_fun; ++rr)
1774 ++t_row_diff_base;
1775
1776 ++t_w;
1777 ++t_coords;
1778 ++t_invH;
1779 ++t_D;
1780 }
1781
1782 MoFEMFunctionReturn(0);
1783}
1784
1785#endif // __HOOKE_ELEMENT_HPP
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)
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)
MAT_TO_DDG
#define MAT_TO_DDG(SM)
Definition HookeElement.hpp:142
setBlocks
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr)
invJac
MatrixDouble invJac
Definition HookeElement.hpp:683
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)
type
std::string type
Definition JsonConfigManager.cpp:28
a
constexpr double a
Definition approx_sphere.cpp:30
FTensor::Tensor2_symmetric
Definition Tensor2_symmetric_value.hpp:14
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
c
const double c
speed of light (cm/ns)
Definition initial_diffusion.cpp:40
D
double D
Definition initial_diffusion.cpp:45
v
const double v
phase velocity of light in medium (cm/ns)
Definition initial_diffusion.cpp:41
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:39
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:23
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:22
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM::Types::VectorInt
UBlasVector< int > VectorInt
Definition Types.hpp:67
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition Templates.hpp:2548
MoFEM::getFTensor2FromMat
auto getFTensor2FromMat(M &data)
Get tensor rank 2 (matrix) form data matrix.
Definition Templates.hpp:542
MoFEM::getFTensor2SymmetricFromMat
static auto getFTensor2SymmetricFromMat(M &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
Definition Templates.hpp:827
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:2276
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:1929
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition Templates.hpp:2228
sdf_wavy_2d.w
int w
Definition sdf_wavy_2d.py:6
sdf.r
int r
Definition sdf.py:205
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
ConvectiveMassElement::BlockData
data for calculation inertia forces
Definition ConvectiveMassElement.hpp:116
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition EntitiesFieldData.hpp:150
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition EntitiesFieldData.hpp:1806
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
auto getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition EntitiesFieldData.hpp:802
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:1600
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
Definition EntitiesFieldData.hpp:1523
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROWCOL
@ OPROWCOL
operator doWork is executed on FE rows &columns
Definition ForcesAndSourcesCore.hpp:604
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition PetscSmartObj.hpp:85
OpAleLhsPre_dX_dx
Definition HookeElement.hpp:536
OpAleLhsPre_dX_dx::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition HookeElement.hpp:544
OpAleLhsPre_dX_dx::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
Definition HookeElement.hpp:1324
OpAleLhsWithDensity_dX_dX
Definition HookeElement.hpp:488
OpAleLhsWithDensity_dX_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for material momentum.
Definition HookeElement.cpp:910
OpAleLhsWithDensity_dX_dX::rhoN
const double rhoN
Definition HookeElement.hpp:492
OpAleLhsWithDensity_dX_dX::rhoAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
Definition HookeElement.hpp:490
OpAleLhsWithDensity_dX_dX::rHo0
const double rHo0
Definition HookeElement.hpp:493
OpAleLhsWithDensity_dX_dX::rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
Definition HookeElement.hpp:491
OpAleLhsWithDensity_dx_dX
Definition HookeElement.hpp:464
OpAleLhsWithDensity_dx_dX::rHo0
const double rHo0
Definition HookeElement.hpp:469
OpAleLhsWithDensity_dx_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for spatial momentum.
Definition HookeElement.cpp:1020
OpAleLhsWithDensity_dx_dX::rhoAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
Definition HookeElement.hpp:466
OpAleLhsWithDensity_dx_dX::rhoN
const double rhoN
Definition HookeElement.hpp:468
OpAleLhsWithDensity_dx_dX::rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
Definition HookeElement.hpp:467
OpAleLhs_dX_dX
Definition HookeElement.hpp:521
OpAleLhs_dX_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for material momentum.
Definition HookeElement.hpp:1155
OpAleLhs_dX_dx
Definition HookeElement.hpp:547
OpAleLhs_dX_dx::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for material momentum.
Definition HookeElement.cpp:800
OpAleLhs_dx_dX
Definition HookeElement.hpp:449
OpAleLhs_dx_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate tangent stiffness for spatial momentum.
Definition HookeElement.hpp:1024
OpAleLhs_dx_dx
Definition HookeElement.hpp:434
OpAleLhs_dx_dx::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate B^T D B operator.
Definition HookeElement.hpp:733
OpAleRhs_dX
Definition HookeElement.hpp:512
OpAleRhs_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data)
Definition HookeElement.cpp:466
OpAleRhs_dx
Definition HookeElement.hpp:425
OpAleRhs_dx::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data)
Definition HookeElement.cpp:408
OpAnalyticalInternalAleStrain_dX
Definition HookeElement.hpp:586
OpAnalyticalInternalAleStrain_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data)
Definition HookeElement.hpp:1628
OpAnalyticalInternalAleStrain_dX::matPosAtPtsPtr
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
Definition HookeElement.hpp:608
OpAnalyticalInternalAleStrain_dX::StrainFunction
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunction
Definition HookeElement.hpp:597
OpAnalyticalInternalAleStrain_dX::strainFun
StrainFunction strainFun
Definition HookeElement.hpp:607
OpAnalyticalInternalAleStrain_dx
Definition HookeElement.hpp:611
OpAnalyticalInternalAleStrain_dx::matPosAtPtsPtr
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
Definition HookeElement.hpp:633
OpAnalyticalInternalAleStrain_dx::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data)
Definition HookeElement.hpp:1715
OpAnalyticalInternalAleStrain_dx::StrainFunction
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunction
Definition HookeElement.hpp:622
OpAnalyticalInternalAleStrain_dx::strainFun
StrainFunction strainFun
Definition HookeElement.hpp:632
OpAnalyticalInternalStrain_dx
Definition HookeElement.hpp:563
OpAnalyticalInternalStrain_dx::StrainFunction
boost::function< FTensor::Tensor2_symmetric< double, 3 >(FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > &t_coords) > StrainFunction
Definition HookeElement.hpp:574
OpAnalyticalInternalStrain_dx::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data)
Definition HookeElement.hpp:1557
OpAnalyticalInternalStrain_dx::strainFun
StrainFunction strainFun
Definition HookeElement.hpp:583
OpAssemble::nbCols
int nbCols
number if dof on column
Definition HookeElement.hpp:374
OpAssemble::iNtegrate
virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Definition HookeElement.cpp:254
OpAssemble::isDiag
bool isDiag
true if this block is on diagonal
Definition HookeElement.hpp:376
OpAssemble::nF
VectorDouble nF
Definition HookeElement.hpp:366
OpAssemble::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Assemble local entity block matrix.
Definition HookeElement.cpp:267
OpAssemble::doWork
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.
Definition HookeElement.cpp:183
OpAssemble::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition HookeElement.hpp:368
OpAssemble::nbIntegrationPts
int nbIntegrationPts
number of integration points
Definition HookeElement.hpp:375
OpAssemble::nbRows
int nbRows
number of dofs on rows
Definition HookeElement.hpp:373
OpAssemble::transK
MatrixDouble transK
Definition HookeElement.hpp:365
OpAssemble::colIndices
VectorInt colIndices
Definition HookeElement.hpp:371
OpAssemble::rowIndices
VectorInt rowIndices
Definition HookeElement.hpp:370
OpAssemble::K
MatrixDouble K
Definition HookeElement.hpp:364
OpCalculateEnergy::ghostVec
SmartPetscObj< Vec > ghostVec
Definition HookeElement.hpp:174
OpCalculateEnergy::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
Definition HookeElement.cpp:83
OpCalculateEnergy::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition HookeElement.hpp:173
OpCalculateEshelbyStress::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
Definition HookeElement.cpp:142
OpCalculateEshelbyStress::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
Definition HookeElement.hpp:186
OpCalculateHomogeneousStiffness::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
Definition HookeElement.hpp:839
OpCalculateHomogeneousStiffness::blockSetsPtr
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
Definition HookeElement.hpp:201
OpCalculateHomogeneousStiffness::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition HookeElement.hpp:203
OpCalculateMassMatrix
Assemble mass matrix for elastic element TODO: CHANGE FORMULA *.
Definition HookeElement.hpp:213
OpCalculateMassMatrix::doWork
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition HookeElement.hpp:231
OpCalculateMassMatrix::massData
MassBlockData & massData
Definition HookeElement.hpp:218
OpCalculateMassMatrix::commonData
boost::shared_ptr< DataAtIntegrationPts > commonData
Definition HookeElement.hpp:220
OpCalculateStiffnessScaledByDensityField::rHo0
const double rHo0
p_0 reference density in E(p) = E * (p / p_0)^n
Definition HookeElement.hpp:327
OpCalculateStiffnessScaledByDensityField::blockSetsPtr
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
Definition HookeElement.hpp:321
OpCalculateStiffnessScaledByDensityField::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition HookeElement.hpp:323
OpCalculateStiffnessScaledByDensityField::rhoAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
Definition HookeElement.hpp:325
OpCalculateStiffnessScaledByDensityField::rhoN
const double rhoN
exponent n in E(p) = E * (p / p_0)^n
Definition HookeElement.hpp:326
OpCalculateStiffnessScaledByDensityField::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
Definition HookeElement.cpp:1128
OpLhs_dx_dx
Definition HookeElement.hpp:410
OpLhs_dx_dx::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate B^T D B operator.
Definition HookeElement.hpp:934
OpPostProcHookeElement::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition HookeElement.hpp:642
OpPostProcHookeElement::blockSetsPtr
map< int, BlockData > & blockSetsPtr
Definition HookeElement.hpp:640
OpPostProcHookeElement::postProcMesh
moab::Interface & postProcMesh
Definition HookeElement.hpp:641
OpPostProcHookeElement::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition HookeElement.hpp:1409
OpPostProcHookeElement::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition HookeElement.hpp:638
OpRhs_dx
Definition HookeElement.hpp:401
OpRhs_dx::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data)
Definition HookeElement.cpp:354

Typedef Documentation

◆ EntData

using EntData = EntitiesFieldData::EntData

Definition at line 74 of file HookeElement.hpp.

◆ MassBlockData

using MassBlockData = ConvectiveMassElement::BlockData

Definition at line 72 of file HookeElement.hpp.

◆ UserDataOperator

using UserDataOperator = ForcesAndSourcesCore::UserDataOperator

Definition at line 75 of file HookeElement.hpp.

◆ VolUserDataOperator

using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator

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
Examples
mofem/users_modules/basic_finite_elements/src/HookeElement.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 
)
static
Examples
mofem/users_modules/basic_finite_elements/src/HookeElement.hpp.

◆ setBlocks()

static MoFEMErrorCode setBlocks ( MoFEM::Interface &  m_field,
boost::shared_ptr< map< int, BlockData > > &  block_sets_ptr 
)
static
Examples
mofem/users_modules/basic_finite_elements/src/HookeElement.hpp.

◆ 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
Examples
mofem/users_modules/basic_finite_elements/src/HookeElement.hpp.

Variable Documentation

◆ invJac

MatrixDouble invJac
private
Examples
mofem/users_modules/basic_finite_elements/src/HookeElement.hpp.

Definition at line 683 of file HookeElement.hpp.

Generated by Doxygen 1.9.8 and hosted at University of Glasgow