v0.14.0
Loading...
Searching...
No Matches
simple_elasticity.cpp
Go to the documentation of this file.
1/** \file simple_elasticity.cpp
2 * \example simple_elasticity.cpp
3
4 The example shows how to solve the linear elastic problem.
5
6*/
7
9
10using namespace boost::numeric;
11using namespace MoFEM;
12
13static char help[] = "-order approximation order\n"
14 "\n";
15
16struct OpK : public VolumeElementForcesAndSourcesCore::UserDataOperator {
17
18 // Finite element stiffness sub-matrix K_ij
20
21 // Elastic stiffness tensor (4th rank tensor with minor and major symmetry)
23
24 // Young's modulus
25 double yOung;
26 // Poisson's ratio
27 double pOisson;
28
29 OpK(bool symm = true)
31 symm) {
32
33 // Evaluation of the elastic stiffness tensor, D
34
35 // hardcoded choice of elastic parameters
36 pOisson = 0.1;
37 yOung = 10;
38
39 // coefficient used in intermediate calculation
40 const double coefficient = yOung / ((1 + pOisson) * (1 - 2 * pOisson));
41
42 FTensor::Index<'i', 3> i;
43 FTensor::Index<'j', 3> j;
44 FTensor::Index<'k', 3> k;
45 FTensor::Index<'l', 3> l;
46
47 tD(i, j, k, l) = 0.;
48
49 tD(0, 0, 0, 0) = 1 - pOisson;
50 tD(1, 1, 1, 1) = 1 - pOisson;
51 tD(2, 2, 2, 2) = 1 - pOisson;
52
53 tD(0, 1, 0, 1) = 0.5 * (1 - 2 * pOisson);
54 tD(0, 2, 0, 2) = 0.5 * (1 - 2 * pOisson);
55 tD(1, 2, 1, 2) = 0.5 * (1 - 2 * pOisson);
56
57 tD(0, 0, 1, 1) = pOisson;
58 tD(1, 1, 0, 0) = pOisson;
59 tD(0, 0, 2, 2) = pOisson;
60 tD(2, 2, 0, 0) = pOisson;
61 tD(1, 1, 2, 2) = pOisson;
62 tD(2, 2, 1, 1) = pOisson;
63
64 tD(i, j, k, l) *= coefficient;
65 }
66
67 /**
68 * \brief Do calculations for give operator
69 * @param row_side row side number (local number) of entity on element
70 * @param col_side column side number (local number) of entity on element
71 * @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
72 * @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
73 * @param row_data data for row
74 * @param col_data data for column
75 * @return error code
76 */
77 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
78 EntityType col_type,
81
83
84 // get number of dofs on row
85 nbRows = row_data.getIndices().size();
86 // if no dofs on row, exit that work, nothing to do here
87 if (!nbRows)
89
90 // get number of dofs on column
91 nbCols = col_data.getIndices().size();
92 // if no dofs on Columbia, exit nothing to do here
93 if (!nbCols)
95
96 // K_ij matrix will have 3 times the number of degrees of freedom of the
97 // i-th entity set (nbRows)
98 // and 3 times the number of degrees of freedom of the j-th entity set
99 // (nbCols)
100 K.resize(nbRows, nbCols, false);
101 K.clear();
102
103 // get number of integration points
104 nbIntegrationPts = getGaussPts().size2();
105 // check if entity block is on matrix diagonal
106 if (row_side == col_side && row_type == col_type) {
107 isDiag = true;
108 } else {
109 isDiag = false;
110 }
111
112 // integrate local matrix for entity block
113 CHKERR iNtegrate(row_data, col_data);
114
115 // assemble local matrix
116 CHKERR aSsemble(row_data, col_data);
117
119 }
120
121protected:
122 int nbRows; ///< number of dofs on rows
123 int nbCols; ///< number if dof on column
124 int nbIntegrationPts; ///< number of integration points
125 bool isDiag; ///< true if this block is on diagonal
126
127 /**
128 * \brief Integrate B^T D B operator
129 * @param row_data row data (consist base functions on row entity)
130 * @param col_data column data (consist base functions on column entity)
131 * @return error code
132 */
134 EntitiesFieldData::EntData &col_data) {
136
137 // get sub-block (3x3) of local stiffens matrix, here represented by second
138 // order tensor
139 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
141 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
142 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
143 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
144 };
145
146 FTensor::Index<'i', 3> i;
147 FTensor::Index<'j', 3> j;
148 FTensor::Index<'k', 3> k;
149 FTensor::Index<'l', 3> l;
150
151 // get element volume
152 double vol = getVolume();
153
154 // get intergration weights
155 auto t_w = getFTensor0IntegrationWeight();
156
157 // get derivatives of base functions on rows
158 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
159 // iterate over integration points
160 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
161
162 // calculate scalar weight times element volume
163 const double a = t_w * vol;
164
165 // iterate over row base functions
166 for (int rr = 0; rr != nbRows / 3; ++rr) {
167
168 // get sub matrix for the row
169 auto t_m = get_tensor2(K, 3 * rr, 0);
170
171 // get derivatives of base functions for columns
172 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
173
174 // iterate column base functions
175 for (int cc = 0; cc != nbCols / 3; ++cc) {
176
177 // integrate block local stiffens matrix
178 t_m(i, k) +=
179 a * (tD(i, j, k, l) * (t_row_diff_base(j) * t_col_diff_base(l)));
180
181 // move to next column base function
182 ++t_col_diff_base;
183
184 // move to next block of local stiffens matrix
185 ++t_m;
186 }
187
188 // move to next row base function
189 ++t_row_diff_base;
190 }
191
192 // move to next integration weight
193 ++t_w;
194 }
195
197 }
198
199 /**
200 * \brief Assemble local entity block matrix
201 * @param row_data row data (consist base functions on row entity)
202 * @param col_data column data (consist base functions on column entity)
203 * @return error code
204 */
206 EntitiesFieldData::EntData &col_data) {
208 // get pointer to first global index on row
209 const int *row_indices = &*row_data.getIndices().data().begin();
210 // get pointer to first global index on column
211 const int *col_indices = &*col_data.getIndices().data().begin();
212 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
213 : getFEMethod()->snes_B;
214 // assemble local matrix
215 CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
216 &*K.data().begin(), ADD_VALUES);
217
218 if (!isDiag && sYmm) {
219 // if not diagonal term and since global matrix is symmetric assemble
220 // transpose term.
221 K = trans(K);
222 CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
223 &*K.data().begin(), ADD_VALUES);
224 }
226 }
227};
228
230
232
233 OpPressure(const double pressure_val = 1)
235 pressureVal(pressure_val) {}
236
237 // vector used to store force vector for each degree of freedom
239
241
244
246 // check that the faces have associated degrees of freedom
247 const int nb_dofs = data.getIndices().size();
248 if (nb_dofs == 0)
250
251 // size of force vector associated to the entity
252 // set equal to the number of degrees of freedom of associated with the
253 // entity
254 nF.resize(nb_dofs, false);
255 nF.clear();
256
257 // get number of gauss points
258 const int nb_gauss_pts = data.getN().size1();
259
260 // create a 3d vector to be used as the normal to the face with length equal
261 // to the face area
262 auto t_normal = getFTensor1Normal();
263
264 // get intergration weights
265 auto t_w = getFTensor0IntegrationWeight();
266
267 // vector of base functions
268 auto t_base = data.getFTensor0N();
269
270 // loop over all gauss points of the face
271 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
272 // weight of gg gauss point
273 double w = 0.5 * t_w;
274
275 // create a vector t_nf whose pointer points an array of 3 pointers
276 // pointing to nF memory location of components
278 &nF[2]);
279 for (int bb = 0; bb != nb_dofs / 3; ++bb) {
280 // scale the three components of t_normal and pass them to the t_nf
281 // (hence to nF)
282 t_nf(i) += (w * pressureVal * t_base) * t_normal(i);
283 // move the pointer to next element of t_nf
284 ++t_nf;
285 // move to next base function
286 ++t_base;
287 }
288
289 // move to next integration weight
290 ++t_w;
291 }
292
293 // add computed values of pressure in the global right hand side vector
294 CHKERR VecSetValues(getFEMethod()->ksp_f, nb_dofs, &data.getIndices()[0],
295 &nF[0], ADD_VALUES);
296
298 }
299};
300
302
304
305 ApplyDirichletBc(const Range &fix_faces, const Range &fix_nodes,
306 const Range &fix_second_node)
307 : MoFEM::FEMethod(), fixFaces(fix_faces), fixNodes(fix_nodes),
308 fixSecondNode(fix_second_node) {
309 // constructor
310 }
311
313
315 std::set<int> set_fix_dofs;
316
318 if (dit->get()->getDofCoeffIdx() == 2) {
319 if (fixFaces.find(dit->get()->getEnt()) != fixFaces.end()) {
320 set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
321 }
322 }
323
324 if (fixSecondNode.find(dit->get()->getEnt()) != fixSecondNode.end()) {
325 if (dit->get()->getDofCoeffIdx() == 1) {
326 set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
327 }
328 }
329
330 if (fixNodes.find(dit->get()->getEnt()) != fixNodes.end()) {
331 set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
332 }
333 }
334
335 std::vector<int> fix_dofs(set_fix_dofs.size());
336
337 std::copy(set_fix_dofs.begin(), set_fix_dofs.end(), fix_dofs.begin());
338
339 CHKERR MatAssemblyBegin(ksp_B, MAT_FINAL_ASSEMBLY);
340 CHKERR MatAssemblyEnd(ksp_B, MAT_FINAL_ASSEMBLY);
341 CHKERR VecAssemblyBegin(ksp_f);
342 CHKERR VecAssemblyEnd(ksp_f);
343
344 Vec x;
345
346 CHKERR VecDuplicate(ksp_f, &x);
347 CHKERR VecZeroEntries(x);
348 CHKERR MatZeroRowsColumns(ksp_B, fix_dofs.size(), &fix_dofs[0], 1, x,
349 ksp_f);
350
351 CHKERR VecDestroy(&x);
352
354 }
355};
356
357/**
358 * \brief Set integration rule to volume elements
359 *
360 * This rule is used to integrate \f$\nabla v \cdot \nabla u\f$, thus
361 * if the approximation field and the testing field are polynomials of order
362 * "p", then the rule for the exact integration is 2*(p-1).
363 *
364 * Integration rule is order of polynomial which is calculated exactly. Finite
365 * element selects integration method based on return of this function.
366 *
367 */
368struct VolRule {
369 int operator()(int, int, int p) const { return 2 * (p - 1); }
370};
371
372/**
373 * \brief Set integration rule to boundary elements
374 *
375 * This rule is used to integrate the work of external forces on a face,
376 * i.e. \f$f \cdot v\f$, where f is the traction vector and v is the test
377 * vector function. The current problem features a Neumann bc with
378 * a pre-defined constant pressure. Therefore, if the test field is
379 * represented by polynomials of order "p", then the rule for the exact
380 * integration is also p.
381 *
382 * Integration rule is order of polynomial which is calculated exactly. Finite
383 * element selects integration method based on return of this function.
384 *
385 */
386struct FaceRule {
387 int operator()(int, int, int p) const { return p; }
388};
389
390int main(int argc, char *argv[]) {
391
392 const string default_options = "-ksp_type fgmres \n"
393 "-pc_type lu \n"
394 "-pc_factor_mat_solver_type mumps \n"
395 "-mat_mumps_icntl_20 0 \n"
396 "-ksp_monitor\n";
397
398 string param_file = "param_file.petsc";
399 if (!static_cast<bool>(ifstream(param_file))) {
400 std::ofstream file(param_file.c_str(), std::ios::ate);
401 if (file.is_open()) {
402 file << default_options;
403 file.close();
404 }
405 }
406
407 MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
408
409 // Create mesh database
410 moab::Core mb_instance; // create database
411 moab::Interface &moab = mb_instance; // create interface to database
412
413 try {
414 // Create MoFEM database and link it to MoAB
415 MoFEM::Core core(moab);
416 MoFEM::Interface &m_field = core;
417
418 CHKERR DMRegister_MoFEM("DMMOFEM");
419
420 // Get command line options
421 int order = 3; // default approximation order
422 PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
423 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "SimpleElasticProblem",
424 "none");
425
426 // Set approximation order
427 CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
428 PETSC_NULL);
429
430 // Set testing (used by CTest)
431 CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
432 &flg_test, PETSC_NULL);
433
434 ierr = PetscOptionsEnd();
435 CHKERRG(ierr);
436
437 Simple *simple_interface = m_field.getInterface<MoFEM::Simple>();
438
439 CHKERR simple_interface->getOptions();
440 CHKERR simple_interface->loadFile();
441
442 Range fix_faces, pressure_faces, fix_nodes, fix_second_node;
443
444 enum MyBcTypes {
445 FIX_BRICK_FACES = 1,
446 FIX_NODES = 2,
447 BRICK_PRESSURE_FACES = 3,
448 FIX_NODES_Y_DIR = 4
449 };
450
452 EntityHandle meshset = bit->getMeshset();
453 int id = bit->getMeshsetId();
454
455 if (id == FIX_BRICK_FACES) { // brick-faces
456
457 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 2,
458 fix_faces, true);
459
460 Range adj_ents;
461 CHKERR m_field.get_moab().get_adjacencies(fix_faces, 0, false, adj_ents,
462 moab::Interface::UNION);
463
464 CHKERR m_field.get_moab().get_adjacencies(fix_faces, 1, false, adj_ents,
465 moab::Interface::UNION);
466 fix_faces.merge(adj_ents);
467 } else if (id == FIX_NODES) { // node(s)
468
469 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 0,
470 fix_nodes, true);
471
472 } else if (id == BRICK_PRESSURE_FACES) { // brick pressure faces
473 CHKERR m_field.get_moab().get_entities_by_dimension(
474 meshset, 2, pressure_faces, true);
475
476 } else if (id ==
477 FIX_NODES_Y_DIR) { // restrained second node in y direction
478 CHKERR m_field.get_moab().get_entities_by_dimension(
479 meshset, 0, fix_second_node, true);
480
481 } else {
482 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Unknown blockset");
483 }
484 }
485
486 CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
487 3);
488 CHKERR simple_interface->setFieldOrder("U", order);
489
490 CHKERR simple_interface->defineFiniteElements();
491
492 // Add pressure element
493 CHKERR m_field.add_finite_element("PRESSURE");
494 CHKERR m_field.modify_finite_element_add_field_row("PRESSURE", "U");
495 CHKERR m_field.modify_finite_element_add_field_col("PRESSURE", "U");
496 CHKERR m_field.modify_finite_element_add_field_data("PRESSURE", "U");
497
498 CHKERR simple_interface->defineProblem();
499
500 DM dm;
501 CHKERR simple_interface->getDM(&dm);
502
503 CHKERR DMMoFEMAddElement(dm, "PRESSURE");
504 CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
505
506 CHKERR simple_interface->buildFields();
507 CHKERR simple_interface->buildFiniteElements();
508
509 CHKERR m_field.add_ents_to_finite_element_by_dim(pressure_faces, 2,
510 "PRESSURE");
511 CHKERR m_field.build_finite_elements("PRESSURE", &pressure_faces);
512
513 CHKERR simple_interface->buildProblem();
514
515 // Create elements instances
516 boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(
518 boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
520
521 // Set integration rule to elements instances
522 elastic_fe->getRuleHook = VolRule();
523 pressure_fe->getRuleHook = FaceRule();
524
525 // Add operators to element instances
526 // push operators to elastic_fe
527 elastic_fe->getOpPtrVector().push_back(new OpK());
528 // push operators to pressure_fe
529 pressure_fe->getOpPtrVector().push_back(new OpPressure());
530
531 boost::shared_ptr<FEMethod> fix_dofs_fe(
532 new ApplyDirichletBc(fix_faces, fix_nodes, fix_second_node));
533
534 boost::shared_ptr<FEMethod> null_fe;
535
536 // Set operators for KSP solver
538 dm, simple_interface->getDomainFEName(), elastic_fe, null_fe, null_fe);
539
540 CHKERR DMMoFEMKSPSetComputeRHS(dm, "PRESSURE", pressure_fe, null_fe,
541 null_fe);
542
543 // initialise matrix A used as the global stiffness matrix
544 Mat A;
545
546 // initialise left hand side vector x and right hand side vector f
547 Vec x, f;
548
549 // allocate memory handled by MoFEM discrete manager for matrix A
550 CHKERR DMCreateMatrix(dm, &A);
551
552 // allocate memory handled by MoFEM discrete manager for vector x
553 CHKERR DMCreateGlobalVector(dm, &x);
554
555 // allocate memory handled by MoFEM discrete manager for vector f of the
556 // same size as x
557 CHKERR VecDuplicate(x, &f);
558
559 // precondition matrix A according to fix_dofs_fe and elastic_fe finite
560 // elements
561 elastic_fe->ksp_B = A;
562 fix_dofs_fe->ksp_B = A;
563
564 // precondition the right hand side vector f according to fix_dofs_fe and
565 // elastic_fe finite elements
566 fix_dofs_fe->ksp_f = f;
567 pressure_fe->ksp_f = f;
568
569 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
570 elastic_fe);
571
572 CHKERR DMoFEMLoopFiniteElements(dm, "PRESSURE", pressure_fe);
573
574 // This is done because only post processor is implemented in the
575 // ApplyDirichletBc struct
576 CHKERR DMoFEMPostProcessFiniteElements(dm, fix_dofs_fe.get());
577
578 // make available a KSP solver
579 KSP solver;
580
581 // make the solver available for parallel computing by determining its MPI
582 // communicator
583 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
584
585 // making available running all options available for KSP solver in running
586 // command
587 CHKERR KSPSetFromOptions(solver);
588
589 // set A matrix with preconditioner
590 CHKERR KSPSetOperators(solver, A, A);
591
592 // set up the solver data strucure for the iterative solver
593 CHKERR KSPSetUp(solver);
594
595 // solve the system of linear equations
596 CHKERR KSPSolve(solver, f, x);
597
598 // destroy solver no needed any more
599 CHKERR KSPDestroy(&solver);
600
601 // make vector x available for parallel computations for visualization
602 // context
603 VecView(x, PETSC_VIEWER_STDOUT_WORLD);
604
605 // save solution in vector x on mesh
606 CHKERR DMoFEMMeshToGlobalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE);
607
608 // Set up post-processor. It is some generic implementation of finite
609 // element
610 auto get_post_proc_ele = [&]() {
611 auto jac_ptr = boost::make_shared<MatrixDouble>();
612 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
613 auto det_ptr = boost::make_shared<VectorDouble>();
614
615 auto post_proc_ele = boost::make_shared<
617 // Add operators to the elements, starting with some generic
618 constexpr int SPACE_DIM = 3;
619 post_proc_ele->getOpPtrVector().push_back(
620 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
621 post_proc_ele->getOpPtrVector().push_back(
622 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
623 post_proc_ele->getOpPtrVector().push_back(
624 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
625
626 auto u_ptr = boost::make_shared<MatrixDouble>();
627 auto grad_ptr = boost::make_shared<MatrixDouble>();
628 post_proc_ele->getOpPtrVector().push_back(
630 post_proc_ele->getOpPtrVector().push_back(
632 grad_ptr));
634
635 post_proc_ele->getOpPtrVector().push_back(
636
637 new OpPPMap(
638
639 post_proc_ele->getPostProcMesh(), post_proc_ele->getMapGaussPts(),
640 {},
641
642 {{"U", u_ptr}},
643
644 {{"GRAD", grad_ptr}},
645
646 {})
647
648 );
649
650 return post_proc_ele;
651 };
652
653 if (auto post_proc = get_post_proc_ele()) {
654 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
655 post_proc);
656 // write output
657 CHKERR post_proc->writeFile("out.h5m");
658 }
659
660 {
661 if (flg_test == PETSC_TRUE) {
662 const double x_vec_norm_const = 0.4;
663 // Check norm_1 value
664 double norm_check;
665 // Takes maximal element of the vector, that should be maximal
666 // displacement at the end of the bar
667 CHKERR VecNorm(x, NORM_INFINITY, &norm_check);
668 if (std::abs(norm_check - x_vec_norm_const) / x_vec_norm_const >
669 1.e-10) {
670 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
671 "test failed (nrm 2 %6.4e)", norm_check);
672 }
673 }
674 }
675
676 // free memory handled by mofem discrete manager for A, x and f
677 CHKERR MatDestroy(&A);
678 CHKERR VecDestroy(&x);
679 CHKERR VecDestroy(&f);
680
681 // free memory allocated for mofem discrete manager
682 CHKERR DMDestroy(&dm);
683
684 // This is a good reference for the future
685 }
687
688 // finish work cleaning memory, getting statistics, etc
690
691 return 0;
692}
static Index< 'p', 3 > p
const std::string default_options
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
int main()
constexpr double a
constexpr int SPACE_DIM
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ H1
continuous field
Definition definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ BLOCKSET
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int order
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1109
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:483
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:542
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:47
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set KSP right hand side evaluation function
Definition DMMoFEM.cpp:623
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:572
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition DMMoFEM.cpp:521
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
Definition DMMoFEM.cpp:664
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_NUMEREDDOF_ROW_FOR_LOOP_(PROBLEMPTR, IT)
use with loops to iterate row DOFs
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr AssemblyType A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
static char help[]
ApplyDirichletBc(const Range &fix_faces, const Range &fix_nodes, const Range &fix_second_node)
MoFEMErrorCode postProcess()
function is run at the end of loop
Set integration rule to boundary elements.
int operator()(int, int, int p) const
const Problem * problemPtr
raw pointer to problem
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)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure for User Loop Methods on finite elements
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEMErrorCode buildProblem()
Build problem.
Definition Simple.cpp:597
MoFEMErrorCode defineFiniteElements()
Define finite elements.
Definition Simple.cpp:380
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:194
MoFEMErrorCode buildFiniteElements()
Build finite elements.
Definition Simple.cpp:571
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition Simple.cpp:264
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:667
MoFEMErrorCode buildFields()
Build fields.
Definition Simple.cpp:482
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:473
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
Definition Simple.cpp:452
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:341
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
double yOung
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate B^T D B operator.
double pOisson
bool isDiag
true if this block is on diagonal
OpK(bool symm=true)
FTensor::Ddg< double, 3, 3 > tD
int nbIntegrationPts
number of integration points
int nbRows
number of dofs on rows
MatrixDouble K
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations for give operator.
int nbCols
number if dof on column
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble local entity block matrix.
VectorDouble nF
FTensor::Index< 'i', 3 > i
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpPressure(const double pressure_val=1)
Set integration rule.
int operator()(int, int, int p) const