v0.14.0 |

COR-6: Solid elasticity

This tutorial's motivation is to show how finite element is implemented using basic blocks, like the Lego block where one can combine them in arbitrarily, restricted only by his or her imagination. That is not short and compact code, where PDE can be in single line, but code which can be modified and changed in a way that was not imagined by developers.

This tutorial presents in detail a number of basic Finite Element processes implemented in MoFEM through a user module implementation that solves a simple linear elastic problem. The aim of this presentation is to facilitate new developers to familiarise themselves with certain practices and formulations regarding finite element processes and mesh information manipulation. Readers are encouraged to read MoFEM Architecture, FUN-2: Hierarchical approximation and FUN-0: Hello world before starting the present tutorial. Initially, generation of a mesh is presented and then the module developed for this tutorial is dissected. Finally, code run and post processing of the results are discussed.

The linear elastic problem modelled here is an orthogonal parallelepiped of dimensions \(1.0\times1.0\times4.0 \,\,\rm{m}\), subjected to a uniformly distributed normal traction \(\bar t\) on one face while its opposite face being simply supported. The geometry of the problem is schematically presented in Figure 1. Boundary conditions are applied on the two faces with both edges equal to \(1\,\,\rm m\). Vertical displacements (along z direction) along the whole area of the bottom face are restricted. To prevent rigid body motion, one of the bottom corners displacements along the x and y direction are restricted by the supports numbered 3 and 2 presented in Figure 1, respectively. Finally, displacements along y direction is restricted at a different corner of the bottom face by the support number 1 to prevent rigid body rotation.

The boundary conditions where chosen to allow the body to deform freely on both horizontal directions. In this way, verification will be possible by checking the match between Poisson ratio recovered by the analysis results with the one used as input as a material parameter. Since only normal tractions are applied on the top face, no bending will occur along the faces parallel to the z axis and therefore no extra support is needed along the bottom edge connecting support 1 and 2.

Particular features associated to mesh input to solve the problem described here will be presented next.

The user module presented here is associated with a specific mesh input. Therefore, mesh generation process is briefly presented focusing only on the steps that are closely linked to the module implementation. There is no intention to provide detail guidelines for using the mesh generator. Mesh generators that can be used to run analyses in MoFEM are Salome, GMsh and Cubit. Mesh features presented here are identical to those found in all three aforementioned software packages.

The mesh presented in Figure 2 corresponds to the geometry of the problem illustrated in Figure 1. However, it should be noted that the mesh generator as well as MoFEM do not require a particular choice of units. Instead, the user is responsible to be consisted with her/his own choice of units for any input parameter.

To prescribe the boundary conditions for the generated mesh presented in Figure 1, a number of blocks must be added to the mesh. To prescribe the restriction on the vertical displacements on the bottom face and the normal traction on the upper face, surface blocks 1 and 3 are added on the corresponding faces as presented in Figure 3 and Figure 4, respectively. Moreover, to assign the displacement restrictions represented by rollers 1, 2 and 3 in Figure 1, vertex blocks have to be added as presented in Figure 5 and Figure 6.

No specific boundary conditions are applied through the mesh generator. To achieve prescription of boundary conditions, certain mesh regions have to be numbered and characterised as ''Blocks''. As explained earlier, block feature exists in Salome, GMsh and Cubit. Therefore, this feature is not restricted to the particular choice of mesh generator. MoFEM user module executable will read the mesh input file and apply boundary conditions (see Sections Mesh block handling and ApplyDirichletBc). As an example, a simple journal script which can be used in Cubit to generate the mesh is found below

reset # Geometry brick x 1 y 1 z 4 # Mesh preparation and generation volume 1 scheme tetmesh volume 1 size auto factor 10 mesh volume 1 block 1 surface 2 block 2 vertex 7 block 3 surface 1 block 4 vertex 6 # Boundary condition # The boundary condition for this problem is defined in the source code # Save the model to *.cub file which will be used in MoFEM code save as "Users/username/mofem_install/um/build/basic_finite_elements/simple_elasticity/cubit.cub" overwrite

The header file included in the user module are presented below.

#include <BasicFiniteElements.hpp>

Initially, MoFEM and PETSc has to be initialised through the code line below

// Initialize MoFEM

and now PETSc objects can now be declared. As discussed in MoFEM Architecture article, PETSc objects are used to solve the large computational problem of the system of linear equations.

To create the mesh database and an interface to interact with it, the following piece of code is included

// Create MoFEM database and link it to MoAB

moab::Core mb_instance;//create database

moab::Interface &moab = mb_instance; //create interface to database

MoFEM database is created along with an interface object that will be used as a mean to interact with MoAB database

MoFEM::Core core(moab); //MoFEM database

MoFEM::Interface &m_field = core; //create interface to database

This interaction is schematically presented as the interaction of level 1 and 2 in Figure 7 on the right hand side.

Moreover, the Discrete Manager that is the PETSc implementation within MoFEM has to be register to PETSc (schematically presented as interaction of levels 1 and 2 on the left hand side in Figure 7)

Interaction of levels 2 and 3 is going to be described in Section Accessing simple interface.

To read data from the command line when executing the program the lines below are introduced

// Get command line options

PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error

// Set approximation order

// Set testing (used by CTest)

CHKERR PetscOptionsEnd();

A pointer to an object of type *simple_interface* is declared.

Simple *simple_interface;

CHKERR m_field.getInterface(simple_interface);

The second line of code creates a link between *simple_interface* and *m_field* interface (schematically presented in Figure 7 as the interaction of levels 2 and 3 by the left hand side arrow).

Interface *simple_interface* has a number of functions implemented that allows the user to build and manipulate PETSc, MoAB core and MoFEM core databases. For instance, using the code below the mesh file is loaded and information of the mesh is passed to MoFEM database.

Fields of user's interest can be introduced to the MoFEM database through *simple_interface* object as presented below for the shape functions field in our 3D domain.

By using MoFEM::Simple::addDomainField function the name of the field chosen by the user can be introduced (in this case the name chosen is U). The space chosen to operate is H1 and the type of shape functions is chosen through the key-word AINSWORTH_LEGENDRE_BASE are Legendre polynomials proposed by Ainsworth and Coyle. Finally, the hard-coded number 3 determines the dimensions of the degrees of freedom related to each shape functions. Here, this number is 3 since it is a displacement field in a 3D space. Should the unknown field be pressure, then the hard-coded number will be 1 instead of 3 since pressure always has one dimension. Furthermore, the order of the approximation field is introduced through function MoFEM::Simple::setFieldOrder. As explained earlier, the order of the approximated field can be chosen and the appropriate Hierarchical shape functions will be taken into account accordingly.

Other types of fields can also be introduced to MoFEM database. For instance, an error field can be introduced to evaluate the error variance of the solution throughout the domain. There will be explicit presentation of introduction of other fields in next tutorials.

- Note
- MoFEM has functionality which allow to set boundary conditions directly in in mesh preprocessor, however here we focus on simple where boundary conditions are handled in direct low level implementation.

The blocks numbered 1, 2, 3 and 4 (defined in Section Mesh input) will now be used to apply boundary conditions. First, four Range objects (that are part of MoAB library) are declared for each block

