v0.14.0
COR-7: Mixed formulation for incompressible elasticity

Introduction

In this tutorial, the MoFEM implemenetation of the mixed formulation for incompressible elasticity problem will be presented. Interested readers are encouraged to have a look at MoFEM Architecture, FUN-2: Hierarchical approximation, FUN-0: Hello world, and COR-2: Solving the Poisson equation to have better ideas how MoFEM works and how a specific problem is implemented in the platform. They may also find this tutorial video helpful.

The remainder of the tutorial is organised as follows. The mathematical and finite element formulation is presented in the next section before the input mesh of the L-shape structure which is used as an example is defined. It is followed by explanation of the source code and then how to run and visualise the results. Remarks on choosing approximation order and its consequence will be discussed at the end of the tutorial.

Mathematical and finite element formulation

A compressible elasticity problem can be formally defined as follows assuming there is no body force.

For given loading \( \bf{f} \), traction \( \bf{t} \), and displacement \( \bf{\bar u} \), find \( \bf{u} \) such that

\[ \begin{align} -{\rm{div}} \boldsymbol{\sigma(u)} &= \boldsymbol{f} \quad \text{in} \quad \Omega, \\ \boldsymbol{u} &= \boldsymbol{\bar u} \quad \text{on} \quad \Gamma_u \subset \partial\Omega, \\ \boldsymbol{\sigma n} &= \boldsymbol{t} \quad \text{on} \quad \Gamma_t \subset \partial\Omega. \end{align} \]

The constitutive equation is given as

\[ \begin{align} \boldsymbol{\sigma} = \lambda \rm{tr}\boldsymbol{(\varepsilon)I} + 2\mu\boldsymbol{\varepsilon}, \end{align} \]

where Lame's constants \( \lambda \) and \( \mu \) are related with Young's modulus \( E \) and Poisson's ratio \( \nu \) as

\[ \begin{align} \lambda &= \frac{{\nu E}}{{\left( {1 + \nu } \right)\left( {1 - 2\nu } \right)}}, \\ \mu &= \frac{E}{{2\left( {1 + \nu } \right)}}. \end{align} \]

The above problem of compressible elasticity can be solved using typical finite element procedure. However, in case of \( \nu=0.5 \), \( \lambda \) becomes infinity implying incompressible state. This causes volumetric locking which makes the approximation converge very slowly or even unable to converge. Therefore, an alternative formulation needs to be considered to bypass this difficulty. This can be done by introducing additional independent unknown called hydrostatic pressure \( p \) where \( p = -\lambda u_{k,k} \). Therefore, the following equation should be added to the formal problem definition presented in Eqs. (1)-(3)

\[ \begin{align} -{\rm{div}} \boldsymbol{u} + \dfrac{1}{\lambda}p &= 0 \quad \text{in} \quad \Omega. \end{align} \]

With two independent unknowns, \( u \) and \( p \), the incompressible elasticity problem is discretised as follows

\[ \begin{align} \left[ {\begin{array}{*{20}{c}} {\bf{K}}&{\bf{G}}\\ {{{\bf{G}}^T}}&{\bf{P}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{\bf{ u}}}\\ {{\bf{ p}}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{\bf{ f}}}\\ {\bf{0}} \end{array}} \right\}, \end{align} \]

where

\[ \begin{align} {\bf{K}} &= \int\limits_\Omega {{{\bf{B}}^T}{{\bf{D}}_d}{\bf{B}}d\Omega } \\ {\bf{G}} &= - \int\limits_\Omega {{{\bf{B}}^T}{\bf m}{{\bf{N}}_p}d\Omega }, \\ {\bf{P}} &= - \int\limits_\Omega {{\bf{N}}_p^T\frac{1}{\lambda }{{\bf{N}}_p}d\Omega }, \\ {\bf{ f}} &= \int\limits_\Gamma {{\bf{N}}_u^T{\bf{t}}d\Gamma }. \end{align} \]

In the platform of MoFEM, while the complete implementation of user data operators (UDO) can be found in elasticity_mixed_formulation.cpp and ElasticityMixedFormulation.hpp, the matrices above are implemented in UDO as follows

It is worth noting that, in order to obtain stable results for the problem of incompressible elasticity, the approximation order of displacement should be higher than that of hydrostatic pressure. The consequences of choosing inappropriate approximation order of the two unknowns will be discussed later in this tutorial.

Input mesh

