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