Range fix_faces, pressure_faces, fix_nodes, fix_second_node;

By using the MoFEM core interface m_field, we loop through all meshSets in our MoFEM m_field through the macro _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit) as

enum MyBcTypes {

FIX_BRICK_FACES = 1,

FIX_NODES = 2,

BRICK_PRESSURE_FACES = 3,

FIX_NODES_Y_DIR = 4

};

EntityHandle meshset = bit->getMeshset();

if (id == FIX_BRICK_FACES) { // brick-faces

fix_faces, true);

Range adj_ents;

moab::Interface::UNION);

moab::Interface::UNION);

fix_faces.merge(adj_ents);

} else if (id == FIX_NODES) { // node(s)

fix_nodes, true);

} else if (id == BRICK_PRESSURE_FACES) { // brick pressure faces

meshset, 2, pressure_faces, true);

} else if (id ==

FIX_NODES_Y_DIR) { // restrained second node in y direction

meshset, 0, fix_second_node, true);

} else {

SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Unknown blockset");

}

}

the identity of each meshSet is found as

EntityHandle meshset = bit->getMeshset();

the *id* values correspond to the block numbering presented in Section Mesh input. When the value of id is equal to 1, 2, 3 or 4 boundary conditions are prescribed following the schematic in Figure 1.

For *id* equal to 1, the *meshset* is copied to object *fix_faces* through the MoAB function

were the hard coded number 2, refers to the dimension of entity under consideration which in this case is a face. Numbers 0, 1 and 3 correspond to vertices, edges and volumes, respectively.

Since we are interested in prescribing all degrees of freedom lying along the faces, we need to copy all corresponding entities i.e. vertices and edges, to the fix_faces object. This can be achieved by using the MoAB function *get_adjacencies* that provides information for the entities adjacent to the input meshset (in this case *fix_faces*). To implement this, another *Range* object with name *adj_ents* is created where all the information of the adjacent entities are going to be copied as presented below.

Range adj_ents;

Similar to the case of function *get_entities_by_dimension*, hard coded integers indicate the type of entities under consideration. Thereafter, all information passed to the new *Range* object adj_ents is copied to *fix_faces* object through the MoAB function merge as

fix_faces.merge(adj_ents);

For the vertex type *meshsets*, no adjacencies are needed since only the degrees of freedom associated to each block is going to be prescribed. Therefore, invoking *get_entities_by_dimension* function is sufficient.

For the case of the block of the upper face (*id* = 3), only the degrees of freedom associated to the face are needed. Since the traction is uniform and normal to the upper surface, it can be implemented as pressure and therefore its implementation requires only degrees of freedom associated to the face as presented in section `‘Application of boundary conditions’'.

Finally, the rest of blocks are of vertex type and there is no need to search for adjacencies since only the degrees of freedom at those points are needed to be prescribed.

All problems tackled by user modules proceed with the same workflow for construction of the FEM solution. The workflow implemented in MoFEM is schematically presented in Figure 8 and consists of Definition,

Definition process encompasses creation of stencils of finite elements, fields and problems that are going to be used for the problem solution. Field stencils include name, dimension and order and are defined as

For the elements, the stencil includes element's name, dimension, type of fields that are going to be approximated in this element and the location of the element that the fields operate (volume, boundary or skeleton). Each kind of element used needs a different kind of stencil to be defined.

CHKERR simple_interface->defineFiniteElements();

For reasons that will be apparent in Section OpPressure, we will apply traction boundary conditions through a new element. Since this element is unknown to the simple interface structure we need to explicitly define the element to be included in *m_field*

// Add pressure element

Initially, the new finite element type with name is added to MoFEM core database through the function *add_finite_element*("PRESSURE"), where the string "PRESSURE" used as input is the name that we chose to call the new element. Then, the row and col field information of the "PRESSURE" element are set to be those of "U" type. Furthermore, the list of fields that are associated with "PRESSURE" element listed in the elements data is then set to be only the "U" type field. In a more complicated case, more than one fields can be added to the data, row and col structures of a finite element. Similarly, definition of each problem stencil consists of problem's name and the stencils of the various elements previously presented and are defined using the code lines below

CHKERR simple_interface->defineProblem();

Before defining a problem, all kinds of elements used have to be previously defined.

The new `‘PRESSURE’' finite element has to be explicitly added to DMMoFEM since it was not created through *simple_interface*. Hence, a DM object is declared and it is set to point to the DMMoFEM create by the interface as

DM dm;

CHKERR simple_interface->getDM(&dm);

and then the `‘PRESSURE’' finite element is going to be added to DMMoFEM as

This is schematically illustrated as the interaction of levels 2 and 3 on the left hand side of Figure 7. Then all information introduced to DMMoFEM manager is partitioned

CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);

This has to be explicitly defined in order to be able to solve the problem in parallel computing where matrix assembly, linear system of equations and mesh is distributed between processors.

Then, fields and finite elements declared through *simple_interface* are built (second stage presented in Figure 8).

First, the dimension and order of fields according to the specific input chosen is added to the field. In this case, the approximation field of Hierarchic shape functions of Ainsworth and Coyle Legendre polynomials type of order 3 and dimension 3 will be set to the field.

Thereafter, finite elements are built i.e. entities of the finite elements are linked to the field components that were previously built and the order of approximation is set to each entity.

CHKERR m_field.add_ents_to_finite_element_by_dim(0, simple_interface->getDim(), simple_interface->getDomainFEName(), true);

CHKERR m_field.build_finite_elements(simple_interface->getDomainFEName());

entities associated with the user defined `‘PRESSURE’' finite element have to be explicitly added to the finite element

Then to build all information regarding the element connectivities, and fields associated to element entities into MoFEM core, the function below is invoked.

CHKERR simple_interface->buildProblem();

MoFEM is designed such that approximation fields definition, FE element declaration, FE element definition and FE implementation is explicitly disentangled. That gives us great flexibility allowing for combining individual parts of the code in a countless number of variations.

Implementation FE element can be done on range of levels of abstractions. Here we present implementation which introduces the concept of *User* *Defined* *Operator*. Using its implementation is not constrained to finite element shape or approximation base of fields. In case of coupled or mix formulation, it enables to break down complex code into smaller easy sub-tasks. Here we only introduce basic concepts of this technology.

To create stiffness matrices and apply boundary conditions we will need to finite element instance and attach *User* *Defined* *Operators* (UDOs). For the stiffness matrix evaluation a new shared pointer *elastic_fe* is instantiated to build elastic finite element instance and a UDO is pushed to the elastic finite elements

boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(new VolumeElementForcesAndSourcesCore(m_field));

elastic_fe->getOpPtrVector().push_back(new OpK());

schematically presented in Figure 9.

It is important to state that each instance of the structure of elastic finite element pointed by shared pointed *elastic_fe* carries with it the operator OpK. Element can have sequence of operators, see for example COR-3: Implementing operators for the Poisson equation.

To apply the traction boundary conditions, shared pointer *pressure_fe* is instantiated and then the operator OpPressure is attached to it to prescribe the traction boundary conditions.

