v0.13.1
Loading...
Searching...
No Matches
plate.cpp
Go to the documentation of this file.
1/**
2 * \file plate.cpp
3 * \example plate.cpp
4 *
5 * Implementation Kirchhoff-Love plate using Discointinous Galerkin (DG-Nitsche
6 * method)
7 */
8
9#include <MoFEM.hpp>
10
11using namespace MoFEM;
12
13static char help[] = "...\n\n";
14
16
17constexpr int BASE_DIM = 1; ///< dimension of base
18constexpr int SPACE_DIM = 2; ///< dimension of space
19constexpr int FIELD_DIM = 1; ///< dimension of approx. field
20
21template <int DIM> struct ElementsAndOps {};
22
23template <> struct ElementsAndOps<2> {
28};
29
34
36
39
40using DomainEleOp =
41 DomainEle::UserDataOperator; ///< Finire element operator type
42using EntData = EntitiesFieldData::EntData; ///< Data on entities
43
46
49 GAUSS>::OpGradGradSymTensorGradGrad<1, 1, SPACE_DIM, 0>;
52
53// Kronecker delta
55
56// material parameters
57constexpr double lambda = 1;
58constexpr double mu = 1; ///< lame parameter
59constexpr double t = 1; ///< plate stiffness
60
65
66static double penalty = 1e6;
67static double phi =
68 -1; // 1 - symmetric Nitsche, 0 - nonsymmetric, -1 antisymmetrica
69static double nitsche = 1;
70static int order = 4; // approximation order
71
72auto source = [](const double x, const double y, const double) {
73 return cos(2 * x * M_PI) * sin(2 * y * M_PI);
74};
75
76/**
77 * @brief get fourth-order constitutive tensor
78 *
79 */
80auto plate_stiffness = []() {
81 constexpr auto a = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
82 auto mat_D_ptr = boost::make_shared<MatrixDouble>(a * a, 1);
83 auto t_D = getFTensor4DdgFromMat<2, 2, 0>(*(mat_D_ptr));
84 constexpr double t3 = t * t * t;
85 constexpr double A = mu * t3 / 12;
86 constexpr double B = lambda * t3 / 12;
87 t_D(i, j, k, l) =
88 2 * B * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) + A * t_kd(i, j) * t_kd(k, l);
89 // t_D(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
90 return mat_D_ptr;
91};
92
94
95// data for skeleton computation
96std::array<std::vector<VectorInt>, 2>
97 indicesSideMap; ///< indices on rows for left hand-side
98std::array<std::vector<MatrixDouble>, 2>
99 diffBaseSideMap; // first derivative of base functions
100std::array<std::vector<MatrixDouble>, 2>
101 diff2BaseSideMap; // second derivative of base functions
102std::array<double, 2> areaMap; // area/volume of elements on the side
103std::array<int, 2> senseMap; // orientaton of local element edge/face in respect
104 // to global orientation of edge/face
105
106/**
107 * @brief Operator tp collect data from elements on the side of Edge/Face
108 *
109 */
110struct OpCalculateSideData : public FaceSideOp {
111
112 OpCalculateSideData(std::string field_name, std::string col_field_name);
113
114 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
115};
116
117/**
118 * @brief Operator the left hand side matrix
119 *
120 */
122public:
123 /**
124 * @brief Construct a new OpH1LhsSkeleton
125 *
126 * @param side_fe_ptr pointer to FE to evaluate side elements
127 */
128 OpH1LhsSkeleton(boost::shared_ptr<FaceSideEle> side_fe_ptr,
129 boost::shared_ptr<MatrixDouble> d_mat_ptr);
130
131 MoFEMErrorCode doWork(int side, EntityType type,
133
134private:
135 boost::shared_ptr<FaceSideEle>
136 sideFEPtr; ///< pointer to element to get data on edge/face sides
137 MatrixDouble locMat; ///< local operator matrix
138 boost::shared_ptr<MatrixDouble> dMatPtr;
139};
140
141struct Plate {
142
143 Plate(MoFEM::Interface &m_field) : mField(m_field) {}
144
146
147private:
154
156};
157
158//! [Read mesh]
161
163 CHKERR simple->getOptions();
164 CHKERR simple->loadFile();
165
167}
168//! [Read mesh]
169
170//! [Set up problem]
173
175
176 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
177 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-penalty", &penalty,
178 PETSC_NULL);
179 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-phi", &phi, PETSC_NULL);
180 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-nitsche", &nitsche,
181 PETSC_NULL);
182
183 MOFEM_LOG("WORLD", Sev::inform) << "Set order: " << order;
184 MOFEM_LOG("WORLD", Sev::inform) << "Set penalty: " << penalty;
185 MOFEM_LOG("WORLD", Sev::inform) << "Set phi: " << phi;
186 MOFEM_LOG("WORLD", Sev::inform) << "Set nitche: " << nitsche;
187
188 CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, FIELD_DIM);
189 CHKERR simple->addSkeletonField("U", H1, AINSWORTH_LEGENDRE_BASE, FIELD_DIM);
190 CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, FIELD_DIM);
191
192 CHKERR simple->setFieldOrder("U", order);
193 CHKERR simple->setUp();
194
196}
197//! [Set up problem]
198
199//! [Boundary condition]
202
203 // get edges and vertices on body skin
204 auto get_skin = [&]() {
205 Range body_ents;
207 mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents));
208 Skinner skin(&mField.get_moab());
209 Range skin_ents;
210 MOAB_THROW(skin.find_skin(0, body_ents, false, skin_ents));
211 Range verts;
212 MOAB_THROW(mField.get_moab().get_connectivity(skin_ents, verts, true));
213 skin_ents.merge(verts);
214 return skin_ents;
215 };
216
217 // remove dofs on skin edges from priblem
218 auto remove_dofs_from_problem = [&](Range &&skin) {
220 auto problem_mng = mField.getInterface<ProblemsManager>();
222 CHKERR problem_mng->removeDofsOnEntities(simple->getProblemName(), "U",
223 skin, 0, 1);
225 };
226
227 // it make plate simply supported
228 CHKERR remove_dofs_from_problem(get_skin());
229
231}
232//! [Boundary condition]
233
234//! [Push operators to pipeline]
237
238 auto pipeline_mng = mField.getInterface<PipelineManager>();
239
240 auto rule_lhs = [](int, int, int p) -> int { return 2 * p; };
241 auto rule_rhs = [](int, int, int p) -> int { return 2 * p; };
242 auto rule_2 = [this](int, int, int) { return 2 * order; };
243
244 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
245 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
246
247 CHKERR pipeline_mng->setSkeletonLhsIntegrationRule(rule_2);
248 CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
249 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(rule_2);
250 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule_2);
251
252 auto det_ptr = boost::make_shared<VectorDouble>();
253 auto jac_ptr = boost::make_shared<MatrixDouble>();
254 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
255 auto base_mass_ptr = boost::make_shared<MatrixDouble>();
256 auto data_l2_ptr = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
257
258 auto mat_D_ptr = plate_stiffness();
259
260 auto push_ho_direcatives = [&](auto &pipeline) {
261 pipeline.push_back(new OpBaseDerivativesMass<BASE_DIM>(
262 base_mass_ptr, data_l2_ptr, AINSWORTH_LEGENDRE_BASE, L2));
263 pipeline.push_back(new OpBaseDerivativesNext<BASE_DIM>(
264 BaseDerivatives::SecondDerivative, base_mass_ptr, data_l2_ptr,
266 };
267
268 /**
269 * @brief calculate jacobian
270 *
271 */
272 auto push_jacobian = [&](auto &pipeline) {
273 pipeline.push_back(new OpSetHOWeightsOnFace());
274 pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
275 pipeline.push_back(
276 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
277 // push first base derivatives tp physical element shape
278 pipeline.push_back(new OpSetInvJacH1ForFace<1>(inv_jac_ptr));
279 // push second base directives tp physical element shape
280 pipeline.push_back(new OpSetInvJacH1ForFace<2>(inv_jac_ptr));
281 };
282
283 push_ho_direcatives(pipeline_mng->getOpDomainLhsPipeline());
284 push_jacobian(pipeline_mng->getOpDomainLhsPipeline());
285
286 pipeline_mng->getOpDomainLhsPipeline().push_back(
287 new OpDomainPlateStiffness("U", "U", mat_D_ptr));
288 // pipeline_mng->getOpDomainLhsPipeline().push_back(new OpDomainGradGrad(
289 // "U", "U", [](const double, const double, const double) { return 1; }));
290
291 pipeline_mng->getOpDomainRhsPipeline().push_back(
292 new OpDomainPlateLoad("U", source));
293
294 // Push operators to the Pipeline for Skeleton
295 auto side_fe_ptr = boost::make_shared<FaceSideEle>(mField);
296 push_ho_direcatives(side_fe_ptr->getOpPtrVector());
297 push_jacobian(side_fe_ptr->getOpPtrVector());
298 side_fe_ptr->getOpPtrVector().push_back(new OpCalculateSideData("U", "U"));
299
300 // Push operators to the Pipeline for Skeleton
301 pipeline_mng->getOpSkeletonLhsPipeline().push_back(
302 new OpH1LhsSkeleton(side_fe_ptr, mat_D_ptr));
303
305}
306//! [Push operators to pipeline]
307
308//! [Solve system]
311
312 auto pipeline_mng = mField.getInterface<PipelineManager>();
314
315 auto ksp_solver = pipeline_mng->createKSP();
316 CHKERR KSPSetFromOptions(ksp_solver);
317 CHKERR KSPSetUp(ksp_solver);
318
319 // Create RHS and solution vectors
320 auto dm = simple->getDM();
321 auto F = smartCreateDMVector(dm);
322 auto D = smartVectorDuplicate(F);
323
324 CHKERR KSPSolve(ksp_solver, F, D);
325
326 // Scatter result data on the mesh
327 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
328 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
329 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
330
332}
333//! [Solve system]
334
335//! [Output results]
338
339 auto pipeline_mng = mField.getInterface<PipelineManager>();
340 pipeline_mng->getDomainLhsFE().reset();
341 pipeline_mng->getSkeletonRhsFE().reset();
342 pipeline_mng->getSkeletonLhsFE().reset();
343 pipeline_mng->getBoundaryRhsFE().reset();
344 pipeline_mng->getBoundaryLhsFE().reset();
345
346 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
347
348 auto u_ptr = boost::make_shared<VectorDouble>();
349 post_proc_fe->getOpPtrVector().push_back(
350 new OpCalculateScalarFieldValues("U", u_ptr));
351
353
354 post_proc_fe->getOpPtrVector().push_back(
355
356 new OpPPMap(post_proc_fe->getPostProcMesh(),
357 post_proc_fe->getMapGaussPts(),
358
359 {{"U", u_ptr}},
360
361 {}, {}, {}
362
363 )
364
365 );
366
367 pipeline_mng->getDomainRhsFE() = post_proc_fe;
368 CHKERR pipeline_mng->loopFiniteElements();
369 CHKERR post_proc_fe->writeFile("out_result.h5m");
370
372}
373//! [Output results]
374
375//! [Run program]
378
384 // CHKERR checkResults();
386
388}
389//! [Run program]
390
391int main(int argc, char *argv[]) {
392
393 // Initialisation of MoFEM/PETSc and MOAB data structures
394 const char param_file[] = "param_file.petsc";
396
397 // Add logging channel for example
398 auto core_log = logging::core::get();
399 core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "PL"));
400 LogManager::setLog("PL");
401 MOFEM_LOG_TAG("PL", "plate");
402
403 try {
404
405 //! [Register MoFEM discrete manager in PETSc]
406 DMType dm_name = "DMMOFEM";
407 CHKERR DMRegister_MoFEM(dm_name);
408 //! [Register MoFEM discrete manager in PETSc
409
410 //! [Create MoAB]
411 moab::Core mb_instance; ///< mesh database
412 moab::Interface &moab = mb_instance; ///< mesh database interface
413 //! [Create MoAB]
414
415 //! [Create MoFEM]
416 MoFEM::Core core(moab); ///< finite element database
417 MoFEM::Interface &m_field = core; ///< finite element database insterface
418 //! [Create MoFEM]
419
420 //! [Plate]
421 Plate ex(m_field);
422 CHKERR ex.runProblem();
423 //! [Plate]
424 }
426
428}
429
431 std::string col_field_name)
432 : FaceSideOp(row_field_name, col_field_name, FaceSideOp::OPROW) {}
433
435 EntData &data) {
437
438 const auto nb_in_loop = getFEMethod()->nInTheLoop;
439
440 auto clear = [](auto nb) {
441 indicesSideMap[nb].clear();
442 diffBaseSideMap[nb].clear();
443 diff2BaseSideMap[nb].clear();
444 };
445
446 if (type == MBVERTEX) {
447 areaMap[nb_in_loop] = getMeasure();
448 senseMap[nb_in_loop] = getEdgeSense();
449 if (!nb_in_loop) {
450 clear(0);
451 clear(1);
452 areaMap[1] = 0;
453 senseMap[1] = 0;
454 }
455 }
456
457 const auto nb_dofs = data.getIndices().size();
458 if (nb_dofs) {
459 indicesSideMap[nb_in_loop].push_back(data.getIndices());
460 diffBaseSideMap[nb_in_loop].push_back(
461 data.getN(BaseDerivatives::FirstDerivative));
462 diff2BaseSideMap[nb_in_loop].push_back(
463 data.getN(BaseDerivatives::SecondDerivative));
464 }
465
467}
468
469template <typename T> inline auto get_ntensor(T &base_mat) {
471 &*base_mat.data().begin());
472};
473
474template <typename T> inline auto get_ntensor(T &base_mat, int gg, int bb) {
475 double *ptr = &base_mat(gg, bb);
477};
478
479template <typename T> inline auto get_diff_ntensor(T &base_mat) {
480 double *ptr = &*base_mat.data().begin();
481 return getFTensor1FromPtr<2>(ptr);
482};
483
484template <typename T>
485inline auto get_diff_ntensor(T &base_mat, int gg, int bb) {
486 double *ptr = &base_mat(gg, 2 * bb);
487 return getFTensor1FromPtr<2>(ptr);
488};
489
490template <typename T> inline auto get_diff2_ntensor(T &base_mat) {
491 double *ptr = &*base_mat.data().begin();
493};
494
495template <typename T>
496inline auto get_diff2_ntensor(T &base_mat, int gg, int bb) {
497 double *ptr = &base_mat(gg, 4 * bb);
499};
500
501/**
502 * @brief Construct a new OpH1LhsSkeleton
503 *
504 * @param side_fe_ptr pointer to FE to evaluate side elements
505 */
506OpH1LhsSkeleton::OpH1LhsSkeleton(boost::shared_ptr<FaceSideEle> side_fe_ptr,
507 boost::shared_ptr<MatrixDouble> mat_d_ptr)
508 : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE), sideFEPtr(side_fe_ptr),
509 dMatPtr(mat_d_ptr) {}
510
514
515 // Collect data from side domian elements
516 CHKERR loopSideFaces("dFE", *sideFEPtr);
517 const auto in_the_loop =
518 sideFEPtr->nInTheLoop; // return number of elements on the side
519
520 // calculate penalty
521 const double s = getMeasure() / (areaMap[0] + areaMap[1]);
522 const double p = penalty * s;
523
524 // get normal of the face or edge
525 auto t_normal = getFTensor1Normal();
526 t_normal.normalize();
527
528 // Elastic stiffness tensor (4th rank tensor with minor and major
529 // symmetry)
530 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*dMatPtr);
531
532 // get number of integration points on face
533 const size_t nb_integration_pts = getGaussPts().size2();
534
535 // beta paramter is zero, when penalty method is used, also, takes value 1,
536 // when boundary edge/face is evaluated, and 2 when is skeleton edge/face.
537 const double beta = static_cast<double>(nitsche) / (in_the_loop + 1);
538
539 auto integrate = [&](auto sense_row, auto &row_ind, auto &row_diff,
540 auto &row_diff2, auto sense_col, auto &col_ind,
541 auto &col_diff, auto &col_diff2) {
543
544 // number of dofs, for homogenous approximation this value is
545 // constant.
546 const auto nb_rows = row_ind.size();
547 const auto nb_cols = col_ind.size();
548
549 const auto nb_row_base_functions = row_diff.size2() / SPACE_DIM;
550
551 if (nb_cols) {
552
553 // resize local element matrix
554 locMat.resize(nb_rows, nb_cols, false);
555 locMat.clear();
556
557 // get base functions, and integration weights
558 auto t_diff_row_base = get_diff_ntensor(row_diff);
559 auto t_diff2_row_base = get_diff2_ntensor(row_diff2);
560
561 auto t_w = getFTensor0IntegrationWeight();
562
563 // iterate integration points on face/edge
564 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
565
566 // t_w is integration weight, and measure is element area, or
567 // volume, depending if problem is in 2d/3d.
568 const double alpha = getMeasure() * t_w;
569 auto t_mat = locMat.data().begin();
570
571 // iterate rows
572 size_t rr = 0;
573 for (; rr != nb_rows; ++rr) {
574
576 t_mv(i, j) = t_D(i, j, k, l) * t_diff2_row_base(k, l);
577
578 // calculate tetting function
580 t_vn_plus(i, j) = beta * (phi * t_mv(i, j) / p);
582 t_vn(i, j) = (t_diff_row_base(j) * (t_normal(i) * sense_row)) -
583 t_vn_plus(i, j);
584
585 // get base functions on columns
586 auto t_diff_col_base = get_diff_ntensor(col_diff, gg, 0);
587 auto t_diff2_col_base = get_diff2_ntensor(col_diff2, gg, 0);
588
589 // iterate columns
590 for (size_t cc = 0; cc != nb_cols; ++cc) {
591
593 t_mu(i, j) = t_D(i, j, k, l) * t_diff2_col_base(k, l);
594
595 // // calculate variance of tested function
597 t_un(i, j) = -p * ((t_diff_col_base(j) * (t_normal(i) * sense_col) -
598 beta * t_mu(i, j) / p));
599
600 // assemble matrix
601 *t_mat -= alpha * (t_vn(i, j) * t_un(i, j));
602 *t_mat -= alpha * (t_vn_plus(i, j) * (beta * t_mu(i, j)));
603
604 // move to next column base and element of matrix
605 ++t_diff_col_base;
606 ++t_diff2_col_base;
607 ++t_mat;
608 }
609
610 // move to next row base
611 ++t_diff_row_base;
612 ++t_diff2_row_base;
613 }
614
615 // this is obsolete for this particular example, we keep it for
616 // generality. in case of multi-physcis problems diffrent fields
617 // can chare hierarchical base but use diffrent approx. order,
618 // so is possible to have more base functions than DOFs on
619 // element.
620 for (; rr < nb_row_base_functions; ++rr) {
621 ++t_diff_row_base;
622 ++t_diff2_row_base;
623 }
624
625 ++t_w;
626 }
627
628 // assemble system
629 CHKERR ::MatSetValues(getKSPB(), nb_rows, &*row_ind.begin(),
630 col_ind.size(), &*col_ind.begin(),
631 &*locMat.data().begin(), ADD_VALUES);
632 }
633
635 };
636
637 // iterate the sides rows
638 for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
639
640 const auto sense_row = senseMap[s0];
641
642 for (auto x0 = 0; x0 != indicesSideMap[s0].size(); ++x0) {
643
644 for (auto s1 : {LEFT_SIDE, RIGHT_SIDE}) {
645 const auto sense_col = senseMap[s1];
646
647 for (auto x1 = 0; x1 != indicesSideMap[s1].size(); ++x1) {
648
649 CHKERR integrate(sense_row, indicesSideMap[s0][x0],
650 diffBaseSideMap[s0][x0], diff2BaseSideMap[s0][x0],
651
652 sense_col, indicesSideMap[s1][x1],
653 diffBaseSideMap[s1][x1], diff2BaseSideMap[s1][x1]
654
655 );
656 }
657 }
658
659 }
660 }
661
663}
static Index< 'p', 3 > p
std::string param_file
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:541
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
ElementsAndOps< SPACE_DIM >::PostProcEle PostProcEle
ElementsAndOps< SPACE_DIM >::FaceSideEle FaceSideEle
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition: helmholtz.cpp:27
double D
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 4 >, 2 > getFTensor2SymmetricLowerFromPtr< 2 >(double *ptr)
Definition: Templates.hpp:947
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr AssemblyType A
Definition: plastic.cpp:46
constexpr double mu
lame parameter
Definition: plate.cpp:58
auto get_diff_ntensor(T &base_mat)
Definition: plate.cpp:479
@ LEFT_SIDE
Definition: plate.cpp:93
@ RIGHT_SIDE
Definition: plate.cpp:93
std::array< std::vector< MatrixDouble >, 2 > diff2BaseSideMap
Definition: plate.cpp:101
FTensor::Index< 'j', SPACE_DIM > j
Definition: plate.cpp:62
static char help[]
Definition: plate.cpp:13
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGradSymTensorGradGrad< 1, 1, SPACE_DIM, 0 > OpDomainPlateStiffness
Definition: plate.cpp:49
std::array< std::vector< MatrixDouble >, 2 > diffBaseSideMap
Definition: plate.cpp:99
FTensor::Index< 'k', SPACE_DIM > k
Definition: plate.cpp:63
constexpr double lambda
Definition: plate.cpp:57
FTensor::Index< 'i', SPACE_DIM > i
Definition: plate.cpp:61
auto get_ntensor(T &base_mat)
Definition: plate.cpp:469
constexpr auto t_kd
Definition: plate.cpp:54
constexpr int SPACE_DIM
dimension of space
Definition: plate.cpp:18
FTensor::Index< 'l', SPACE_DIM > l
Definition: plate.cpp:64
auto source
Definition: plate.cpp:72
static double nitsche
Definition: plate.cpp:69
constexpr int FIELD_DIM
dimension of approx. field
Definition: plate.cpp:19
auto plate_stiffness
get fourth-order constitutive tensor
Definition: plate.cpp:80
std::array< double, 2 > areaMap
Definition: plate.cpp:102
std::array< std::vector< VectorInt >, 2 > indicesSideMap
indices on rows for left hand-side
Definition: plate.cpp:97
std::array< int, 2 > senseMap
Definition: plate.cpp:103
constexpr int BASE_DIM
dimension of base
Definition: plate.cpp:17
constexpr double t
plate stiffness
Definition: plate.cpp:59
static double penalty
Definition: plate.cpp:66
static int order
Definition: plate.cpp:70
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpDomainPlateLoad
Definition: plate.cpp:51
static double phi
Definition: plate.cpp:67
auto get_diff2_ntensor(T &base_mat)
Definition: plate.cpp:490
constexpr auto field_name
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Base face element used to integrate on skeleton.
auto getFTensor0IntegrationWeight()
Get integration weights.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Operator tp collect data from elements on the side of Edge/Face.
OpCalculateSideData(std::string field_name, std::string col_field_name, UserDataOperator::OpType type)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator the left hand side matrix.
Definition: plate.cpp:121
boost::shared_ptr< MatrixDouble > dMatPtr
Definition: plate.cpp:138
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
Definition: plate.cpp:136
MatrixDouble locMat
local operator matrix
Definition: plate.cpp:137
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: plate.cpp:511
OpH1LhsSkeleton(boost::shared_ptr< FaceSideEle > side_fe_ptr, boost::shared_ptr< MatrixDouble > d_mat_ptr)
Construct a new OpH1LhsSkeleton.
Definition: plate.cpp:506
Definition: plate.cpp:141
MoFEMErrorCode setupProblem()
[Read mesh]
Definition: plate.cpp:171
MoFEM::Interface & mField
Definition: plate.cpp:155
MoFEMErrorCode runProblem()
[Output results]
Definition: plate.cpp:376
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: plate.cpp:200
MoFEMErrorCode outputResults()
[Solve system]
Definition: plate.cpp:336
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
Definition: plate.cpp:309
MoFEMErrorCode assembleSystem()
[Boundary condition]
Definition: plate.cpp:235
Plate(MoFEM::Interface &m_field)
Definition: plate.cpp:143
MoFEMErrorCode readMesh()
[Read mesh]
Definition: plate.cpp:159