As an example, an L-shape structure being clamped at the bottom surface and pressurised with unit magnitude at the top right surface will be used for the analysis of incompressible elasticity. The 3D model and mesh which shown in Figure 1 are created by Cubit using a simple journal script below

reset
set duplicate block elements on
brick x 1 y 2 z 0.5
brick x 2 y 1 z 0.5
move curve 23 midpoint location curve 11 include_merged
unite volume all
Sideset 1 surface 3
Sideset 2 surface 12
Sideset 100 curve all
nodeset 101 vertex all
Sideset 102 surface all
block 1 volume all
block 1 name 'MAT_ELASTIC'
block 1 attribute count 2
block 1 attribute index 1 {young_modulus}
block 1 attribute index 2 {poisson_ratio}
create displacement on surface 3 dof 1 dof 2 dof 3 fix 0
create pressure on surface 12 magnitude 1
volume all scheme tetmesh
volume all size auto factor 7
mesh volume all
save as "/Users/username/mofem_install/um/build/basic_finite_elements/elasticity_mixed_formulation/LShape_incompressible.cub" overwrite

Figure 1: Element mesh.

Code dissection

Initialisation

  • There are two header files in which the first one for basic finite elements and the second one contains the implementation of operators for finite elements.

Solving the problem

  • Stiffness matrix, vector of DOFs, and right-hand-side (RHS) vector are declared and associated to the Discrete Manager
    Mat Aij; // Stiffness matrix
    Vec d, F_ext; // Vector of DOFs and the RHS
    {
    CHKERR VecZeroEntries(d);
    CHKERR VecDuplicate(d, &F_ext);
    CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
    CHKERR MatZeroEntries(Aij);
    }
  • Finite element instances are pointed to the matrix and vector
    // Assign global matrix/vector
    feLhs->ksp_B = Aij;
    feLhs->ksp_f = F_ext;
  • Loop over element to compute and assemble stiffness matrix, vectors of DOFs and RHS
    boost::shared_ptr<DirichletDisplacementBc> dirichlet_bc_ptr(
    new DirichletDisplacementBc(m_field, "U", Aij, d, F_ext));
    dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
    dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
    CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
    CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", feLhs);
    // Assemble pressure and traction forces.
    boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
    F_ext, "U");
    {
    boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
    neumann_forces.begin();
    for (; mit != neumann_forces.end(); mit++) {
    CHKERR DMoFEMLoopFiniteElements(dm, mit->first.c_str(),
    &mit->second->getLoopFe());
    }
    }
    // Assemble forces applied to nodes, see implementation in NodalForce
    boost::ptr_map<std::string, NodalForce> nodal_forces;
    CHKERR MetaNodalForces::setOperators(m_field, nodal_forces, F_ext, "U");
    {
    boost::ptr_map<std::string, NodalForce>::iterator fit =
    nodal_forces.begin();
    for (; fit != nodal_forces.end(); fit++) {
    CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
    &fit->second->getLoopFe());
    }
    }
    // Assemble edge forces
    boost::ptr_map<std::string, EdgeForce> edge_forces;
    CHKERR MetaEdgeForces::setOperators(m_field, edge_forces, F_ext, "U");
    {
    boost::ptr_map<std::string, EdgeForce>::iterator fit =
    edge_forces.begin();
    for (; fit != edge_forces.end(); fit++) {
    CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
    &fit->second->getLoopFe());
    }
    }
    CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
    CHKERR DMMoFEMKSPSetComputeOperators(dm, "ELASTIC", feLhs, nullFE, nullFE);
    CHKERR VecAssemblyBegin(F_ext);
    CHKERR VecAssemblyEnd(F_ext);
  • The matrix and vectors have been assembled and now the system of equation is solved using KSP solver
    KSP solver;
    CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
    CHKERR KSPSetFromOptions(solver);
    CHKERR KSPSetOperators(solver, Aij, Aij);
    CHKERR KSPSetUp(solver);
    CHKERR KSPSolve(solver, F_ext, d);
  • Map the solutions to mesh entities
    CHKERR DMoFEMMeshToGlobalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);