// push operators to elastic_fe

boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(new FaceElementForcesAndSourcesCore(m_field));

pressure_fe->getOpPtrVector().push_back(new OpPressure());

Finally, to prescribe the Dirichlet boundary conditions, the shared pointer fix_dofs_fe is instantiated as a ApplyDirichletBc that is inherited from MoFEM::FEMethod type (which lies within the class that are in the finite elements family within MoFEM) as

boost::shared_ptr<FEMethod> fix_dofs_fe(new ApplyDirichletBc(fix_faces, fix_nodes, fix_second_node));

The Ranges of *fix_faces*, *fix_nodes* and *fix_second_node* are given as input to the constructor of ApplyDirichletBc to provide the mesh information regarding the location of the degrees of freedom that will be prescribed.

Implementation of OpK, OpPressure and ApplyDirichletBc is presented in Section User Defined Operators (UDOs). A more general implementation for application of boundary conditions will be presented in future tutorials.

In this example the KSP solver provided by PETSc is used. Initially, finite elements are pushed to the discrete manager and set KSP operators for assembling the stiffness matrix.

// Set operators for KSP solver

dm, simple_interface->getDomainFEName(), elastic_fe, nullFE, nullFE);

Then, right hand side vector is evaluated and KSP operators are set to handle it.

dm, "PRESSURE", pressure_fe, nullFE, nullFE);

Thereafter, stiffness matrix *A*, vector of degrees of freedom *x* and right hand side vector *f* are declared and associated to the Discrete Manager DM

//initialise matrix A used as the global stiffness matrix

Mat A;

//initialise left hand side vector x and right hand side vector f

//allocate memory handled by MoFEM discrete manager for matrix A

//allocate memory handled by MoFEM discrete manager for vector x

CHKERR DMCreateGlobalVector(dm, &x);

//allocate memory handled by MoFEM discrete manager for vector f of the same size as x

Thereafter, members of finite elements instances are set to point to the aforementioned objects

//precondition matrix A according to fix_dofs_fe and elastic_fe finite elements

elastic_fe->ksp_B = A;

fix_dofs_fe->ksp_B = A;

//precondition the right hand side vector f according to fix_dofs_fe and elastic_fe finite elements

fix_dofs_fe->ksp_f = f;

pressure_fe->ksp_f = f;

Now all finite elements are prepared to be iterated (all operators have been pushed to the elements and their members point to the structural problem A, f and x objects) in order to Assembly the problem, i.e. the third part of the work flow presented in Figure 8. Now, each element created is iterated to compute and assemble the stiffness matrix and left hand and right hand vectors as

CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(), elastic_fe);

Furthermore, the part of the prescribed degrees of freedom are assembled in the left hand side vector through

//This is done because only post processor is implemented in the ApplyDirichletBc struct

CHKERR DMoFEMPostProcessFiniteElements(dm, fix_dofs_fe.get());

Implementation of UDOs associated with *pressure_fe* and *fix_dof_fe* are going to be presented in Section User Defined Operators (UDOs).

At this point matrix and two vectors have been assembled and we can now solve our system of linear equations. Solver is declared first and then connected to the MoFEM core as

//make available a KSP solver

KSP solver;

//make the solver available for parallel computing by determining its MPI communicator

CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);

The solver object operators have to be prepared

//making available running all options available for KSP solver in running command

CHKERR KSPSetFromOptions(solver);

//set A matrix to KSP solver and preconditioner

//set up the solver data structure for the iterative solver

CHKERR KSPSetUp(solver);

and finally the problem is solved

the results are written in left hand side vector x

//make vector x available for parallel computations for visualization context

VecView(x, PETSC_VIEWER_STDOUT_WORLD);

and all values of degrees of freedom in x are mapped onto mesh entities

// save solution in vector x on mesh

CHKERR DMoFEMMeshToGlobalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE);

To post process numerical results, initially a post processor object of type PostProcVolumeOnRefinedMesh has to be created based on MoAB database

// Set up post-processor. It is some generic implementation of finite

// element

PostProcVolumeOnRefinedMesh post_proc(m_field);

Thereafter, information of the problem solution is passed to the post processor as

// Add operators to the elements, strating with some generic

auto u_ptr = boost::make_shared<MatrixDouble>();

auto grad_ptr = boost::make_shared<MatrixDouble>();

post_proc_ele->getOpPtrVector().push_back(

new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));

post_proc_ele->getOpPtrVector().push_back(

grad_ptr));

using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;

post_proc_ele->getOpPtrVector().push_back(

new OpPPMap(

post_proc_ele->getPostProcMesh(), post_proc_ele->getMapGaussPts(),

{},

{{"U", u_ptr}},

{{"GRAD", grad_ptr}},

{})

);

CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName().c_str(), &post_proc);

where initially a reference mesh is created is created, then field values are added and finally all information related to finite elements such ass stain and stresses are added into the post processor.

Finally, all information introduced to the post processor object is written to an output file (in this case named out.h5m) as presented below

//write output

CHKERR post_proc.writeFile("out.h5m");

All dynamic memory objects A, x, f and DM that were created have to be destroyed as presented below

//free memory handled by mofem discrete manager for A, x and f

CHKERR VecDestroy(&x);

//free memory allocated for mofem discrete manager

CHKERR DMDestroy(&dm);

to free the heap memory and prevent any memory leaks.

In structural computational mechanics, the relationship to evaluate the element stiffness matrix \(\mathbf{K}\) is given by

\[ {}^{(r,c)}_eK^{\alpha\beta}_{ik} = \int_{\Omega^e} \varphi_{,j}^\alpha D_{ijkl} \phi_{,k}^\beta \, {d\Omega^e} \label{eq:StiffnessIntegral} \]

where \({\Omega^e}\) is the element volume, \({\mathbf {D}}\) is the elastic material stiffness tensor. The row base funcrions and column base functions are represented by \(\pmb\varphi^\alpha\) and \(\pmb\phi^\beta\), respectively. \(\alpha=1,\dots,N^\textrm{row}\) and \(\beta= 1,\dots,N^\textrm{col}\) are indices of base function on row entities and column entities, respectively. \(N^\textrm{row}\) and \(N^\textrm{col}\) are numbers of base funcions on row entity and collumn entity. The differentation is indicated by subscript and index after the comma, as follows

\[ \varphi^\alpha_{,j} := \frac{\partial \varphi^\alpha}{\partial x_j}. \]

where

\[ i,j = 0,1,2 \]

The integration of the element stiffness matrix is broken down into blocks, utilizing the construction of hierarchical approximation basis. For a hierarchical approximation, base DOFs are grouped on vertices, edges, faces and in the volume of the element, see for details to tutorial FUN-0: Hello world. In this tutorial matrix \({}^{(r,c)}_e\mathbf{K}\) is integrated by running over a combination of block, that is indicated by superscripts \((r,c)\) on the left to the matrix \(\mathbf{K}\), where left-subscripts takes values

\[ r,c := \{\textrm{ver, edge, face, vol} \} \label{eq:setI} \]

The element matrix broked down on block subentities cab be represented as follows

