v0.14.0
Loading...
Searching...
No Matches
SolidShellPrismElement.hpp
Go to the documentation of this file.
1/** \file SolidShellPrismElement.hpp
2 * \ingroup solid_shell_prism_element
3 * \brief Implementation of solid shell prism element
4 *
5 */
6
7/*
8 * This file is part of MoFEM.
9 * MoFEM is free software: you can redistribute it and/or modify it under
10 * the terms of the GNU Lesser General Public License as published by the
11 * Free Software Foundation, either version 3 of the License, or (at your
12 * option) any later version.
13 *
14 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
15 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 * License for more details.
18 *
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
21
22#ifndef __SOLIDSHELLPRISMELEMENT_HPP__
23#define __SOLIDSHELLPRISMELEMENT_HPP__
24
25namespace SolidShellModule {
26
27/** \brief Main class for solid shell element
28 \ingroup solid_shell_prism_element
29*/
31
32 /** \brief Common data for solid shell element
33 \ingroup solid_shell_prism_element
34 */
35 struct CommonData {
36
37 string meshNodalPosition; ///< Name of field of positions in reference
38 ///< configuration
39 string displacementFieldName; ///< Name of displacement field
40 string zetaFieldName; ///< Name of zeta field
41 string ksiEtaFieldName; ///< In plane HO approximation
44
45 vector<ublas::matrix<double, ublas::column_major>>
46 jacobianInvAtGaussPts; ///< For element local to global
47 MatrixDouble positionsAtGaussPts; ///< Reference gauss points coordinates
48 MatrixDouble directorVectorAtGaussPts; ///< Approximation of director vector
49 ///< at integration Pts
50 vector<MatrixDouble>
51 diffDirectorVectorAtGaussPts; ///< Director vector derivatives at
52 ///< integration pts
53 MatrixDouble tangent0AtGaussPts; ///< Approximation of tangent 0 vector at
54 ///< integration Pts
55 MatrixDouble tangent1AtGaussPts; ///< Approximation of tangent 1 vector at
56 ///< integration Pts
57 MatrixDouble
58 localDisplacementsAtGaussPts; ///< Local displacement at integration Pts
59 MatrixDouble
60 globalDisplacentAtGaussPts; ///< Global displacement at integration ptrs
61 VectorDouble jacDeterminantAtGaussPts; ///< Determinant of jacobian
62 map<EntityType, map<int, MatrixDouble>> diffGlobN;
63
64 MatrixDouble sTrain; ///< Strain in Voigt
65 ///< (https://en.wikipedia.org/wiki/Voigt_notation)
66 ///< notation (User coord. sys.)
67 MatrixDouble sTress; ///< Stress in Voigt
68 ///< (https://en.wikipedia.org/wiki/Voigt_notation)
69 ///< notation (User coors. sys.)
70 MatrixDouble
71 materialStiffness; ///< Material stiffness matrix in Voigt notation
72
73 double tHickness; ///< Shell thickness
74 double sHift; ///< Shift elements in normal direction
75 VectorDouble
76 surfaceParametrisation; ///< User direction of surface parametrization
77 // vector<MatrixDouble> rOtationUser; ///< Rotate from local user
78 // coord. sys. to global
79 vector<ublas::symmetric_matrix<double, ublas::lower>>
80 localStress; ///< Stress in local user coordinate system (Tensorial
81 ///< Notation)
82
83 vector<MatrixDouble> mGrad1,
84 mGrad2; ///< Mid-surface gradient of deformation
85
86 // Note: "_" before index indicate lower index, if no "_" it is upper index
87 vector<MatrixDouble3by3> mE_I_A, mEAI, mGAB, mG_A_B;
88 vector<MatrixDouble3by3> globalPiolaStressIi_I;
89 vector<MatrixDouble3by3> localPiolaStressIA_a;
90 vector<MatrixDouble3by3> globalPiolaStressIiBTangent;
91 vector<MatrixDouble3by3> localPiolaStressIaBTangent;
92
94 vector<MatrixDouble> tangentOfPiolaStressLocalLocal;
95 vector<MatrixDouble> tangentOfPiolaStressLocalGlobal;
96 vector<MatrixDouble> tangentOfPiolaStressGlobalLocal;
97
98 // CoVectors
99 map<EntityHandle, int> shellElementsMap;
100 map<int, map<UId, double>> coVector0DofsMap;
101 map<int, map<UId, double>> coVector1DofsMap;
102
103 // This is to evaluate error
106
107#ifdef ADOLC_ADOLC_H
108
109 template <typename TYPE> struct ADouble {
110 TYPE dEt, tR;
111 // Acrive Variables
112 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>>
113 mGrad1;
114 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>>
115 mGrad2;
116
117 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>> mMF,
118 mF;
119 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>> gMab,
120 gM_a_b;
121 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>>
122 eM_i_a, eMai;
123
124 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>>
125 mCA_B;
126 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>>
127 mEA_B;
128 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>> mCAB;
129 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>> mEAB;
130
131 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>>
132 piolaStressIIA_B;
133 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>>
134 globalPiolaStressIi_I;
135 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>>
136 localPiolaStressIA_a;
137 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>>
138 localPiolaStressIa_A;
139
140 TYPE t0, t1, t2, t3, t4;
141 };
142
143 ADouble<adouble> aDouble;
144 ADouble<double> dDouble;
145
146#endif // ADOLC_ADOLC_H
147
148#ifdef __ARC_LENGTH_TOOLS_HPP__
149 vector<VectorDouble3> incrementsOfDisplacements;
150 // vector<MatrixDouble3> incrementsOfGradOfDisplacements;
151
152#endif
153
155 };
156
158
160 SolidShell(MoFEM::Interface &m_field, CommonData &common_data)
162 commonData(common_data) {}
164 Tag th;
165
166 rval = mField.get_moab().tag_get_handle("ADAPT_ORDER", th);
167 if (rval == MB_SUCCESS) {
168 order = 2;
169 int last_ent = 0;
170 for (_IT_FENUMEREDDOF_BY_NAME_ROW_FOR_LOOP_(numeredEntFiniteElementPtr,
171 "DISPLACEMENT", dof)) {
172 EntityHandle ent = dof->get()->getEnt();
173 if (last_ent == ent)
174 continue;
175 last_ent = ent;
176 int ent_order;
177 rval = mField.get_moab().tag_get_data(th, &ent, 1, &ent_order);
178 MOAB_THROW(rval);
179 order = (order > ent_order) ? order : ent_order;
180 }
181 }
182 return 2 * (order + 1) + 1;
183 };
185 int rule = 2 * order + 1;
186 return (rule <= 1) ? 2 : rule;
187 }
188 PetscErrorCode preProcess();
189 PetscErrorCode postProcess();
190 };
191
192 struct SolidShellError : public SolidShell {
194 : SolidShell(m_field, common_data) {}
195
196 PetscErrorCode preProcess();
197 PetscErrorCode postProcess();
198 };
199
200 /** \brief Calculate strain element operators
201 \ingroup solid_shell_prism_element
202 */
203 struct MakeB {
204 /** \brief Calculate strain operator for displacements DoFs
205
206 Note: Displacements DoFs are in global coordinate system
207
208 */
209 PetscErrorCode makeBDisp(const MatrixAdaptor &diffN, const MatrixDouble &T,
210 MatrixDouble &B);
211
212 /** \brief Calculate strain operator for displacement thought shell
213 thickness
214
215 Note: Zeta displacements are in local coordinate system
216
217 */
218 PetscErrorCode makeBZeta(const MatrixAdaptor &diffN, MatrixDouble &B);
219
220 /** \brief Calculate strain operator for non-plannar displacements
221 through shell thickness
222
223 Note: Ksi and Eta displacements are in local coordinate system
224
225 */
226 PetscErrorCode makeBKsiEta(const MatrixAdaptor &diffN, MatrixDouble &B);
227
228 /** \brief Gradient of deformation
229
230 Values of tensor are arranged in columns (row order)
231
232 */
233 PetscErrorCode makeFDisp(const MatrixAdaptor &diffN, MatrixDouble &F);
234
235 PetscErrorCode makeFZeta(const MatrixAdaptor &diffN, MatrixDouble &F);
236
237 PetscErrorCode makeFKsiEta(const MatrixAdaptor &diffN, MatrixDouble &F);
238
239 template <typename TYPE>
240 PetscErrorCode dEterminatnt(
241 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>> &a,
242 TYPE &det) {
243 PetscFunctionBegin;
244 // a11a22a33
245 //+a21a32a13
246 //+a31a12a23
247 //-a11a32a23
248 //-a31a22a13
249 //-a21a12a33
250 // http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
251 // http://mathworld.wolfram.com/MatrixInverse.html
252 det = a(0, 0) * a(1, 1) * a(2, 2) + a(1, 0) * a(2, 1) * a(0, 2) +
253 a(2, 0) * a(0, 1) * a(1, 2) - a(0, 0) * a(2, 1) * a(1, 2) -
254 a(2, 0) * a(1, 1) * a(0, 2) - a(1, 0) * a(0, 1) * a(2, 2);
255 PetscFunctionReturn(0);
256 }
257
258 /** \brief Calculate inverse of 3x3 matrix
259 */
260 template <typename TYPE>
261 PetscErrorCode iNvert(
262 TYPE det,
263 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>> &a,
264 ublas::matrix<TYPE, ublas::row_major, ublas::bounded_array<TYPE, 9>>
265 &inv_a) {
266 PetscFunctionBegin;
267 //
268 inv_a.resize(3, 3, false);
269 // http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
270 // http://mathworld.wolfram.com/MatrixInverse.html
271 inv_a(0, 0) = a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1);
272 inv_a(0, 1) = a(0, 2) * a(2, 1) - a(0, 1) * a(2, 2);
273 inv_a(0, 2) = a(0, 1) * a(1, 2) - a(0, 2) * a(1, 1);
274 inv_a(1, 0) = a(1, 2) * a(2, 0) - a(1, 0) * a(2, 2);
275 inv_a(1, 1) = a(0, 0) * a(2, 2) - a(0, 2) * a(2, 0);
276 inv_a(1, 2) = a(0, 2) * a(1, 0) - a(0, 0) * a(1, 2);
277 inv_a(2, 0) = a(1, 0) * a(2, 1) - a(1, 1) * a(2, 0);
278 inv_a(2, 1) = a(0, 1) * a(2, 0) - a(0, 0) * a(2, 1);
279 inv_a(2, 2) = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
280 inv_a /= det;
281 PetscFunctionReturn(0);
282 }
283 };
284
285 /** \brief Calculate displacements at integration points
286 \ingroup solid_shell_prism_element
287 */
292 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
293 common_data.displacementFieldName, UserDataOperator::OPROW),
294 commonData(common_data) {}
295 PetscErrorCode doWork(int side, EntityType type,
296 DataForcesAndSourcesCore::EntData &data);
297 };
298
303 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
304 common_data.directorFieldName, UserDataOperator::OPROW),
305 commonData(common_data) {}
306 PetscErrorCode doWork(int side, EntityType type,
307 DataForcesAndSourcesCore::EntData &data);
308 };
309
313 VectorDouble coVector0Data;
314 VectorDouble coVector1Data;
316 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
317 common_data.coVectorFieldName, UserDataOperator::OPROW),
318 commonData(common_data) {}
319 PetscErrorCode doWork(int side, EntityType type,
320 DataForcesAndSourcesCore::EntData &data);
321 };
322
323 /** \brief Calculate rotation matrices to local (user) element coordinate
324 system \ingroup solid_shell_prism_element
325
326 Material physical equation is defined and evaluated in local element
327 coordinate system.
328
329 */
332 MakeB {
335 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
336 common_data.meshNodalPosition, UserDataOperator::OPROW),
337 commonData(common_data) {}
338 PetscErrorCode doWork(int side, EntityType type,
339 DataForcesAndSourcesCore::EntData &data);
340 };
341
342 /** \brief Calculate jacobian
343 \ingroup solid_shell_prism_element
344 */
349 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
350 common_data.meshNodalPosition, UserDataOperator::OPROW),
351 commonData(common_data) {}
352 PetscErrorCode doWork(int side, EntityType type,
353 DataForcesAndSourcesCore::EntData &data);
354 };
355
356 /** \brief Calculate inverse of Jacobian
357 \ingroup solid_shell_prism_element
358
359 In addition all derivatives of shape functions are transformed to
360 local element coordinate system
361
362 */
367 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
368 common_data.displacementFieldName, UserDataOperator::OPROW),
369 commonData(common_data) {}
370 PetscErrorCode doWork(int side, EntityType type,
371 DataForcesAndSourcesCore::EntData &data);
372 };
373
374 /** \brief Get strain form displacements at integration points
375 \ingroup solid_shell_prism_element
376 */
379 MakeB {
381 MatrixDouble mB;
382 VectorDouble sTrain;
384 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
385 common_data.displacementFieldName, UserDataOperator::OPROW),
386 commonData(common_data) {}
387 PetscErrorCode doWork(int side, EntityType type,
388 DataForcesAndSourcesCore::EntData &data);
389 };
390
391 /** \brief Get strain form local through thickness displacements at
392 integration points \ingroup solid_shell_prism_element
393 */
396 MakeB {
398 MatrixDouble mB;
399 VectorDouble sTrain;
400 OpGetStrainLocal(CommonData &common_data, const string &field)
401 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
402 field, UserDataOperator::OPROW),
403 commonData(common_data) {}
404 PetscErrorCode doWork(int side, EntityType type,
405 DataForcesAndSourcesCore::EntData &data);
406 };
407
408 /** \brief Evaluate stress at integration points
409 \ingroup solid_shell_prism_element
410
411 Stress is expressed in local coordinate system
412 */
417 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
418 common_data.displacementFieldName, UserDataOperator::OPROW),
419 commonData(common_data) {}
420 PetscErrorCode doWork(int side, EntityType type,
421 DataForcesAndSourcesCore::EntData &data);
422 };
423
424 /** \brief Assemble internal forces for displacement DoFs
425 \ingroup solid_shell_prism_element
426 */
429 MakeB {
431 VectorDouble nF;
432 MatrixDouble B;
433
435 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
436 common_data.displacementFieldName, UserDataOperator::OPROW),
437 commonData(common_data) {}
438 PetscErrorCode doWork(int side, EntityType type,
439 DataForcesAndSourcesCore::EntData &data);
440 };
441
442 /** \brief Assemble internal forces for local thorough thickness DoFs
443 \ingroup solid_shell_prism_element
444 */
447 MakeB {
449 VectorDouble nF;
450 MatrixDouble mB;
451
452 OpAssembleRhsLocal(CommonData &common_data, const string &field)
453 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
454 field, UserDataOperator::OPROW),
455 commonData(common_data) {}
456 PetscErrorCode doWork(int side, EntityType type,
457 DataForcesAndSourcesCore::EntData &data);
458 };
459
460 /** \brief Assemble stiffness matrix
461 */
464 MakeB {
466 MatrixDouble mK, mTransK;
467 MatrixDouble mRowB, mColB;
468 MatrixDouble mCB;
469
470 OpAssembleLhsMix(CommonData &common_data, const string &row_field_name,
471 const string &col_field_name, const bool symm = false)
472 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
473 row_field_name, col_field_name, UserDataOperator::OPROWCOL),
474 commonData(common_data) {
475 sYmm = symm;
476 }
477 PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
478 EntityType col_type,
479 DataForcesAndSourcesCore::EntData &row_data,
480 DataForcesAndSourcesCore::EntData &col_data);
481 };
482
483 /** \brief Post process strain and stresses on post-processing mesh
484 \ingroup solid_shell_prism_element
485 */
488 moab::Interface &postProcMesh;
489 map<EntityHandle, EntityHandle> &elementsMap;
490 vector<EntityHandle> &mapGaussPts;
492 OpPostProcStressAndStrain(moab::Interface &post_proc_mesh,
493 map<EntityHandle, EntityHandle> &elements_map,
494 vector<EntityHandle> &map_gauss_pts,
495 CommonData &common_data)
496 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
497 common_data.displacementFieldName, UserDataOperator::OPROW),
498 postProcMesh(post_proc_mesh), elementsMap(elements_map),
499 mapGaussPts(map_gauss_pts), commonData(common_data) {}
500 PetscErrorCode doWork(int side, EntityType type,
501 DataForcesAndSourcesCore::EntData &data);
502 };
503
504 /** \brief Post process displacements on post-processing mesh
505 \ingroup solid_shell_prism_element
506 */
509 moab::Interface &postProcMesh;
510 vector<EntityHandle> &mapGaussPts;
512 OpPostProcMeshNodeDispalcements(moab::Interface &post_proc_mesh,
513 vector<EntityHandle> &map_gauss_pts,
514 CommonData &common_data)
515 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
516 common_data.displacementFieldName, UserDataOperator::OPROW),
517 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
518 commonData(common_data) {}
519 MatrixDouble meshNodePositions;
520 PetscErrorCode doWork(int side, EntityType type,
521 DataForcesAndSourcesCore::EntData &data);
522 };
523
524 /** \brief Post process geometry (mesh nodal positions) on post-processing
525 mesh \ingroup solid_shell_prism_element
526 */
529 moab::Interface &postProcMesh;
530 vector<EntityHandle> &mapGaussPts;
532 OpPostProcMeshNodePositions(moab::Interface &post_proc_mesh,
533 vector<EntityHandle> &map_gauss_pts,
534 CommonData &common_data)
535 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
536 common_data.meshNodalPosition, UserDataOperator::OPROW),
537 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
538 commonData(common_data) {}
539 MatrixDouble meshNodePositions;
540 PetscErrorCode doWork(int side, EntityType type,
541 DataForcesAndSourcesCore::EntData &data);
542 };
543
544#ifndef __POSTPROC_ON_REF_MESH_HPP
545#error "Include file PostProcOnRefMesh.hpp before this header"
546#endif //__POSTPROC_ON_REF_MESH_HPP
547
548 /** \brief Generate post-processing triangular mesh for solid-shell element
549 \ingroup solid_shell_prism_element
550 */
554 MoFEM::Interface &m_field, bool ten_nodes_post_proc_tets = true)
555 : PostProcFatPrismOnRefinedMesh(m_field, ten_nodes_post_proc_tets) {}
556 int getRuleTrianglesOnly(int order) { return -1; };
557 int getRuleThroughThickness(int order) { return -1; };
558 PetscErrorCode setGaussPtsTrianglesOnly(int order_triangles_only);
559 PetscErrorCode setGaussPtsThroughThickness(int order_thickness);
560 PetscErrorCode generateReferenceElementMesh();
561 PetscErrorCode postProcess();
562 };
563
564 /** \brief Post process axial and shear forces
565 \ingroup solid_shell_prism_element
566 */
569 moab::Interface &postProcMesh;
570 vector<EntityHandle> &mapGaussPts;
572 OpGetAxialForces(moab::Interface &post_proc_mesh,
573 vector<EntityHandle> &map_gauss_pts,
574 CommonData &common_data)
575 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
576 common_data.displacementFieldName, UserDataOperator::OPROW),
577 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
578 commonData(common_data) {}
579 enum AxialForces { NX = 0, NY, NXY };
580 enum ShearForces { TZ = 0, TXZ, TYZ };
581 MatrixDouble aXial;
582 MatrixDouble sHear;
583 PetscErrorCode doWork(int side, EntityType type,
584 DataForcesAndSourcesCore::EntData &data);
585 };
586
587 /** \brief Post process moments forces
588 \ingroup solid_shell_prism_element
589 */
592 moab::Interface &postProcMesh;
593 vector<EntityHandle> &mapGaussPts;
595 OpGetMoment(moab::Interface &post_proc_mesh,
596 vector<EntityHandle> &map_gauss_pts, CommonData &common_data)
597 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
598 common_data.displacementFieldName, UserDataOperator::OPROW),
599 postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
600 commonData(common_data) {}
601 enum Moments { MX = 0, MY, MXY };
602 MatrixDouble mOment;
603 PetscErrorCode doWork(int side, EntityType type,
604 DataForcesAndSourcesCore::EntData &data);
605 };
606
609 MakeB {
611 MatrixDouble mF;
613 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
614 common_data.displacementFieldName, UserDataOperator::OPROW),
615 commonData(common_data) {}
616 PetscErrorCode doWork(int side, EntityType type,
617 DataForcesAndSourcesCore::EntData &data);
618 };
619
622 MakeB {
624 MatrixDouble mF;
625 bool zEro;
626 OpGetFLocal(CommonData &common_data, const string &field_name, bool zero)
627 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
629 commonData(common_data), zEro(zero) {}
630 PetscErrorCode doWork(int side, EntityType type,
631 DataForcesAndSourcesCore::EntData &data);
632 };
633
634#ifdef ADOLC_ADOLC_H
635
636 /** \brief Evaluate physical equation and tangent using ADOL-C
637 */
638 struct EvaluateForward : public MakeB {
639
640 CommonData &commonData;
641 double lAmbda, mU;
642 bool convectiveLike; //< current base of shell CS is fixed
643 EvaluateForward(CommonData &common_data, double lambda, double mu)
644 : commonData(common_data), lAmbda(lambda), mU(mu),
645 convectiveLike(true) {}
646
647 template <typename TYPE>
648 PetscErrorCode initializeData(CommonData::ADouble<TYPE> &data) {
649 PetscFunctionBegin;
650 if (data.mGrad1.size1() != 3) {
651 // get current base
652 data.mGrad1.resize(3, 3, false);
653 data.mGrad2.resize(3, 3, false);
654 data.mMF.resize(3, 3, false);
655 data.mF.resize(3, 3, false);
656 data.eM_i_a.resize(3, 3, false);
657 data.eMai.resize(3, 3, false);
658 data.gMab.resize(3, 3, false);
659 data.gM_a_b.resize(3, 3, false);
660 data.mCA_B.resize(3, 3, false);
661 data.mEA_B.resize(3, 3, false);
662 data.mCAB.resize(3, 3, false);
663 data.mEAB.resize(3, 3, false);
664 data.piolaStressIIA_B.resize(3, 3, false);
665 data.localPiolaStressIA_a.resize(3, 3, false);
666 data.localPiolaStressIa_A.resize(3, 3, false);
667 data.globalPiolaStressIi_I.resize(3, 3, false);
668 }
669 PetscFunctionReturn(0);
670 }
671
672 /** \brief Calculate Piola-Kirchhoff II Stress
673
674 \[f
675 S^A_B = \lambda E^A_A + 2\mu E^A_B
676 \f]
677
678 */
679 template <typename TYPE>
680 PetscErrorCode physicalEquation(CommonData::ADouble<TYPE> &data, int ggf,
681 int ggt, int gg, bool global_row,
682 bool global_col) {
683 PetscFunctionBegin;
684 try {
685
686 // Calculate material strain tensor
687 noalias(data.mEA_B) = data.mCA_B;
688 for (int AA = 0; AA != 3; AA++) {
689 int BB = AA;
690 data.mEA_B(AA, BB) -= 1;
691 }
692 data.mEA_B *= 0.5;
693
694 // Calculate stress (physical equation here)
695 data.piolaStressIIA_B = 2 * mU * data.mEA_B;
696 data.tR = 0;
697 for (int AA = 0; AA != 3; AA++) {
698 int BB = AA;
699 data.tR += data.mEA_B(AA, BB);
700 }
701 for (int AA = 0; AA != 3; AA++) {
702 int BB = AA;
703 data.piolaStressIIA_B(AA, BB) += data.tR * lAmbda;
704 }
705
706 /* Piola Local Coords
707 P^c_C = F^c_B S^B_C \\
708 P^A_a = g_{ac} G^{AC} F^c_B S^B_C = g_{ac} G^{AC} P^c_C \\
709 P^A_a = g_{ac} G^{AC} F^c_B S^B_C
710 */
711 data.localPiolaStressIa_A.clear();
712 for (int cc = 0; cc != 3; cc++) {
713 for (int CC = 0; CC != 3; CC++) {
714 TYPE &localPiolaStressIa_A = data.localPiolaStressIa_A(CC, cc);
715 for (int BB = 0; BB != 3; BB++) {
716 localPiolaStressIa_A +=
717 data.mF(cc, BB) * data.piolaStressIIA_B(BB, CC);
718 }
719 }
720 }
721 data.localPiolaStressIA_a.clear();
722 for (int aa = 0; aa != 3; aa++) {
723 for (int AA = 0; AA != 3; AA++) {
724 TYPE &localPiolaStressIA_a = data.localPiolaStressIA_a(aa, AA);
725 for (int cc = 0; cc != 3; cc++) {
726 TYPE &gM_a_b = data.gM_a_b(aa, cc);
727 for (int CC = 0; CC != 3; CC++) {
728 data.t0 = gM_a_b * commonData.mGAB[ggf](AA, CC);
729 data.t1 = data.t0 * data.localPiolaStressIa_A(CC, cc);
730 localPiolaStressIA_a += data.t1;
731 }
732 }
733 }
734 }
735
736 } catch (const std::exception &ex) {
737 ostringstream ss;
738 ss << "throw in method: " << ex.what() << endl;
739 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
740 }
741 PetscFunctionReturn(0);
742 }
743
744 template <typename TYPE>
745 PetscErrorCode evaluateForward(CommonData::ADouble<TYPE> &data, int ggf,
746 int ggt, int gg, bool global_row,
747 bool global_col) {
748 PetscFunctionBegin;
749
750 // FIXME: This is not optimized,
751
752 try {
753
754 if (convectiveLike) {
755
756 noalias(data.mMF) = data.mGrad1;
757 for (int II = 0; II != 3; II++) {
758 int ii = II;
759 data.mMF(ii, II) += 1;
760 }
761
762 data.eMai.clear();
763 for (int aa = 0; aa != 3; aa++) {
764 int AA = aa;
765 for (int ii = 0; ii != 3; ii++) {
766 TYPE &eMai = data.eMai(aa, ii);
767 for (int II = 0; II != 3; II++) {
768 eMai += data.mMF(ii, II) * commonData.mEAI[ggf](AA, II);
769 }
770 }
771 }
772
773 CHKERR dEterminatnt(data.eMai, data.dEt);
774 CHKERR iNvert(data.dEt, data.eMai, data.eM_i_a);
775
776 data.gMab.clear();
777 for (int aa = 0; aa != 3; aa++) {
778 for (int bb = 0; bb != 3; bb++) {
779 TYPE &gMab = data.gM_a_b(aa, bb);
780 for (int ii = 0; ii != 3; ii++) {
781 data.t0 = data.eMai(aa, ii) * data.eMai(bb, ii);
782 gMab += data.t0;
783 }
784 }
785 }
786 data.gM_a_b.clear();
787 for (int aa = 0; aa != 3; aa++) {
788 for (int bb = 0; bb != 3; bb++) {
789 TYPE &gM_a_b = data.gM_a_b(aa, bb);
790 for (int ii = 0; ii != 3; ii++) {
791 data.t0 = data.eM_i_a(ii, aa) * data.eM_i_a(ii, bb);
792 gM_a_b += data.t0;
793 }
794 }
795 }
796
797 } else {
798
799 noalias(data.eMai) = commonData.mEAI[ggf];
800 noalias(data.eM_i_a) = commonData.mE_I_A[ggf];
801 noalias(data.gMab) = commonData.mGAB[ggf];
802 noalias(data.gM_a_b) = commonData.mG_A_B[ggf];
803 }
804
805 // Gradient of defamation (curvilinear system)
806 // Local current coordinate system is convected by
807 // global DoFs
808 data.mF.clear();
809 for (int aa = 0; aa != 3; aa++) {
810 for (int AA = 0; AA != 3; AA++) {
811 TYPE &mF = data.mF(aa, AA);
812 for (int ii = 0; ii != 3; ii++) {
813 int II = ii;
814 mF += data.eMai(aa, ii) * commonData.mE_I_A[ggf](II, AA);
815 }
816 }
817 }
818 if (global_col) {
819 for (int aa = 0; aa != 3; aa++) {
820 for (int AA = 0; AA != 3; AA++) {
821 TYPE &mF = data.mF(aa, AA);
822 for (int ii = 0; ii != 3; ii++) {
823 for (int II = 0; II != 3; II++) {
824 data.t0 = data.eMai(aa, ii) * commonData.mE_I_A[ggf](II, AA);
825 mF += data.t0 * data.mGrad1(ii, II);
826 }
827 }
828 }
829 }
830 noalias(data.mF) += commonData.mGrad2[gg];
831 } else {
832 for (int aa = 0; aa != 3; aa++) {
833 for (int AA = 0; AA != 3; AA++) {
834 TYPE &mF = data.mF(aa, AA);
835 for (int ii = 0; ii != 3; ii++) {
836 for (int II = 0; II != 3; II++) {
837 data.t0 = data.eMai(aa, ii) * commonData.mE_I_A[ggf](II, AA);
838 mF += data.t0 * commonData.mGrad1[gg](ii, II);
839 }
840 }
841 }
842 }
843 noalias(data.mF) += data.mGrad2;
844 }
845
846 // Calculate material deformation tensor
847 // C^A_B = g_{ab}G^{AC} F^b_C F^a_B
848 data.mCA_B.clear();
849 TYPE &mGAB_AA_CCmF_bb_CC = data.t0;
850 TYPE &mGAB_AA_CCmF_bb_CCgM_aa_ba = data.t1;
851 for (int AA = 0; AA != 3; AA++) {
852 for (int CC = 0; CC != 3; CC++) {
853 double mGAB_AA_CC = commonData.mGAB[ggf](AA, CC);
854 for (int bb = 0; bb != 3; bb++) {
855 mGAB_AA_CCmF_bb_CC = mGAB_AA_CC * data.mF(bb, CC);
856 for (int aa = 0; aa != 3; aa++) {
857 mGAB_AA_CCmF_bb_CCgM_aa_ba =
858 mGAB_AA_CCmF_bb_CC * data.gM_a_b(aa, bb);
859 for (int BB = 0; BB != 3; BB++) {
860 data.t3 = mGAB_AA_CCmF_bb_CCgM_aa_ba * data.mF(aa, BB);
861 data.mCA_B(AA, BB) += data.t3;
862 // data.gM_a_b(aa,bb)*
863 // commonData.mGAB[ggf](AA,CC)*
864 // data.mF(bb,CC)*
865 // data.mF(aa,BB);
866 }
867 }
868 }
869 }
870 }
871
872 CHKERR
873 physicalEquation<TYPE>(data, ggf, ggt, gg, global_row, global_col);
874
875 if (global_row) {
876 // Piola Global Coords
877 // P^I_i = E_A^I e^a_i P^A_a
878 data.globalPiolaStressIi_I.clear();
879 for (int ii = 0; ii != 3; ii++) {
880 for (int II = 0; II != 3; II++) {
881 TYPE &globalPiolaStressIi_I = data.globalPiolaStressIi_I(ii, II);
882 for (int aa = 0; aa != 3; aa++) {
883 for (int AA = 0; AA != 3; AA++) {
884 data.t0 = data.eMai(aa, ii) * commonData.mE_I_A[ggf](II, AA);
885 data.t1 = data.t0 * data.localPiolaStressIA_a(aa, AA);
886 globalPiolaStressIi_I += data.t1;
887 }
888 }
889 }
890 }
891 }
892
893 } catch (const std::exception &ex) {
894 ostringstream ss;
895 ss << "throw in method: " << ex.what() << endl;
896 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
897 }
898
899 PetscFunctionReturn(0);
900 }
901 };
902
903 struct OpEvaluatePiolaStress
905 EvaluateForward {
906 int orderThickness;
907 bool getTangent;
908 int tAg;
909 OpEvaluatePiolaStress(CommonData &common_data, double lambda, double mu,
910 int tag, int order_thickness, bool get_tangent)
911 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
912 common_data.displacementFieldName, UserDataOperator::OPROW),
913 EvaluateForward(common_data, lambda, mu),
914 orderThickness(order_thickness), getTangent(get_tangent), tAg(tag) {}
915 PetscErrorCode doWork(int side, EntityType type,
916 DataForcesAndSourcesCore::EntData &data);
917 };
918
919 struct OpAssembleLhsPiolaStress
921 MakeB {
922 CommonData &commonData;
923 MatrixDouble mK, mTransK;
924 MatrixDouble mRowF, mColF;
925 MatrixDouble dP;
926 bool assembleTranspose;
927 OpAssembleLhsPiolaStress(CommonData &common_data,
928 const string &row_field_name,
929 const string &col_field_name,
930 const bool symm = false,
931 const bool assemble_transpose = true)
932 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
933 row_field_name, col_field_name, UserDataOperator::OPROWCOL),
934 commonData(common_data), assembleTranspose(assemble_transpose) {
935 sYmm = symm;
936 }
937 PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
938 EntityType col_type,
939 DataForcesAndSourcesCore::EntData &row_data,
940 DataForcesAndSourcesCore::EntData &col_data);
941 };
942
943 struct OpAssembleRhsPiolaStress
945 MakeB {
946 CommonData &commonData;
947 MatrixDouble mF;
948 VectorDouble vF;
949 OpAssembleRhsPiolaStress(CommonData &common_data, const string &field_name)
950 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
952 commonData(common_data) {}
953 PetscErrorCode doWork(int side, EntityType type,
954 DataForcesAndSourcesCore::EntData &data);
955 };
956
957#endif // ADOLC_ADOLC_H
958
959 /** \brief Evaluate normals and tangents on shell surface
960 \ingroup solid_shell_prism_element
961
962 It is used to create prism geometry for geometry of triangles and
963 making prism by moving nodes of triangles in normal direction (extruding)
964
965 */
967 moab::Interface &mOab;
968 double aNgle;
969 vector<VectorDouble> rEsult;
971 NormalApprox(moab::Interface &moab) : mOab(moab), aNgle(0) {}
972 NormalApprox(moab::Interface &moab, const double angle)
973 : mOab(moab), aNgle(angle) {}
974
975 MatrixDouble rotMatrix;
976 PetscErrorCode getRotMatrix(VectorDouble &u, double phi);
977
978 vector<VectorDouble> &operator()(EntityHandle ent, double x, double y,
979 double z, VectorDouble normal,
980 VectorDouble tangent1,
981 VectorDouble tangent2);
982 };
983
984#ifdef __ARC_LENGTH_TOOLS_HPP__
985
986 struct OpGetIncrementOfDisplacements
988 CommonData &commonData;
989 OpGetIncrementOfDisplacements(CommonData &common_data)
990 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
991 common_data.displacementFieldName, UserDataOperator::OPROW),
992 commonData(common_data) {}
993 PetscErrorCode doWork(int side, EntityType type,
994 DataForcesAndSourcesCore::EntData &data);
995 };
996
997 struct OpNumericalViscousDampingRhs
999 CommonData &commonData;
1000 double &aLpha;
1001 VectorDouble vC;
1002 OpNumericalViscousDampingRhs(CommonData &common_data, double &alpha)
1003 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
1004 common_data.displacementFieldName, UserDataOperator::OPROW),
1005 commonData(common_data), aLpha(alpha) {}
1006 PetscErrorCode doWork(int side, EntityType type,
1007 DataForcesAndSourcesCore::EntData &data);
1008 };
1009
1010 struct OpNumericalViscousDampingLhs
1012 CommonData &commonData;
1013 double &aLpha;
1014 MatrixDouble mC, mCTrans;
1015 OpNumericalViscousDampingLhs(CommonData &common_data, double &alpha)
1016 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
1017 common_data.displacementFieldName, UserDataOperator::OPROWCOL),
1018 commonData(common_data), aLpha(alpha) {
1019 sYmm = true;
1020 }
1021 PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
1022 EntityType col_type,
1023 DataForcesAndSourcesCore::EntData &row_data,
1024 DataForcesAndSourcesCore::EntData &col_data);
1025 };
1026
1027#endif //__ARC_LENGTH_TOOLS_HPP__
1028
1029 /**
1030 * \brief Calculate semi norm
1031 */
1034 MakeB {
1038 : MoFEM::FatPrismElementForcesAndSourcesCore::UserDataOperator(
1039 "DISPLACEMENT", UserDataOperator::OPROW),
1040 commonData(common_data) {
1041 commonData.saveOnTag = true;
1042 }
1043 PetscErrorCode doWork(int side, EntityType type,
1044 DataForcesAndSourcesCore::EntData &data);
1045 };
1046};
1047
1048} // namespace SolidShellModule
1049
1050#endif //__SOLIDSHELLPRISMELEMENT_HPP__
1051
1052/***************************************************************************/ /**
1053 * \defgroup solid_shell_prism_element Solid Shell Prism Element
1054 * \ingroup user_modules
1055 ******************************************************************************/
1056
1057/***************************************************************************/ /**
1058 * \defgroup user_modules User modules
1059 ******************************************************************************/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
#define CHKERR
Inline error check.
constexpr int order
static double phi
@ F
static double lambda
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
constexpr auto field_name
virtual moab::Interface & get_moab()=0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
Postprocess on prism.
vector< MatrixDouble > mGrad2
Mid-surface gradient of deformation.
MatrixDouble materialStiffness
Material stiffness matrix in Voigt notation.
MatrixDouble positionsAtGaussPts
Reference gauss points coordinates.
VectorDouble jacDeterminantAtGaussPts
Determinant of jacobian.
MatrixDouble localDisplacementsAtGaussPts
Local displacement at integration Pts.
VectorDouble surfaceParametrisation
User direction of surface parametrization.
vector< ublas::matrix< double, ublas::column_major > > jacobianInvAtGaussPts
For element local to global.
MatrixDouble globalDisplacentAtGaussPts
Global displacement at integration ptrs.
map< EntityType, map< int, MatrixDouble > > diffGlobN
vector< ublas::symmetric_matrix< double, ublas::lower > > localStress
PetscErrorCode makeFKsiEta(const MatrixAdaptor &diffN, MatrixDouble &F)
PetscErrorCode makeBZeta(const MatrixAdaptor &diffN, MatrixDouble &B)
Calculate strain operator for displacement thought shell thickness.
PetscErrorCode iNvert(TYPE det, ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > &a, ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > &inv_a)
Calculate inverse of 3x3 matrix.
PetscErrorCode makeBDisp(const MatrixAdaptor &diffN, const MatrixDouble &T, MatrixDouble &B)
Calculate strain operator for displacements DoFs.
PetscErrorCode makeBKsiEta(const MatrixAdaptor &diffN, MatrixDouble &B)
Calculate strain operator for non-plannar displacements through shell thickness.
PetscErrorCode makeFDisp(const MatrixAdaptor &diffN, MatrixDouble &F)
Gradient of deformation.
PetscErrorCode dEterminatnt(ublas::matrix< TYPE, ublas::row_major, ublas::bounded_array< TYPE, 9 > > &a, TYPE &det)
PetscErrorCode makeFZeta(const MatrixAdaptor &diffN, MatrixDouble &F)
Evaluate normals and tangents on shell surface.
vector< VectorDouble > & operator()(EntityHandle ent, double x, double y, double z, VectorDouble normal, VectorDouble tangent1, VectorDouble tangent2)
NormalApprox(moab::Interface &moab, const double angle)
PetscErrorCode getRotMatrix(VectorDouble &u, double phi)
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpAssembleLhsMix(CommonData &common_data, const string &row_field_name, const string &col_field_name, const bool symm=false)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Assemble internal forces for local thorough thickness DoFs.
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpAssembleRhsLocal(CommonData &common_data, const string &field)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetAxialForces(moab::Interface &post_proc_mesh, vector< EntityHandle > &map_gauss_pts, CommonData &common_data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetFLocal(CommonData &common_data, const string &field_name, bool zero)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetMoment(moab::Interface &post_proc_mesh, vector< EntityHandle > &map_gauss_pts, CommonData &common_data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Calculate rotation matrices to local (user) element coordinate system.
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Get strain form displacements at integration points.
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Get strain form local through thickness displacements at integration points.
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetStrainLocal(CommonData &common_data, const string &field)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Post process displacements on post-processing mesh.
vector< EntityHandle > & mapGaussPts
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
CommonData & commonData
MatrixDouble meshNodePositions
moab::Interface & postProcMesh
OpPostProcMeshNodeDispalcements(moab::Interface &post_proc_mesh, vector< EntityHandle > &map_gauss_pts, CommonData &common_data)
Post process geometry (mesh nodal positions) on post-processing mesh.
MatrixDouble meshNodePositions
CommonData & commonData
vector< EntityHandle > & mapGaussPts
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpPostProcMeshNodePositions(moab::Interface &post_proc_mesh, vector< EntityHandle > &map_gauss_pts, CommonData &common_data)
moab::Interface & postProcMesh
Post process strain and stresses on post-processing mesh.
moab::Interface & postProcMesh
OpPostProcStressAndStrain(moab::Interface &post_proc_mesh, map< EntityHandle, EntityHandle > &elements_map, vector< EntityHandle > &map_gauss_pts, CommonData &common_data)
map< EntityHandle, EntityHandle > & elementsMap
CommonData & commonData
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
vector< EntityHandle > & mapGaussPts
Generate post-processing triangular mesh for solid-shell element.
PostProcFatPrismOnTriangleOnRefinedMesh(MoFEM::Interface &m_field, bool ten_nodes_post_proc_tets=true)
PetscErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
PetscErrorCode postProcess()
function is run at the end of loop
PetscErrorCode generateReferenceElementMesh()
PetscErrorCode setGaussPtsThroughThickness(int order_thickness)
int getRuleThroughThickness(int order)
int getRuleTrianglesOnly(int order)
PetscErrorCode preProcess()
function is run at the beginning of loop
PetscErrorCode postProcess()
function is run at the end of loop
SolidShellError(MoFEM::Interface &m_field, CommonData &common_data)
PetscErrorCode postProcess()
function is run at the end of loop
PetscErrorCode preProcess()
function is run at the beginning of loop
SolidShell(MoFEM::Interface &m_field, CommonData &common_data)
Main class for solid shell element.