v0.9.2
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.

Introduction

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 Basic design, Basics of hierarchical approximation and 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 problem

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.

Problem.png
Figure 1: Geometry of the problem under consideration: orthogonal parallelepiped with dimensions of 1 m along x and y axis and 4 m along z axis that is simply supported on the bottom face and uniformly distributed normal tractions along the upper face.

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.

Mesh input

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.

CubitMesh1_1.png
Figure 2: Mesh geometry of an orthogonal parallelepiped of dimensions of 1x1x4 along x, y and z axis, respectively.

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.

Block1_1.png
Figure 3: Face block added to the bottom face of the orthogonal parallelepiped that will be used to restrain vertical displacements along that face.

Block3_1.png
Figure 4: Face block added to the top face of the orthogonal parallelepiped that will be used to apply uniformly distributed normal tractions along that face.

Block2_1.png
Figure 5: Vertex block added to one of the corners of the bottom face of the orthogonal parallelepiped that will be used to prescribe the boundary conditions represented by the rollers 2 and 3 presented in Figure 1.

Block4_1.png
Figure 6: Vertex block added to one of the corners of the bottom face of the orthogonal parallelepiped that will be used to prescribe the boundary conditions represented by the roller 1 presented in Figure 1.

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

Code dissection

Header files

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

Initialisation

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

// Initialize MoFEM
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);

and now PETSc objects can now be declared. As discussed in Basic design 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.

Interface.png
Figure 7: Database interaction using interface objects Discrete Manager, m_field and Simple interface.

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
int order = 3; // default approximation order
PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "SimpleElasticProblem", "none");
// Set approximation order
CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order, PETSC_NULL);
// Set testing (used by CTest)
CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test, &flg_test, PETSC_NULL);
CHKERR PetscOptionsEnd();

Accessing simple interface

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.

CHKERR simple_interface->getOptions();
CHKERR simple_interface->loadFile();

Accessing simple interface

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.

CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 3);
CHKERR simple_interface->setFieldOrder("U", order); // to approximate function

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.

Mesh block handling

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();
int id = bit->getMeshsetId();
if (id == FIX_BRICK_FACES) { // brick-faces
CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 2,
fix_faces, true);
Range adj_ents;
CHKERR m_field.get_moab().get_adjacencies(fix_faces, 0, false, adj_ents,
moab::Interface::UNION);
CHKERR m_field.get_moab().get_adjacencies(fix_faces, 1, false, adj_ents,
moab::Interface::UNION);
fix_faces.merge(adj_ents);
} else if (id == FIX_NODES) { // node(s)
CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 0,
fix_nodes, true);
} else if (id == BRICK_PRESSURE_FACES) { // brick pressure faces
CHKERR m_field.get_moab().get_entities_by_dimension(
meshset, 2, pressure_faces, true);
} else if (id ==
FIX_NODES_Y_DIR) { // restrained second node in y direction
CHKERR m_field.get_moab().get_entities_by_dimension(
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();
int id = bit->getMeshsetId();

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

CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 2, fix_faces, true);

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;
CHKERR m_field.get_moab().get_adjacencies(fix_faces, 0, false, adj_ents, moab::Interface::UNION);
CHKERR m_field.get_moab().get_adjacencies(fix_faces, 1, false, adj_ents, moab::Interface::UNION);

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.

Mesh block handling

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,

Build and Assemble processes.

WorkFlow.png
Figure 8: MoFEM work flow.

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

CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 3);
CHKERR simple_interface->setFieldOrder("U", order);

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
CHKERR m_field.add_finite_element("PRESSURE");
CHKERR m_field.modify_finite_element_add_field_col("PRESSURE", "U");

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.

Introduction of user defined finite element

DMMoFEM

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

CHKERR DMMoFEMAddElement(dm, "PRESSURE");

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

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.

Build

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

CHKERR simple_interface->buildFields();
CHKERR simple_interface->buildFiniteElements();

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

CHKERR m_field.add_ents_to_finite_element_by_dim(pressure_faces, 2, "PRESSURE");
CHKERR m_field.build_finite_elements("PRESSURE", &pressure_faces);

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();

Finite element instances

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.

OperatorPush.png
Figure 9: Schematic representation inclusion process of User Defined Operators to finite element data structure.

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 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.