\[ {}_e\mathbf{K} = \left[ \begin{array}{cccc} {}^{\textrm {ver ver}}\mathbf{K} & {}^{\textrm {ver edge}}\mathbf{K} & {}^{\textrm {ver face}}\mathbf{K} & {}^{\textrm {ver vol}}\mathbf{K}\\ & {}^{\textrm {edge edge}}\mathbf{K} & {}^{\textrm {edge face}}\mathbf{K} & {}^{\textrm {edge vol}}\mathbf{K}\\ && {}^{\textrm {face face}}\mathbf{K} & {}^{\textrm {face vol}}\mathbf{K}\\ &&& {}^{\textrm {vol vol}}\mathbf{K}\\ \end{array} \right] \label{eq:StiffStiff} \]

Note that only symmetric part is calculated since the element stiffness matrix is symmetric.

The operator called OpK is implemented as a class inheriting class MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator

Struct's public members are

// Finite element stiffness sub-matrix K

// Elastic stiffness tensor (4th rank tensor with minor and major symmetry)

//Young's modulus

double yOung;

//Poisson's ratio

double pOisson;

The first four lines encompass the declaration of matrices used for the calculation of the elements stiffness matrix as presented in Section Element stiffness matrix.

The next line variables are associated to the material stiffness tensor where D is forth order tensor symmetric on two fisrt and two last indices, *yOung* is the Young's modulus and *pOisson* is Poisson's ratio. OpK constructor is implemented as follows

OpK(bool symm = true) : VolumeElementForcesAndSourcesCore::UserDataOperator("U", "U", OPROWCOL, symm)

{

// Evaluation of the elastic stiffness tensor, D, in the Voigt notation is

// done in the constructor

// hardcoded choice of elastic parameters

pOisson = 0.1;

yOung = 10;

// coefficient used in intermediate calculation

const double coefficient = yOung / ((1 + pOisson) * (1 - 2 * pOisson));

tD(0, 0, 0, 0) = 1 - pOisson;

tD(1, 1, 1, 1) = 1 - pOisson;

tD(2, 2, 2, 2) = 1 - pOisson;

tD(0, 1, 0, 1) = 0.5 * (1 - 2 * pOisson);

tD(0, 2, 0, 2) = 0.5 * (1 - 2 * pOisson);

tD(1, 2, 1, 2) = 0.5 * (1 - 2 * pOisson);

tD(0, 0, 1, 1) = pOisson;

tD(1, 1, 0, 0) = pOisson;

tD(0, 0, 2, 2) = pOisson;

tD(2, 2, 0, 0) = pOisson;

tD(1, 1, 2, 2) = pOisson;

tD(2, 2, 1, 1) = pOisson;

}

where the material stiffness tensor for Hooke material is calculated.

When solving the problem, as presented in Section Solving the problem, all finite elements are looped where this loop invokes the OpK::doWork method of each operator of the finite element. The OpK::doWork method of operator OpK is implemented as follows

MoFEMErrorCode doWork(

int row_side, int col_side,

EntityType row_type, EntityType col_type,

EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)

{

// get number of dofs on row

nbRows = row_data.getIndices().size();

// if no dofs on row, exit that work, nothing to do here

if (!nbRows)

// get number of dofs on column

nbCols = col_data.getIndices().size();

// if no dofs on Columbia, exit nothing to do here

if (!nbCols)

// K_ij matrix will have 3 times the number of degrees of freedom of the

// i-th entity set (nbRows)

// and 3 times the number of degrees of freedom of the j-th entity set

// (nbCols)

K.resize(nbRows, nbCols, false);

K.clear();

// get number of integration points

nbIntegrationPts = getGaussPts().size2();

// check if entity block is on matrix diagonal

if (row_side == col_side && row_type == col_type) {

isDiag = true;

} else {

isDiag = false;

}

// integrate local matrix for entity block

CHKERR iNtegrate(row_data, col_data);

// assemble local matrix

CHKERR aSsemble(row_data, col_data);

}

The element's column and row data structures as described in `‘Hierarchic approximation’' tutorial are passed to this function as *col_type* and *row_type* respectively.

The number of degrees of freedom associated to each structure (which are equal to the length of the column and row) is then passed to variables *nbCols* and *nbRows* and checked that they are non-zero

// get number of dofs on row

nbRows = row_data.getIndices().size();

// if no dofs on row, exit that work, nothing to do here

if (!nbRows)

// get number of dofs on column

nbCols = col_data.getIndices().size();

// if no dofs on columns, exit nothing to do here

if (!nbCols)

after the element's lengths of row and columns are passed to nbRows and nbCols and checked to be non-zero, memory for the element's stiffness matrix K that is allocated;

Next, integration of element block stiffnes matrxi is executed, explained in \(\eqref {eq:StiffnessIntegral}\), and finally local matrix is assmble into global stiffnes matrix as follows

// integrate local matrix for entity block

CHKERR iNtegrate(row_data, col_data);

// assemble local matrix

CHKERR aSsemble(row_data, col_data);

Lets first see the OpK::iNtegrate implementation that performs the two inner loops described in Section Entity approximation functions of FUN-2: Hierarchical approximation tutorial

MoFEMErrorCode iNtegrate(

EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)

{

// get sub-block (3x3) of local stiffens matrix, here represented by second

// order tensor

return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(

};

// get element volume

double vol = getVolume();

// get intergrayion weights

auto t_w = getFTensor0IntegrationWeight();

// get derivatives of base functions on rows

auto t_row_diff_base = row_data.getFTensor1DiffN<3>();

// iterate over integration points

for (int gg = 0; gg != nbIntegrationPts; ++gg) {

// calculate scalar weight times element volume

// iterate over row base functions

for (int rr = 0; rr != nbRows / 3; ++rr) {

// get sub matrix for the row

auto t_m = get_tensor2(K, 3 * rr, 0);

// get derivatives of base functions for columns

auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);

// iterate column base functions

for (int cc = 0; cc != nbCols / 3;++cc) {

// integrate block local stiffens matrix

// move to next column base function

++t_col_diff_base;

// move to next block of local stiffens matrix

++t_m;

}

// move to next row base function

++t_row_diff_base;

}

// move to next integration weight

++t_w;

}

}

The code iterates over integration points, base functions on rows and base functions on columns. The essential part of the code is integration itself,

Note that is equivalent of sub integral (eq:StiffnessIntegral ). Here we group operation in blocks using brackets, that first directives of vase functions on rows and columns are multiplied yielding tensor of rank 2

and then the tenor of rank two is multiplied by the tenor of rank two, yielding tensor of rank two. Resultant is multiplied by integration weight times volume of the element and added to local element stiffnes matrix. For more information about tenorial operations see FTensor library

Finally, block local stiffness matrix is assembled into global stiffness matrix as follows

/**

* \brief Assemble local entity block matrix

* @param row_data row data (consist base functions on row entity)

* @param col_data column data (consist base functions on column entity)

* @return error code

*/

MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data,

