No Matches

The example shows how to solve the linear elastic problem.

/** \file simple_elasticity.cpp
* \example simple_elasticity.cpp
The example shows how to solve the linear elastic problem.
using namespace boost::numeric;
using namespace MoFEM;
static char help[] = "-order approximation order\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);
// 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);
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
// move to next block of local stiffens matrix
// move to next row base function
// move to next integration weight
* \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
// 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);
// 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
// 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
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
// move to next base function
// move to next integration weight
// 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);
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()) {
if (fixSecondNode.find(dit->get()->getEnt()) != fixSecondNode.end()) {
if (dit->get()->getDofCoeffIdx() == 1) {
if (fixNodes.find(dit->get()->getEnt()) != fixNodes.end()) {
std::vector<int> fix_dofs(set_fix_dofs.size());
std::copy(set_fix_dofs.begin(), set_fix_dofs.end(), fix_dofs.begin());
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,
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_type mumps \n"
"-mat_mumps_icntl_20 0 \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;
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",
// Set approximation order
CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
// 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 {
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,
CHKERR m_field.get_moab().get_adjacencies(fix_faces, 1, false, 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 {
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 simple_interface->buildFields();
CHKERR simple_interface->buildFiniteElements();
CHKERR m_field.add_ents_to_finite_element_by_dim(pressure_faces, 2,
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,
// 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;
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
// 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
// save solution in vector x on mesh
// 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;
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
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()) {
// 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) {
"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;
static Index< 'p', 3 > p
const std::string default_options
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static char help[]
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
static PetscErrorCode ierr
constexpr int SPACE_DIM
Catch errors.
Definition: definitions.h:372
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
Definition: definitions.h:148
Definition: definitions.h:40
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1070
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:450
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:503
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set KSP right hand side evaluation function
Definition: DMMMoFEM.cpp:584
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:533
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:482
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
Definition: DMMMoFEM.cpp:625
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
use with loops to iterate row DOFs
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
auto f
Definition: HenckyOps.hpp:5
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
double w(const double g, const double t)
constexpr AssemblyType A
Definition: plastic.cpp:35
MoFEMErrorCode postProcess()
function is run at the end of loop
Set integration rule to boundary elements.
int operator()(int, int, int) const
const Problem * problemPtr
raw pointer to problem
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure for User Loop Methods on finite elements
auto getFTensor0IntegrationWeight()
Get integration weights.
operator doWork function is executed on FE rows
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
MoFEMErrorCode buildProblem()
Build problem.
Definition: Simple.cpp:709
MoFEMErrorCode defineFiniteElements()
Define finite elements.
Definition: Simple.cpp:492
MoFEMErrorCode buildFiniteElements()
Build finite elements.
Definition: Simple.cpp:683
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:373
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:289
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:780
MoFEMErrorCode buildFields()
Build fields.
Definition: Simple.cpp:596
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:303
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:588
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
Definition: Simple.cpp:564
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:313
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
double yOung
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate B^T D B operator.
double pOisson
bool isDiag
true if this block is on diagonal
FTensor::Ddg< double, 3, 3 > tD
int nbIntegrationPts
number of integration points
int nbRows
number of dofs on rows
MatrixDouble K
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations for give operator.
int nbCols
number if dof on column
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble local entity block matrix.
VectorDouble nF
FTensor::Index< 'i', 3 > i
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Set integration rule.
int operator()(int, int, int) const