Postprocessing

  • The solutions are passed to the postprocessor and information is written to the output file
    CHKERR post_proc.generateReferenceElementMesh();
    CHKERR post_proc.addFieldValuesPostProc("U");
    CHKERR post_proc.addFieldValuesPostProc("P");
    CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
    PetscPrintf(PETSC_COMM_WORLD, "Output file: %s\n", "out.h5m");
    CHKERR post_proc.postProcMesh.write_file("out.h5m", "MOAB",
    "PARALLEL=WRITE_PART");
  • Finally, all dynamic objects Aij, d, F_ext, and Discrete Manager are destroyed to release memory
    CHKERR MatDestroy(&Aij);
    CHKERR VecDestroy(&d);
    CHKERR VecDestroy(&F_ext);
    CHKERR DMDestroy(&dm);

Running the program

In order to run the program, one should first go to the directory where the problem is located, compile the code and then run the executable file. Typically, this can be done as follows

cd mofem_install/um/build/basic_finite_elements/elasticity_mixed_formulation
make -j2
mpirun -np 2 ./elasticity_mixed_formulation -my_file LShape_incompressible.cub -my_order_p 2 -my_order_u 3 -ksp_type gmres -pc_type lu -pc_factor_mat_solver_type mumps -ksp_monitor

The options are explained as follows

  • mpirun -np 2: two processors are used to run the analysis
  • ./elasticity_mixed_formulation: path and name of the executable file, dot means current directory
  • -my_file LShape_incompressible.cub: name of the mesh file
  • -my_order_p 2: second order approximation for hydrostatic pressure p
  • -my_order_u 3: third order approximation for displacement u
  • -ksp_type gmres: Generalised Minimal Residual is used as an iterative Krylov subspace method
  • -pc_type lu: Direct solver based on LU factorisation as a preconditioner
  • -pc_factor_mat_solver_type mumps: MUMPS is used as matrix solver
  • -ksp_monitor: function to be called at every iteration to monitor the residual/error

Output dissection