EntitiesFieldData::EntData &col_data) {

// get pointer to first global index on row

const int *row_indices = &*row_data.getIndices().data().begin();

// get pointer to first global index on column

const int *col_indices = &*col_data.getIndices().data().begin();

Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B

: getFEMethod()->snes_B;

// assemble local matrix

CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,

&*K.data().begin(), ADD_VALUES);

if (!isDiag && sYmm) {

// if not diagonal term and since global matrix is symmetric assemble

// transpose term.

CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,

&*K.data().begin(), ADD_VALUES);

}

}

The operator created for the "PRESSURE" finite element that is used to evaluate and assign Neumann boundary conditions according to (12) presented in FUN-2: Hierarchical approximation will be presented in detail. The operator is the struct presented below

pressureVal(pressure_val)

{

}

//vector used to store force vector for each degree of freedom

VectorDouble nF;

{

// check that the faces have associated degrees of freedom

const int nb_dofs = data.getIndices().size();

if (nb_dofs == 0)

//size of force vector associated to the entity

//set equal to the number of degrees of freedom of associated with the entity

nF.resize(nb_dofs, false);

nF.clear();

// get number of gauss points

const int nb_gauss_pts = data.getN().size1();

// create a 3d vector to be used as the normal to the face with length equal

// to the face area

auto t_normal = getFTensor1Normal();

//vector of base functions

auto t_base = data.getFTensor0N();

// get intergrayion weights

auto t_w = getFTensor0IntegrationWeight();

// loop over all gauss points of the face

for (int gg = 0; gg != nb_gauss_pts; ++gg) {

// weight of gg gauss point

double w = 0.5 * t_w;

//create a vector t_nf whose pointer points an array of 3 pointers pointing to nF

//memory location of components

FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_nf(&nF[0], &nF[1], &nF[2]);

for (int bb = 0; bb != nb_dofs / 3; ++bb) {

//scale the three components of t_normal and pass them to the t_nf (hence to nF)

//move the pointer to next element of t_nf

++t_nf;

//move to next base function

++t_base;

}

// move to next integration weight

++t_w;

}

getFEMethod()->ksp_f, nb_dofs, &data.getIndices()[0], &nF[0], ADD_VALUES);

}

};

Initially, vector *nF* that will contain the boundary condition forces is declared and index i is defined to be used to perform compact operations later discussed

//vector used to store force vector for each degree of freedom

VectorDouble nF;

The evaluation of boundary conditions for each element that are applied is performed within the OpPressure::doWork function. The number of degrees of freedom is initially passed to *nb_dofs* that is used to appropriately resize vector *nF*.

//size of force vector associated to the entity

//set equal to the number of degrees of freedom of associated with the entity

nF.resize(nb_dofs, false);

nF.clear();

Thereafter, the number of gauss points associated to the triangular element is assigned to variable *nb_gauss_pts*

// get number of gauss points

const int nb_gauss_pts = data.getN().size1();

the normal to the face is passed to *t_normal* variable

// create a 3d vector to be used as the normal to the face with length equal

// to the face area

auto t_normal = getFTensor1Normal();

Vector *t_normal* has 3 component (one for each global axis). Its magnitude is equal to the double of the area of the triangle face since FaceElementForcesAndSourcesCore:UserDataOperator:getTensor1Norma() function returns the vector evaluated by the cross product of the vectors having a common origin on one of the faces vertex and their non-coinciding ends located at the other two vertices. This product is the area defined by the parallelogram rule.

Furthermore, *t_base* vector containing the element's shape functions is initialised

//vector of base functions

auto t_base = data.getFTensor0N();

Then, the loop over all gauss points is performed to evaluate the area integral

// get intergrayion weights

auto t_w = getFTensor0IntegrationWeight();

// loop over all gauss points of the face

for (int gg = 0; gg != nb_gauss_pts; ++gg) {

// weight of gg gauss point

double w = 0.5 * t_w;

//create a vector t_nf whose pointer points an array of 3 pointers pointing to nF

//memory location of components

auto t_nf(&nF[0], &nF[1], &nF[2], 3);

for (int bb = 0; bb != nb_dofs / 3; ++bb) {

//scale the three components of t_normal and pass them to the t_nf (hence to nF)

//move the pointer to next element of t_nf

++t_nf;

//move to next base function

++t_base;

}

// move to next integration weight

++t_w;

}

In the loop, initially the weight of the gauss point is assigned to w

double w = getGaussPts()(2, gg);

then *t_nF* that is a pointer to a vector is initialised in a fashion that points to a piece of memory containing three contiguous double values that are in this case the first three elements of the boundary condition vector *nF*.

Since the *t_nF* presents this structure of three contiguous elements pointed by each pointer value, we only need to loop over a third of the number of degrees of freedom.

for (int bb = 0; bb != nb_dofs / 3; bb++)

and the within the loop for each *t_nF* value, the corresponding block of three members of nF can be evaluated as

Index i in the last line operates the for each part of the three elements of *nF*. So in the first instance, where bb = 0, *nF*[0], *nF*[1] and *nF*[2] are evaluated.

Thereafter, pointers *t_nf* and *t_base* are increased, and now they are pointing to the next corresponding bit of memory.

++t_nf;

++t_base;

In the second instance, where *bb* = 1, *t_nf(i)* operates in the next three elements of *nF*, i.e. *nF*[3], *nF*[4], *nF*[5] are evaluated and so on.

Each member of *nF* contains the surface integral presented in (12) of `‘Hierarchical Approximation’' associated to the the corresponding to gauss point gg and degree of freedom 3*bb. The product that is assigned to *t_nf(i)* involves multiplication of weight w of gauss point gg, the uniformly distributed traction is a scalar pressureVal times the shape function *t_base* of degrees of freedom encapsulated in bb, times half the vector normal to the surface. The halving of the vector is done to recover the face area that is the double of the vector magnitude as describe previously.

Finally, the values passed to vector *nF* are then copied to the right hand side vector \r ksp_f that points to vector f presented in Section Solving the problem used for the solution of the system of linear equations.

ierr = VecSetValues(getFEMethod()->ksp_f, nb_dofs, &data.getIndices()[0], &nF[0], ADD_VALUES);

The class implemented to apply homogeneous Dirichlet boundary conditions named ApplyDirichletBc is presented below

const Range &fix_second_node)

fixSecondNode(fix_second_node) {

// constructor

}

std::set<int> set_fix_dofs;

if (dit->get()->getDofCoeffIdx() == 2) {

set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());

}

}

if (dit->get()->getDofCoeffIdx() == 1) {

set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());

}

}

set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());

}

}

std::vector<int> fix_dofs(set_fix_dofs.size());

std::copy(set_fix_dofs.begin(), set_fix_dofs.end(), fix_dofs.begin());

}

};

Entities that were previously handled in Section Mesh block handling are assigned through the constructor in public variables *fixFaces*, *fixNodes*, *fixSecondNode*

Range fixFaces, fixNodes, fixSecondNode;

const Range &fix_second_node)

fixSecondNode(fix_second_node) {

// constructor

}

The homogeneous boundary conditions are then implemented in postProcess() function. Initially, a loop over degrees of freedom is performed through macro `‘for (\_IT\_NUMEREDDOF\_ROW\_FOR\_LOOP\_(problemPtr, dit))’'. Within the loop, the degrees of freedom corresponding to the entities passed to the public variables of the struct are inserted to vector of integers *set_fix_dofs* through insert function

