v0.15.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
8#include <MoFEM.hpp>
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_NULLPTR ? 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
242 MoFEMErrorCode doWork(int side, EntityType type,
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 PetscOptionsBegin(PETSC_COMM_WORLD, "", "SimpleElasticProblem",
424 "none");
425
426 // Set approximation order
427 CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
428 PETSC_NULLPTR);
429
430 // Set testing (used by CTest)
431 CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
432 &flg_test, PETSC_NULLPTR);
433
434 PetscOptionsEnd();
435
436 Simple *simple_interface = m_field.getInterface<MoFEM::Simple>();
437
438 CHKERR simple_interface->getOptions();
439 CHKERR simple_interface->loadFile();
440
441 Range fix_faces, pressure_faces, fix_nodes, fix_second_node;
442
443 enum MyBcTypes {
444 FIX_BRICK_FACES = 1,
445 FIX_NODES = 2,
446 BRICK_PRESSURE_FACES = 3,
447 FIX_NODES_Y_DIR = 4
448 };
449
451 EntityHandle meshset = bit->getMeshset();
452 int id = bit->getMeshsetId();
453
454 if (id == FIX_BRICK_FACES) { // brick-faces
455
456 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 2,
457 fix_faces, true);
458
459 Range adj_ents;
460 CHKERR m_field.get_moab().get_adjacencies(fix_faces, 0, false, adj_ents,
461 moab::Interface::UNION);
462
463 CHKERR m_field.get_moab().get_adjacencies(fix_faces, 1, false, adj_ents,
464 moab::Interface::UNION);
465 fix_faces.merge(adj_ents);
466 } else if (id == FIX_NODES) { // node(s)
467
468 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 0,
469 fix_nodes, true);
470
471 } else if (id == BRICK_PRESSURE_FACES) { // brick pressure faces
472 CHKERR m_field.get_moab().get_entities_by_dimension(
473 meshset, 2, pressure_faces, true);
474
475 } else if (id ==
476 FIX_NODES_Y_DIR) { // restrained second node in y direction
477 CHKERR m_field.get_moab().get_entities_by_dimension(
478 meshset, 0, fix_second_node, true);
479
480 } else {
481 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Unknown blockset");
482 }
483 }
484
485 CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
486 3);
487 CHKERR simple_interface->setFieldOrder("U", order);
488
489 CHKERR simple_interface->defineFiniteElements();
490
491 // Add pressure element
492 CHKERR m_field.add_finite_element("PRESSURE");
493 CHKERR m_field.modify_finite_element_add_field_row("PRESSURE", "U");
494 CHKERR m_field.modify_finite_element_add_field_col("PRESSURE", "U");
495 CHKERR m_field.modify_finite_element_add_field_data("PRESSURE", "U");
496
497 CHKERR simple_interface->defineProblem();
498
499 DM dm;
500 CHKERR simple_interface->getDM(&dm);
501
502 CHKERR DMMoFEMAddElement(dm, "PRESSURE");
503 CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
504
505 CHKERR simple_interface->buildFields();
506 CHKERR simple_interface->buildFiniteElements();
507
508 CHKERR m_field.add_ents_to_finite_element_by_dim(pressure_faces, 2,
509 "PRESSURE");
510 CHKERR m_field.build_finite_elements("PRESSURE", &pressure_faces);
511
512 CHKERR simple_interface->buildProblem();
513
514 // Create elements instances
515 boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(
517 boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
519
520 // Set integration rule to elements instances
521 elastic_fe->getRuleHook = VolRule();
522 pressure_fe->getRuleHook = FaceRule();
523
524 // Add operators to element instances
525 // push operators to elastic_fe
526 elastic_fe->getOpPtrVector().push_back(new OpK());
527 // push operators to pressure_fe
528 pressure_fe->getOpPtrVector().push_back(new OpPressure());
529
530 boost::shared_ptr<FEMethod> fix_dofs_fe(
531 new ApplyDirichletBc(fix_faces, fix_nodes, fix_second_node));
532
533 boost::shared_ptr<FEMethod> null_fe;
534
535 // Set operators for KSP solver
537 dm, simple_interface->getDomainFEName(), elastic_fe, null_fe, null_fe);
538
539 CHKERR DMMoFEMKSPSetComputeRHS(dm, "PRESSURE", pressure_fe, null_fe,
540 null_fe);
541
542 // initialise matrix A used as the global stiffness matrix
543 Mat A;
544
545 // initialise left hand side vector x and right hand side vector f
546 Vec x, f;
547
548 // allocate memory handled by MoFEM discrete manager for matrix A
549 CHKERR DMCreateMatrix(dm, &A);
550
551 // allocate memory handled by MoFEM discrete manager for vector x
552 CHKERR DMCreateGlobalVector(dm, &x);
553
554 // allocate memory handled by MoFEM discrete manager for vector f of the
555 // same size as x
556 CHKERR VecDuplicate(x, &f);
557
558 // precondition matrix A according to fix_dofs_fe and elastic_fe finite
559 // elements
560 elastic_fe->ksp_B = A;
561 fix_dofs_fe->ksp_B = A;
562
563 // precondition the right hand side vector f according to fix_dofs_fe and
564 // elastic_fe finite elements
565 fix_dofs_fe->ksp_f = f;
566 pressure_fe->ksp_f = f;
567
568 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
569 elastic_fe);
570
571 CHKERR DMoFEMLoopFiniteElements(dm, "PRESSURE", pressure_fe);
572
573 // This is done because only post processor is implemented in the
574 // ApplyDirichletBc struct
575 CHKERR DMoFEMPostProcessFiniteElements(dm, fix_dofs_fe.get());
576
577 // make available a KSP solver
578 KSP solver;
579
580 // make the solver available for parallel computing by determining its MPI
581 // communicator
582 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
583
584 // making available running all options available for KSP solver in running
585 // command
586 CHKERR KSPSetFromOptions(solver);
587
588 // set A matrix with preconditioner
589 CHKERR KSPSetOperators(solver, A, A);
590
591 // set up the solver data strucure for the iterative solver
592 CHKERR KSPSetUp(solver);
593
594 // solve the system of linear equations
595 CHKERR KSPSolve(solver, f, x);
596
597 // destroy solver no needed any more
598 CHKERR KSPDestroy(&solver);
599
600 // make vector x available for parallel computations for visualization
601 // context
602 VecView(x, PETSC_VIEWER_STDOUT_WORLD);
603
604 // save solution in vector x on mesh
605 CHKERR DMoFEMMeshToGlobalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE);
606
607 // Set up post-processor. It is some generic implementation of finite
608 // element
609 auto get_post_proc_ele = [&]() {
610 auto jac_ptr = boost::make_shared<MatrixDouble>();
611 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
612 auto det_ptr = boost::make_shared<VectorDouble>();
613
614 auto post_proc_ele = boost::make_shared<
616 // Add operators to the elements, starting with some generic
617 constexpr int SPACE_DIM = 3;
618 post_proc_ele->getOpPtrVector().push_back(
619 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
620 post_proc_ele->getOpPtrVector().push_back(
621 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
622 post_proc_ele->getOpPtrVector().push_back(
623 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
624
625 auto u_ptr = boost::make_shared<MatrixDouble>();
626 auto grad_ptr = boost::make_shared<MatrixDouble>();
627 post_proc_ele->getOpPtrVector().push_back(
629 post_proc_ele->getOpPtrVector().push_back(
631 grad_ptr));
633
634 post_proc_ele->getOpPtrVector().push_back(
635
636 new OpPPMap(
637
638 post_proc_ele->getPostProcMesh(), post_proc_ele->getMapGaussPts(),
639 {},
640
641 {{"U", u_ptr}},
642
643 {{"GRAD", grad_ptr}},
644
645 {})
646
647 );
648
649 return post_proc_ele;
650 };
651
652 if (auto post_proc = get_post_proc_ele()) {
653 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
654 post_proc);
655 // write output
656 CHKERR post_proc->writeFile("out.h5m");
657 }
658
659 {
660 if (flg_test == PETSC_TRUE) {
661 const double x_vec_norm_const = 0.4;
662 // Check norm_1 value
663 double norm_check;
664 // Takes maximal element of the vector, that should be maximal
665 // displacement at the end of the bar
666 CHKERR VecNorm(x, NORM_INFINITY, &norm_check);
667 if (std::abs(norm_check - x_vec_norm_const) / x_vec_norm_const >
668 1.e-10) {
669 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
670 "test failed (nrm 2 %6.4e)", norm_check);
671 }
672 }
673 }
674
675 // free memory handled by mofem discrete manager for A, x and f
676 CHKERR MatDestroy(&A);
677 CHKERR VecDestroy(&x);
678 CHKERR VecDestroy(&f);
679
680 // free memory allocated for mofem discrete manager
681 CHKERR DMDestroy(&dm);
682
683 // This is a good reference for the future
684 }
686
687 // finish work cleaning memory, getting statistics, etc
689
690 return 0;
691}
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 nme:nme847.
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 ...
@ 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
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1113
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:546
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set compute operator for KSP solver via sub-matrix and IS.
Definition DMMoFEM.cpp:627
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:576
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition DMMoFEM.cpp:525
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:668
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_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
#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
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
FTensor::Index< 'm', 3 > m
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:118
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 field rank 0, i.e. vector field.
Get values at integration pts for tensor field 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:722
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:261
MoFEMErrorCode defineFiniteElements()
Define finite elements.
Definition Simple.cpp:471
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:191
MoFEMErrorCode buildFiniteElements()
Build finite elements.
Definition Simple.cpp:660
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:800
MoFEMErrorCode buildFields()
Build fields.
Definition Simple.cpp:584
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:575
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
Definition Simple.cpp:551
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:389
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference 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 to volume elements.
int operator()(int, int, int p) const