Initialisation and Setup

  • Cubit mesh being read
    read cubit meshset 12682136550675316776 type BLOCKSET MAT_ELASTICSET msId 1 name MAT_ELASTIC block header: blockCol = 4294967295 blockMat = 0 blockDimension = 3
    read cubit meshset 12682136550675316777 type NODESET DISPLACEMENTSET msId 1
    read cubit meshset 12682136550675316778 type NODESET msId 101
    read cubit meshset 12682136550675316779 type SIDESET msId 1
    read cubit meshset 12682136550675316780 type SIDESET msId 2
    read cubit meshset 12682136550675316781 type SIDESET PRESSURESET msId 3
    read cubit meshset 12682136550675316782 type SIDESET msId 100
    read cubit meshset 12682136550675316783 type SIDESET msId 102
  • MoFEM version and commit ID
    MoFEM version 0.8.15 (MOAB 5.0.2 Petsc Release Version 3.9.3, Jul, 02, 2018 )
    git commit id cdfa37ee13d92b014cc527b5ab7b844baf73f98d
  • Boundary conditions from the mesh
    D i s p l a c e m e n t
    Flag for X-Translation (0/1): 1
    Flag for Y-Translation (0/1): 1
    Flag for Z-Translation (0/1): 1
    Flag for X-Rotation (0/1): 0
    Flag for Y-Rotation (0/1): 0
    Flag for Z-Rotation (0/1): 0
    Displacement magnitude (X-Translation): 0
    Displacement magnitude (Y-Translation): 0
    Displacement magnitude (Z-Translation): 0
    Displacement magnitude (X-Rotation): N/A
    Displacement magnitude (Y-Rotation): N/A
    Displacement magnitude (Z-Rotation): N/A
  • Material used in the mesh
    meshset 12682136550675316776 type BLOCKSET MAT_ELASTICSET msId 1 name MAT_ELASTIC block header: blockCol = 4294967295 blockMat = 0 blockDimension = 3
    Material Properties
    -------------------
    Young's modulus = 1
    Poisson's ratio = 0.5
    Thermal expansion = 0
    User attribute 1 = 0
    User attribute 2 = 0
    User attribute 3 = 0
    User attribute 4 = 0
    User attribute 5 = 0
    User attribute 6 = 0
    User attribute 7 = 0
    MAT_ELATIC msId 1 nb. tets 471
  • Three fields are added to the database. Fields are then built showing the number of active degrees of freedom (DOFs) for entities including vertices, edges, triangles (faces).
    add: name MESH_NODE_POSITIONS BitFieldId 1 bit number 1 space H1 approximation base AINSWORTH_LEGENDRE_BASE rank 3 meshset 12682136550675316786
    add: name U BitFieldId 2 bit number 2 space H1 approximation base AINSWORTH_LEGENDRE_BASE rank 3 meshset 12682136550675316787
    add: name P BitFieldId 4 bit number 3 space H1 approximation base
    AINSWORTH_LEGENDRE_BASE rank 1 meshset 12682136550675316788
    Build Field MESH_NODE_POSITIONS (rank 0)
    nb added dofs (vertices) 465 (inactive 0)
    nb added dofs 465 (number of inactive dofs 0)
    Build Field U (rank 0)
    nb added dofs (vertices) 465 (inactive 0)
    nb added dofs (edges) 4482 (inactive 0)
    nb added dofs (triangles) 3192 (inactive 0)
    nb added dofs 8139 (number of inactive dofs 0)
    Build Field P (rank 0)
    nb added dofs (vertices) 155 (inactive 0)
    nb added dofs (edges) 747 (inactive 0)
    nb added dofs 902 (number of inactive dofs 0)
    Nb. dofs 9506
    Build Field MESH_NODE_POSITIONS (rank 1)
    nb added dofs (vertices) 465 (inactive 0)
    nb added dofs 465 (number of inactive dofs 0)
    Build Field U (rank 1)
    nb added dofs (vertices) 465 (inactive 0)
    nb added dofs (edges) 4482 (inactive 0)
    nb added dofs (triangles) 3192 (inactive 0)
    nb added dofs 8139 (number of inactive dofs 0)
    Build Field P (rank 1)
    nb added dofs (vertices) 155 (inactive 0)
    nb added dofs (edges) 747 (inactive 0)
    nb added dofs 902 (number of inactive dofs 0)
    Nb. dofs 9506
  • Three elements are added to the database. Finite elements are then built.
    add finite element: FORCE_FE
    add finite element: PRESSURE_FE
    add finite element: ELASTIC
    Build Finite Elements FORCE_FE
    Build Finite Elements PRESSURE_FE
    Build Finite Elements ELASTIC
    Nb. FEs 485
    Nb. FEs 485
    id 00000000000000000000000000000001 name FORCE_FE f_id_row 00000000000000000000000000000010 f_id_col 00000000000000000000000000000010 BitFEId_data 00000000000000000000000000000011 Nb. FEs 0
    id 00000000000000000000000000000001 name FORCE_FE f_id_row 00000000000000000000000000000010 f_id_col 00000000000000000000000000000010 BitFEId_data 00000000000000000000000000000011 Nb. FEs 0
    id 00000000000000000000000000000010 name PRESSURE_FE f_id_row 00000000000000000000000000000010 f_id_col 00000000000000000000000000000010 BitFEId_data 00000000000000000000000000000011 Nb. FEs 14
    id 00000000000000000000000000000010 name PRESSURE_FE f_id_row 00000000000000000000000000000010 f_id_col 00000000000000000000000000000010 BitFEId_data 00000000000000000000000000000011 Nb. FEs 14
    id 00000000000000000000000000000100 name ELASTIC f_id_row 00000000000000000000000000000110 f_id_col 00000000000000000000000000000110 BitFEId_data 00000000000000000000000000000111 Nb. FEs 471
    id 00000000000000000000000000000100 name ELASTIC f_id_row 00000000000000000000000000000110 f_id_col 00000000000000000000000000000110 BitFEId_data 00000000000000000000000000000111 Nb. FEs 471
    Nb. entFEAdjacencies 13328
    Nb. entFEAdjacencies 13328
  • Problem is added and then built
    add problem: DM_ELASTIC_MIX
    Problem DM_ELASTIC_MIX Nb. rows 9041 Nb. cols 9041
    Problem DM_ELASTIC_MIX Nb. rows 9041 Nb. cols 9041
    Partition problem DM_ELASTIC_MIX
    create_Mat: row lower 0 row upper 4521
    create_Mat: row lower 4521 row upper 9041
    partition_problem: rank = 0 FEs row ghost dofs problem id 00000000000000000000000000000001 FiniteElement id 00000000000000000000000000000111 name DM_ELASTIC_MIX Nb. local dof 4653 nb global row dofs 9041
    partition_problem: rank = 0 FEs col ghost dofs problem id 00000000000000000000000000000001 FiniteElement id 00000000000000000000000000000111 name DM_ELASTIC_MIX Nb. local dof 4653 nb global col dofs 9041
    partition_problem: rank = 1 FEs row ghost dofs problem id 00000000000000000000000000000001 FiniteElement id 00000000000000000000000000000111 name DM_ELASTIC_MIX Nb. local dof 4388 nb global row dofs 9041
    partition_problem: rank = 1 FEs col ghost dofs problem id 00000000000000000000000000000001 FiniteElement id 00000000000000000000000000000111 name DM_ELASTIC_MIX Nb. local dof 4388 nb global col dofs 9041
    problem id 00000000000000000000000000000001 FiniteElement id 00000000000000000000000000000111 name DM_ELASTIC_MIX Nb. elems 244 on proc 0
    problem id 00000000000000000000000000000001 FiniteElement id 00000000000000000000000000000111 name DM_ELASTIC_MIX Nb. elems 241 on proc 1
    partition_ghost_col_dofs: rank = 0 FEs col ghost dofs problem id 00000000000000000000000000000001 FiniteElement id 00000000000000000000000000000111 name DM_ELASTIC_MIX Nb. col ghost dof 221 Nb. local dof 4653
    partition_ghost_row_dofs: rank = 0 FEs row ghost dofs problem id 00000000000000000000000000000001 FiniteElement id 00000000000000000000000000000111 name DM_ELASTIC_MIX Nb. row ghost dof 221 Nb. local dof 4653
    partition_ghost_col_dofs: rank = 1 FEs col ghost dofs problem id 00000000000000000000000000000001 FiniteElement id 00000000000000000000000000000111 name DM_ELASTIC_MIX Nb. col ghost dof 181 Nb. local dof 4388
    partition_ghost_row_dofs: rank = 1 FEs row ghost dofs problem id 00000000000000000000000000000001 FiniteElement id 00000000000000000000000000000111 name DM_ELASTIC_MIX Nb. row ghost dof 181 Nb. local dof 4388