std::set<int> set_fix_dofs;

for (_IT_NUMEREDDOF_ROW_FOR_LOOP_(problemPtr, dit)) {

if (dit->get()->getDofCoeffIdx() == 2) {

if (fixFaces.find(dit->get()->getEnt()) != fixFaces.end()) {

set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());

}

}

if (fixSecondNode.find(dit->get()->getEnt()) != fixSecondNode.end()) {

if (dit->get()->getDofCoeffIdx() == 1) {

printf("The extra node \n");

set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());

}

}

if (fixNodes.find(dit->get()->getEnt()) != fixNodes.end()) {

set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());

}

}

then this data is pass to a new vector of integers *fix_dofs*

std::vector<int> fix_dofs(set_fix_dofs.size());

std::copy(set_fix_dofs.begin(), set_fix_dofs.end(), fix_dofs.begin());

Thereafter, the stiffness matrix *ksp_B* and right hand side vector ksp_f structures are constructed

CHKERR MatAssemblyBegin(ksp_B, MAT_FINAL_ASSEMBLY);

CHKERR MatAssemblyEnd(ksp_B, MAT_FINAL_ASSEMBLY);

CHKERR VecAssemblyBegin(ksp_f);

CHKERR VecAssemblyEnd(ksp_f);

and the left hand side vector *x* is constructed by copying *ksp_f* to *x* and subsequently zeroing it

Invoking function MatZeroRowsColumns the columns and rows corresponding to the prescribed degrees of freedom stored in fix_dofs and the diagonal element is set to unity and the corresponding degrees of freedom in vector x are set to zero.

CHKERR MatZeroRowsColumns(ksp_B, fix_dofs.size(), &fix_dofs[0], 1, x, ksp_f);

finally the memory allocated for vector *x* is freed

CHKERR VecDestroy(&x);

To run the code one should change directory to

cd $HOME/mofem_installation/users_modules/basic_finite_elements/simple_elasticity

Then run the command

./simple_elasticity -file_name simple_elasticity_part.h5m -my_order 2

To visualise results first convert the out put file from .h5m format to a vtk one

mbconvert out.h5m out.vtk

then open the .vtk file in ParaView http://www.paraview.org. For example on macOS

open out.vtk

Figure 10 shows the ParaView visualisation of the result of the considering problem.

- Note
- In the following example solution is homogenous. Thus linear approximation is sufficient to get the exact the solution. Thus setting second order approximation does not improve solution. Note that all DOFs on edges for this case have zero values.

- Change boundary conditions to have fix cantilever and apply tractions to bend beam in
*x*or*y*direction? Check how displacements would converge with increasing approximation order? - Modify code to apply torsion to the beam.
- Write finite element to calculate strain energy of the beam. With that at hand plot convergence for beam subjected bending and torsion.

If you have problems with understanding this tutorial or exercises, please contact us on Q&A forum, see link.

The plain program is located in *users_modules/basic_finite_elements/simple_elasticity/simple_elasticity.cpp*

/** \file simple_elasticity.cpp

* \example simple_elasticity.cpp

The example shows how to solve the linear elastic problem.

*/

#include <BasicFiniteElements.hpp>

using namespace boost::numeric;

using namespace MoFEM;

"\n";

// Finite element stiffness sub-matrix K_ij

// Elastic stiffness tensor (4th rank tensor with minor and major symmetry)

// Young's modulus

double yOung;

// Poisson's ratio

double pOisson;

symm) {

// Evaluation of the elastic stiffness tensor, D

// hardcoded choice of elastic parameters

pOisson = 0.1;

yOung = 10;

// coefficient used in intermediate calculation

const double coefficient = yOung / ((1 + pOisson) * (1 - 2 * pOisson));

tD(0, 0, 0, 0) = 1 - pOisson;

tD(1, 1, 1, 1) = 1 - pOisson;

tD(2, 2, 2, 2) = 1 - pOisson;

tD(0, 1, 0, 1) = 0.5 * (1 - 2 * pOisson);

tD(0, 2, 0, 2) = 0.5 * (1 - 2 * pOisson);

tD(1, 2, 1, 2) = 0.5 * (1 - 2 * pOisson);

tD(0, 0, 1, 1) = pOisson;

tD(1, 1, 0, 0) = pOisson;

tD(0, 0, 2, 2) = pOisson;

tD(2, 2, 0, 0) = pOisson;

tD(1, 1, 2, 2) = pOisson;

tD(2, 2, 1, 1) = pOisson;

}

/**

* \brief Do calculations for give operator

* @param row_side row side number (local number) of entity on element

* @param col_side column side number (local number) of entity on element

* @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET

* @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET

* @param row_data data for row

* @param col_data data for column

* @return error code

*/

EntityType col_type,

EntitiesFieldData::EntData &row_data,

EntitiesFieldData::EntData &col_data) {

// get number of dofs on row

nbRows = row_data.getIndices().size();

// if no dofs on row, exit that work, nothing to do here

if (!nbRows)

// get number of dofs on column

nbCols = col_data.getIndices().size();

// if no dofs on Columbia, exit nothing to do here

if (!nbCols)

// K_ij matrix will have 3 times the number of degrees of freedom of the

// i-th entity set (nbRows)

// and 3 times the number of degrees of freedom of the j-th entity set

// (nbCols)

K.resize(nbRows, nbCols, false);

K.clear();

// get number of integration points

nbIntegrationPts = getGaussPts().size2();

// check if entity block is on matrix diagonal

if (row_side == col_side && row_type == col_type) {

isDiag = true;

} else {

isDiag = false;

}

// integrate local matrix for entity block

CHKERR iNtegrate(row_data, col_data);

// assemble local matrix

CHKERR aSsemble(row_data, col_data);

}

protected:

int nbRows; ///< number of dofs on rows

int nbCols; ///< number if dof on column

int nbIntegrationPts; ///< number of integration points

bool isDiag; ///< true if this block is on diagonal

/**

* \brief Integrate B^T D B operator

* @param row_data row data (consist base functions on row entity)

* @param col_data column data (consist base functions on column entity)

* @return error code

*/

MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,

EntitiesFieldData::EntData &col_data) {

// get sub-block (3x3) of local stiffens matrix, here represented by second

// order tensor

return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(

};

// get element volume

double vol = getVolume();

// get intergration weights

auto t_w = getFTensor0IntegrationWeight();

// get derivatives of base functions on rows

auto t_row_diff_base = row_data.getFTensor1DiffN<3>();

// iterate over integration points

for (int gg = 0; gg != nbIntegrationPts; ++gg) {

// calculate scalar weight times element volume

// iterate over row base functions

for (int rr = 0; rr != nbRows / 3; ++rr) {

// get sub matrix for the row

auto t_m = get_tensor2(K, 3 * rr, 0);

// get derivatives of base functions for columns

auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);

// iterate column base functions

for (int cc = 0; cc != nbCols / 3; ++cc) {

// integrate block local stiffens matrix

// move to next column base function

++t_col_diff_base;

// move to next block of local stiffens matrix

++t_m;

}

// move to next row base function

++t_row_diff_base;

}

// move to next integration weight

++t_w;

}

}