Solving the problem

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
Vec x, f;
//allocate memory handled by MoFEM discrete manager for matrix A
CHKERR DMCreateMatrix(dm, &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
CHKERR VecDuplicate(x, &f);

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);
CHKERR DMoFEMLoopFiniteElements(dm, "PRESSURE", pressure_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
CHKERR KSPSetOperators(solver, A, A);
//set up the solver data structure for the iterative solver
CHKERR KSPSetUp(solver);

and finally the problem is solved

//solve the system of linear equations
CHKERR KSPSolve(solver, f, x);

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);

Post processing

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
CHKERR post_proc.generateReferenceElementMesh();
CHKERR post_proc.addFieldValuesPostProc("U");
CHKERR post_proc.addFieldValuesGradientPostProc("U");
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");

Clean memory

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 MatDestroy(&A);
CHKERR VecDestroy(&x);
CHKERR VecDestroy(&f);
//free memory allocated for mofem discrete manager
CHKERR DMDestroy(&dm);

to free the heap memory and prevent any memory leaks.

User Defined Operators (UDOs)

Element stiffness matrix

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 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.

OpK

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(i, j, k, l) = 0.;
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;
tD(i, j, k, l) *= coefficient;
}

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

int row_side, int col_side,
EntityType row_type, EntityType col_type,
{
// 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;

K.resize(nbRows, nbCols, false);
K.clear();

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 Basics of hierarchical approximation tutorial

{
// get sub-block (3x3) of local stiffens matrix, here represented by second
// order tensor
auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
&m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
&m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
&m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
};
// 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
const double a = t_w * vol;
// 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
t_m(i, k) +=
a * (tD(i, j, k, l) * (t_row_diff_base(j) * t_col_diff_base(l)));
// 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,

t_m(i, k) += a * (tD(i, j, k, l) * (t_row_diff_base(j) * t_col_diff_base(l)));

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

t_row_diff_base(j) * t_col_diff_base(l)

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
*/
// 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.
K = trans(K);
CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
&*K.data().begin(), ADD_VALUES);
}
}

OpPressure

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

OpPressure(const double pressure_val = 1) :
pressureVal(pressure_val)
{
}
//vector used to store force vector for each degree of freedom
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
{
// 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)
t_nf(i) += (w * pressureVal * t_base) * t_normal(i);
//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

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)
t_nf(i) += (w * pressureVal * t_base) * t_normal(i);
//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

t_nf(i) += (w * pressureVal * t_base) * 0.5 * t_normal(i);

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);

ApplyDirichletBc

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

ApplyDirichletBc(const Range &fix_faces, const Range &fix_nodes,
const Range &fix_second_node)
: MoFEM::FEMethod(), fixFaces(fix_faces), fixNodes(fix_nodes),
fixSecondNode(fix_second_node) {
// constructor
}
std::set<int> set_fix_dofs;
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);
}
};

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;
ApplyDirichletBc(const Range &fix_faces, const Range &fix_nodes,
const Range &fix_second_node)
: MoFEM::FEMethod(), fixFaces(fix_faces), fixNodes(fix_nodes),
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

Vec x;
CHKERR VecDuplicate(ksp_f, &x);
CHKERR VecZeroEntries(x);

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);

Running the program

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.

simple_elasticity_result.png
Figure 10: Result visualisation in Paraview
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.

