v0.14.0
Loading...
Searching...
No Matches
simple_elasticity.cpp

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"
"\n";
struct OpK : public VolumeElementForcesAndSourcesCore::UserDataOperator {
// Finite element stiffness sub-matrix K_ij
MatrixDouble K;
// 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)
: VolumeElementForcesAndSourcesCore::UserDataOperator("U", "U", OPROWCOL,
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));
FTensor::Index<'i', 3> i;
FTensor::Index<'j', 3> j;
FTensor::Index<'k', 3> k;
FTensor::Index<'l', 3> l;
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));
};
FTensor::Index<'i', 3> i;
FTensor::Index<'j', 3> j;
FTensor::Index<'k', 3> k;
FTensor::Index<'l', 3> l;
// 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
VectorDouble nF;
FTensor::Index<'i', 3> i;
// 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
// 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);
}
};
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,
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"
"-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(
new VolumeElementForcesAndSourcesCore(m_field));
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
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(
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 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()
constexpr double a
static PetscErrorCode ierr
constexpr int SPACE_DIM
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
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()
@ 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 ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ BLOCKSET
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int order
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1109
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:483
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:542
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.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 DMMoFEM.cpp:623
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:572
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition DMMoFEM.cpp:521
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
Definition DMMoFEM.cpp:664
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
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
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_NUMEREDDOF_ROW_FOR_LOOP_(PROBLEMPTR, IT)
use with loops to iterate row DOFs
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
double w(const double g, const double t)
constexpr AssemblyType A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
ApplyDirichletBc(const Range &fix_faces, const Range &fix_nodes, const Range &fix_second_node)
MoFEMErrorCode postProcess()
function is run at the end of loop
Set integration rule to boundary elements.
int operator()(int, int, int) 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.
@ OPROW
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:27
MoFEMErrorCode buildProblem()
Build problem.
Definition Simple.cpp:597
MoFEMErrorCode defineFiniteElements()
Define finite elements.
Definition Simple.cpp:380
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:194
MoFEMErrorCode buildFiniteElements()
Build finite elements.
Definition Simple.cpp:571
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:264
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:667
MoFEMErrorCode buildFields()
Build fields.
Definition Simple.cpp:482
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:473
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
Definition Simple.cpp:452
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:341
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.
OpPressure(const double pressure_val=1)
Set integration rule.
int operator()(int, int, int) const