/**

* \brief Assemble local entity block matrix

* @param row_data row data (consist base functions on row entity)

* @param col_data column data (consist base functions on column entity)

* @return error code

*/

MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data,

EntitiesFieldData::EntData &col_data) {

// get pointer to first global index on row

// get pointer to first global index on column

Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B

: getFEMethod()->snes_B;

// assemble local matrix

CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,

&*K.data().begin(), ADD_VALUES);

if (!isDiag && sYmm) {

// if not diagonal term and since global matrix is symmetric assemble

// transpose term.

CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,

&*K.data().begin(), ADD_VALUES);

}

}

};

double pressureVal;

pressureVal(pressure_val) {}

// vector used to store force vector for each degree of freedom

VectorDouble nF;

EntitiesFieldData::EntData &data) {

// check that the faces have associated degrees of freedom

if (nb_dofs == 0)

// size of force vector associated to the entity

// set equal to the number of degrees of freedom of associated with the

// entity

nF.resize(nb_dofs, false);

nF.clear();

// get number of gauss points

// create a 3d vector to be used as the normal to the face with length equal

// to the face area

auto t_normal = getFTensor1Normal();

// get intergration weights

auto t_w = getFTensor0IntegrationWeight();

// vector of base functions

auto t_base = data.getFTensor0N();

// loop over all gauss points of the face

for (int gg = 0; gg != nb_gauss_pts; ++gg) {

// weight of gg gauss point

double w = 0.5 * t_w;

// create a vector t_nf whose pointer points an array of 3 pointers

// pointing to nF memory location of components

FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_nf(&nF[0], &nF[1],

&nF[2]);

for (int bb = 0; bb != nb_dofs / 3; ++bb) {

// scale the three components of t_normal and pass them to the t_nf

// (hence to nF)

// move the pointer to next element of t_nf

++t_nf;

// move to next base function

++t_base;

}

// move to next integration weight

++t_w;

}

// add computed values of pressure in the global right hand side vector

&nF[0], ADD_VALUES);

}

};

Range fixFaces, fixNodes, fixSecondNode;

const Range &fix_second_node)

fixSecondNode(fix_second_node) {

// constructor

}

MoFEMErrorCode postProcess() {

std::set<int> set_fix_dofs;

for (_IT_NUMEREDDOF_ROW_FOR_LOOP_(problemPtr, dit)) {

if (dit->get()->getDofCoeffIdx() == 2) {

if (fixFaces.find(dit->get()->getEnt()) != fixFaces.end()) {

set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());

}

}

if (fixSecondNode.find(dit->get()->getEnt()) != fixSecondNode.end()) {

if (dit->get()->getDofCoeffIdx() == 1) {

set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());

}

}

if (fixNodes.find(dit->get()->getEnt()) != fixNodes.end()) {

set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());

}

}

std::vector<int> fix_dofs(set_fix_dofs.size());

std::copy(set_fix_dofs.begin(), set_fix_dofs.end(), fix_dofs.begin());

CHKERR MatAssemblyBegin(ksp_B, MAT_FINAL_ASSEMBLY);

CHKERR MatAssemblyEnd(ksp_B, MAT_FINAL_ASSEMBLY);

CHKERR VecAssemblyBegin(ksp_f);

CHKERR VecAssemblyEnd(ksp_f);

Vec x;

CHKERR VecDuplicate(ksp_f, &x);

CHKERR VecZeroEntries(x);

CHKERR MatZeroRowsColumns(ksp_B, fix_dofs.size(), &fix_dofs[0], 1, x,

ksp_f);

CHKERR VecDestroy(&x);

}

};

/**

* \brief Set integration rule to volume elements

*

* This rule is used to integrate \f$\nabla v \cdot \nabla u\f$, thus

* if the approximation field and the testing field are polynomials of order

* "p", then the rule for the exact integration is 2*(p-1).

*

* Integration rule is order of polynomial which is calculated exactly. Finite

* element selects integration method based on return of this function.

*

*/

struct VolRule {

int operator()(int, int, int p) const { return 2 * (p - 1); }

};

/**

* \brief Set integration rule to boundary elements

*

* This rule is used to integrate the work of external forces on a face,

* i.e. \f$f \cdot v\f$, where f is the traction vector and v is the test

* vector function. The current problem features a Neumann bc with

* a pre-defined constant pressure. Therefore, if the test field is

* represented by polynomials of order "p", then the rule for the exact

* integration is also p.

*

* Integration rule is order of polynomial which is calculated exactly. Finite

* element selects integration method based on return of this function.

*

*/

struct FaceRule {

int operator()(int, int, int p) const { return p; }

};

const string default_options = "-ksp_type fgmres \n"

"-pc_type lu \n"

"-pc_factor_mat_solver_type mumps \n"

"-mat_mumps_icntl_20 0 \n"

"-ksp_monitor\n";

string param_file = "param_file.petsc";

if (!static_cast<bool>(ifstream(param_file))) {

std::ofstream file(param_file.c_str(), std::ios::ate);

if (file.is_open()) {

file << default_options;

file.close();

}

}

MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);

// Create mesh database

moab::Core mb_instance; // create database

moab::Interface &moab = mb_instance; // create interface to database

try {

// Create MoFEM database and link it to MoAB

MoFEM::Core core(moab);

MoFEM::Interface &m_field = core;

// Get command line options

PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error

"none");

// Set approximation order

PETSC_NULL);

// Set testing (used by CTest)

&flg_test, PETSC_NULL);

ierr = PetscOptionsEnd();

CHKERR simple_interface->getOptions();

Range fix_faces, pressure_faces, fix_nodes, fix_second_node;

enum MyBcTypes {

FIX_BRICK_FACES = 1,

FIX_NODES = 2,

BRICK_PRESSURE_FACES = 3,

FIX_NODES_Y_DIR = 4

};

EntityHandle meshset = bit->getMeshset();

if (id == FIX_BRICK_FACES) { // brick-faces

fix_faces, true);

Range adj_ents;

moab::Interface::UNION);

moab::Interface::UNION);

fix_faces.merge(adj_ents);

} else if (id == FIX_NODES) { // node(s)

fix_nodes, true);

} else if (id == BRICK_PRESSURE_FACES) { // brick pressure faces

meshset, 2, pressure_faces, true);

} else if (id ==

FIX_NODES_Y_DIR) { // restrained second node in y direction

meshset, 0, fix_second_node, true);

} else {

SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Unknown blockset");

}

}

3);

CHKERR simple_interface->defineFiniteElements();

// Add pressure element

CHKERR simple_interface->defineProblem();

DM dm;

CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);

CHKERR simple_interface->buildFields();

CHKERR simple_interface->buildFiniteElements();

CHKERR m_field.add_ents_to_finite_element_by_dim(pressure_faces, 2,

"PRESSURE");

CHKERR simple_interface->buildProblem();