Solution

  • KSP solver output
    0 KSP Residual norm 2.422191940287e+02
    1 KSP Residual norm 3.461359964609e-11
  • Output file
    Output file: out.h5m

Visualisation

Once the solution is obtained, one can convert the output file which has the extension of h5m to a vtk one before opening it in ParaView.

mbconvert out.h5m out.vtk
open out.vtk

The visualisation of the results can then be viewed as shown in Figure 2 and Figure 3 for displacement and hydrostatic pressure, respectively.

Figure 2: Visualisation of displacement result.

Figure 3: Visualisation of hydrostatic pressure result.

Remark

As mentioned earlier, the choices of approximation order for hydrostatic pressure and displacement are important. Practically, the approximation order of displacement should be higher than that of hydrostatic pressure. If the same order is assigned for displacement and hydrostatic pressure, oscillating results are expected as shown in Figure 4.

Figure 4: Output of hydrostatic pressure when second order approximation is chosen for both variables.
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
SIDESET
@ SIDESET
Definition: definitions.h:160
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
surface
Definition: surface.py:1
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
PRESSURESET
@ PRESSURESET
Definition: definitions.h:165
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
ElasticityMixedFormulation.hpp
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
MoFEM::Types::BitFieldId
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:42
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MoFEM::DMMoFEMKSPSetComputeOperators
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:678
DirichletDisplacementBc
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
MoFEM::CoreInterface::add_ents_to_field_by_type
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MetaNodalForces::setOperators
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NodalForce > &nodal_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: NodalForce.hpp:128
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1167
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
NODESET
@ NODESET
Definition: definitions.h:159
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1197
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::CoreInterface::add_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
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MetaNeumannForces::setMomentumFluxOperators
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
Definition: SurfacePressure.cpp:2069
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:172
DISPLACEMENTSET
@ DISPLACEMENTSET
Definition: definitions.h:163
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:197
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:535
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:556
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
MoFEM::DMMoFEMCreateMoFEM
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:114
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
convert.n
n
Definition: convert.py:82
N
const int N
Definition: speed_test.cpp:3
MoFEM::CoreTmp< 0 >::Initialize
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
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::CoreInterface::check_finite_element
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
MoFEM::DMoFEMLoopFiniteElements
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:586
MoFEM::CoreInterface::add_field
virtual MoFEMErrorCode add_field(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_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
EshelbianPlasticity::U
@ U
Definition: EshelbianContact.cpp:197
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
MetaEdgeForces::setOperators
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name, std::string mesh_node_positions="MESH_NODE_POSITIONS")
Set integration point operators.
Definition: EdgeForce.hpp:97