Exercises

  • 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

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.
*/
/* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
using namespace boost::numeric;
using namespace MoFEM;
static char help[] = "-order approximation order\n"
"\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;
OpK(bool symm = true)
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(i, j, k, l) = 0.;
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;
tD(i, j, k, l) *= coefficient;
}
/**
* \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
*/
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
// 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
*/
// get sub-block (3x3) of local stiffens matrix, here represented by second
// order tensor
auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
&m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
&m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
&m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
};
// 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
const double a = t_w * vol;
// 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
t_m(i, k) +=
a * (tD(i, j, k, l) * (t_row_diff_base(j) * t_col_diff_base(l)));
// 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
*/
// 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.
K = trans(K);
CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
&*K.data().begin(), ADD_VALUES);
}
}
};
double pressureVal;
OpPressure(const double pressure_val = 1)
pressureVal(pressure_val) {}
// vector used to store force vector for each degree of freedom
MoFEMErrorCode doWork(int side, EntityType type,
// 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();
// 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
&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)
t_nf(i) += (w * pressureVal * t_base) * t_normal(i);
// 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
CHKERR VecSetValues(getFEMethod()->ksp_f, nb_dofs, &data.getIndices()[0],
&nF[0], ADD_VALUES);
}
};
Range fixFaces, fixNodes, fixSecondNode;
ApplyDirichletBc(const Range &fix_faces, const Range &fix_nodes,
const Range &fix_second_node)
: MoFEM::FEMethod(), fixFaces(fix_faces), fixNodes(fix_nodes),
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;
}
};
int main(int argc, char *argv[]) {
const string default_options = "-ksp_type fgmres \n"
"-pc_type lu \n"
"-pc_factor_mat_solver_package mumps \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
int order = 3; // default approximation order
PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "SimpleElasticProblem",
"none");
// Set approximation order
CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
PETSC_NULL);
// Set testing (used by CTest)
CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
&flg_test, PETSC_NULL);
ierr = PetscOptionsEnd();
Simple *simple_interface = m_field.getInterface<MoFEM::Simple>();
CHKERR simple_interface->getOptions();
CHKERR simple_interface->loadFile();
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();
int id = bit->getMeshsetId();
if (id == FIX_BRICK_FACES) { // brick-faces
CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 2,
fix_faces, true);
Range adj_ents;
CHKERR m_field.get_moab().get_adjacencies(fix_faces, 0, false, adj_ents,
moab::Interface::UNION);
CHKERR m_field.get_moab().get_adjacencies(fix_faces, 1, false, adj_ents,
moab::Interface::UNION);
fix_faces.merge(adj_ents);
} else if (id == FIX_NODES) { // node(s)
CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 0,
fix_nodes, true);
} else if (id == BRICK_PRESSURE_FACES) { // brick pressure faces
CHKERR m_field.get_moab().get_entities_by_dimension(
meshset, 2, pressure_faces, true);
} else if (id ==
FIX_NODES_Y_DIR) { // restrained second node in y direction
CHKERR m_field.get_moab().get_entities_by_dimension(
meshset, 0, fix_second_node, true);
} else {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Unknown blockset");
}
}
3);
CHKERR simple_interface->setFieldOrder("U", order);
CHKERR simple_interface->defineFiniteElements();
// Add pressure element
CHKERR m_field.add_finite_element("PRESSURE");
CHKERR m_field.modify_finite_element_add_field_row("PRESSURE", "U");
CHKERR m_field.modify_finite_element_add_field_col("PRESSURE", "U");
CHKERR m_field.modify_finite_element_add_field_data("PRESSURE", "U");
CHKERR simple_interface->defineProblem();
DM dm;
CHKERR simple_interface->getDM(&dm);
CHKERR DMMoFEMAddElement(dm, "PRESSURE");
CHKERR simple_interface->buildFields();
CHKERR simple_interface->buildFiniteElements();
CHKERR m_field.add_ents_to_finite_element_by_dim(pressure_faces, 2,
"PRESSURE");
CHKERR m_field.build_finite_elements("PRESSURE", &pressure_faces);
CHKERR simple_interface->buildProblem();
// Create elements instances
boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(
boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
// 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);
CHKERR DMMoFEMKSPSetComputeRHS(dm, "PRESSURE", pressure_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
Vec x, f;
// allocate memory handled by MoFEM discrete manager for matrix A
CHKERR DMCreateMatrix(dm, &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
CHKERR VecDuplicate(x, &f);
// 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);
CHKERR DMoFEMLoopFiniteElements(dm, "PRESSURE", pressure_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
CHKERR KSPSetOperators(solver, A, A);
// set up the solver data strucure for the iterative solver
CHKERR KSPSetUp(solver);
// solve the system of linear equations
CHKERR KSPSolve(solver, f, x);
// 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
PostProcVolumeOnRefinedMesh post_proc(m_field);
// Add operators to the elements, starting with some generic
CHKERR post_proc.generateReferenceElementMesh();
CHKERR post_proc.addFieldValuesPostProc("U");
CHKERR post_proc.addFieldValuesGradientPostProc("U");
dm, simple_interface->getDomainFEName().c_str(), &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 MatDestroy(&A);
CHKERR VecDestroy(&x);
CHKERR VecDestroy(&f);
// 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;
}