// Create elements instances

boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(

new VolumeElementForcesAndSourcesCore(m_field));

boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(

new FaceElementForcesAndSourcesCore(m_field));

// Set integration rule to elements instances

elastic_fe->getRuleHook = VolRule();

pressure_fe->getRuleHook = FaceRule();

// Add operators to element instances

// push operators to elastic_fe

elastic_fe->getOpPtrVector().push_back(new OpK());

// push operators to pressure_fe

pressure_fe->getOpPtrVector().push_back(new OpPressure());

boost::shared_ptr<FEMethod> fix_dofs_fe(

new ApplyDirichletBc(fix_faces, fix_nodes, fix_second_node));

boost::shared_ptr<FEMethod> null_fe;

// Set operators for KSP solver

dm, simple_interface->getDomainFEName(), elastic_fe, null_fe, null_fe);

null_fe);

// initialise matrix A used as the global stiffness matrix

Mat A;

// initialise left hand side vector x and right hand side vector f

// allocate memory handled by MoFEM discrete manager for matrix A

// allocate memory handled by MoFEM discrete manager for vector x

CHKERR DMCreateGlobalVector(dm, &x);

// allocate memory handled by MoFEM discrete manager for vector f of the

// same size as x

// precondition matrix A according to fix_dofs_fe and elastic_fe finite

// elements

elastic_fe->ksp_B = A;

fix_dofs_fe->ksp_B = A;

// precondition the right hand side vector f according to fix_dofs_fe and

// elastic_fe finite elements

fix_dofs_fe->ksp_f = f;

pressure_fe->ksp_f = f;

elastic_fe);

// This is done because only post processor is implemented in the

// ApplyDirichletBc struct

CHKERR DMoFEMPostProcessFiniteElements(dm, fix_dofs_fe.get());

// make available a KSP solver

KSP solver;

// make the solver available for parallel computing by determining its MPI

// communicator

CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);

// making available running all options available for KSP solver in running

// command

CHKERR KSPSetFromOptions(solver);

// set A matrix with preconditioner

// set up the solver data strucure for the iterative solver

CHKERR KSPSetUp(solver);

// solve the system of linear equations

// destroy solver no needed any more

CHKERR KSPDestroy(&solver);

// make vector x available for parallel computations for visualization

// context

VecView(x, PETSC_VIEWER_STDOUT_WORLD);

// save solution in vector x on mesh

CHKERR DMoFEMMeshToGlobalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE);

// Set up post-processor. It is some generic implementation of finite

// element

auto get_post_proc_ele = [&]() {

auto jac_ptr = boost::make_shared<MatrixDouble>();

auto inv_jac_ptr = boost::make_shared<MatrixDouble>();

auto det_ptr = boost::make_shared<VectorDouble>();

auto post_proc_ele = boost::make_shared<

// Add operators to the elements, starting with some generic

constexpr int SPACE_DIM = 3;

post_proc_ele->getOpPtrVector().push_back(

new OpCalculateHOJac<SPACE_DIM>(jac_ptr));

post_proc_ele->getOpPtrVector().push_back(

new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));

post_proc_ele->getOpPtrVector().push_back(

auto u_ptr = boost::make_shared<MatrixDouble>();

auto grad_ptr = boost::make_shared<MatrixDouble>();

post_proc_ele->getOpPtrVector().push_back(

post_proc_ele->getOpPtrVector().push_back(

grad_ptr));

post_proc_ele->getOpPtrVector().push_back(

new OpPPMap(

post_proc_ele->getPostProcMesh(), post_proc_ele->getMapGaussPts(),

{},

{{"U", u_ptr}},

{{"GRAD", grad_ptr}},

{})

);

return post_proc_ele;

};

if (auto post_proc = get_post_proc_ele()) {

post_proc);

// write output

CHKERR post_proc->writeFile("out.h5m");

}

{

if (flg_test == PETSC_TRUE) {

const double x_vec_norm_const = 0.4;

// Check norm_1 value

double norm_check;

// Takes maximal element of the vector, that should be maximal

// displacement at the end of the bar

CHKERR VecNorm(x, NORM_INFINITY, &norm_check);

if (std::abs(norm_check - x_vec_norm_const) / x_vec_norm_const >

1.e-10) {

SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,

"test failed (nrm 2 %6.4e)", norm_check);

}

}

}

// free memory handled by mofem discrete manager for A, x and f

CHKERR VecDestroy(&x);

// free memory allocated for mofem discrete manager

CHKERR DMDestroy(&dm);

// This is a good reference for the future

}

// finish work cleaning memory, getting statistics, etc

return 0;

}

#define MoFEMFunctionReturnHot(a)

Last executable line of each PETSc function used for error handling. Replaces return()

MoFEMErrorCode getInterface(IFACE *&iface) const

Get interface reference to pointer of interface.

FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradSymTensorGrad< BASE_DIM, SPACE_DIM, SPACE_DIM, 0 > OpK

[Define entities]

Data on single entity (This is passed as argument to DataOperator::doWork)

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

MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)

Assemble PETSc matrix.

double w(double eqiv, double dot_tau, double f, double sigma_y, double sigma_Y)

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

Get values at integration pts for tensor filed rank 1, i.e. vector field.

#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)

Iterator that loops over a specific Cubit MeshSet in a moFEM field.

MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)

Load mesh file.

static MoFEMErrorCode Finalize()

Checks for options to be called at the conclusion of the program.

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.

MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)

Assemble PETSc vector.

PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)

add element to dm

Deprecated interface functions.

ApplyDirichletBc(const Range &fix_faces, const Range &fix_nodes, const Range &fix_second_node)

FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)

Get base function as Tensor0.

virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0

add finite element

virtual moab::Interface & get_moab()=0

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 build_finite_elements(int verb=DEFAULT_VERBOSITY)=0

Build finite elements.

default operator for TRI element

Set inverse jacobian to base functions.

PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)

set ghosted vector values on all existing mesh entities

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.

const VectorInt & getIndices() const

Get global indices of dofs on entity.

PetscErrorCode DMRegister_MoFEM(const char sname[])

Register MoFEM problem.

Face finite element.

PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)

execute finite element method for each element in dm (problem)

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

Volume finite element base.

MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)

Set field order.

Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.

#define _IT_NUMEREDDOF_ROW_FOR_LOOP_(PROBLEMPTR, IT)

use with loops to iterate row DOFs

static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])

Initializes the MoFEM database PETSc, MOAB and MPI.

FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)

Get derivatives of base functions.

MatrixDouble & getN(const FieldApproximationBase base)

get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....

ForcesAndSourcesCore::UserDataOperator UserDataOperator

static MoFEMErrorCodeGeneric< PetscErrorCode > ierr

virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0

set finite element field data

#define MoFEMFunctionBeginHot

First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...

MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)

define problem

PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())

Executes FEMethod for finite elements in DM.

#define MoFEMFunctionReturn(a)

Last executable line of each PETSc function used for error handling. Replaces return()

#define MoFEMFunctionBegin

First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...

PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)

Post post-proc data at points from hash maps.

Generated by Doxygen 1.8.17 and hosted at