18static char help[] = 
"...\n\n";
 
   20template <
typename TYPE>
 
   37      boost::shared_ptr<const NumeredEntFiniteElement> fe_ptr) {
 
   44      for (
int rr = 0; rr < 3; rr++) {
 
   45        for (
int cc = 0; cc < 3; cc++) {
 
   50    if (
D_mu.size1() == 0) {
 
   53      for (
int rr = 0; rr < 6; rr++) {
 
   54        D_mu(rr, rr) = rr < 3 ? 2 : 1;
 
   65      sTrain[1] = this->
F(1, 1) - 1;
 
   66      sTrain[2] = this->
F(2, 2) - 1;
 
   67      sTrain[3] = this->
F(0, 1) + this->
F(1, 0);
 
   68      sTrain[4] = this->
F(1, 2) + this->
F(2, 1);
 
   69      sTrain[5] = this->
F(0, 2) + this->
F(2, 0);
 
   72      for (
int ii = 0; ii != 6; ++ii)
 
   73        for (
int jj = 0; jj != 6; ++jj)
 
   76      this->
P(0, 0) = sTress[0];
 
   77      this->
P(1, 1) = sTress[1];
 
   78      this->
P(2, 2) = sTress[2];
 
   79      this->
P(0, 1) = this->
P(1, 0) = sTress[3];
 
   80      this->
P(1, 2) = this->
P(2, 1) = sTress[4];
 
   81      this->
P(0, 2) = this->
P(2, 0) = sTress[5];
 
   88      for (
int ii = 0; ii != 6; ++ii)
 
   89        for (
int jj = 0; jj != 6; ++jj)
 
  104      this->
t_P(i, 
j) *= 
J;
 
 
 
  125      this->
sTrain0[3] <<= (G0(1, 0) + G0(0, 1));
 
  126      this->
sTrain0[4] <<= (G0(2, 1) + G0(1, 2));
 
  127      this->
sTrain0[5] <<= (G0(2, 0) + G0(0, 2));
 
  129      nb_active_variables += 6;
 
  131    } 
catch (
const std::exception &ex) {
 
  132      std::ostringstream ss;
 
  133      ss << 
"throw in method: " << ex.what() << std::endl;
 
  134      SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
 
 
  148      active_variables[shift + 0] = G0(0, 0);
 
  149      active_variables[shift + 1] = G0(1, 1);
 
  150      active_variables[shift + 2] = G0(2, 2);
 
  151      active_variables[shift + 3] = G0(0, 1) + G0(1, 0);
 
  152      active_variables[shift + 4] = G0(1, 2) + G0(2, 1);
 
  153      active_variables[shift + 5] = G0(0, 2) + G0(2, 0);
 
  155    } 
catch (
const std::exception &ex) {
 
  156      std::ostringstream ss;
 
  157      ss << 
"throw in method: " << ex.what() << std::endl;
 
  158      SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
 
 
 
  165int main(
int argc, 
char *argv[]) {
 
  168  const string default_options = 
"-ksp_type fgmres \n" 
  170                                 "-pc_factor_mat_solver_type mumps \n" 
  171                                 "-mat_mumps_icntl_20 0 \n" 
  175                                 "-snes_type newtonls \n" 
  176                                 "-snes_linesearch_type basic \n" 
  177                                 "-snes_max_it 100 \n" 
  183  string param_file = 
"param_file.petsc";
 
  184  if (!
static_cast<bool>(ifstream(param_file))) {
 
  185    std::ofstream file(param_file.c_str(), std::ios::ate);
 
  186    if (file.is_open()) {
 
  187      file << default_options;
 
  192  SlepcInitialize(&argc, &argv, param_file.c_str(), 
help);
 
  197    moab::Core mb_instance;
 
  198    moab::Interface &moab = mb_instance;
 
  199    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, 
MYPCOMM_INDEX);
 
  200    auto moab_comm_wrap =
 
  201        boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, 
false);
 
  203      pcomm = 
new ParallelComm(&moab, moab_comm_wrap->get_comm());
 
  205    PetscBool flg = PETSC_TRUE;
 
  209    if (flg != PETSC_TRUE) {
 
  210      SETERRQ(PETSC_COMM_SELF, 1, 
"*** ERROR -my_file (MESH FILE NEEDED)");
 
  215    PetscBool is_partitioned = PETSC_FALSE;
 
  217                               &is_partitioned, &flg);
 
  219    if (is_partitioned == PETSC_TRUE) {
 
  222      option = 
"PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=" 
  223               "PARALLEL_PARTITION;";
 
  225      CHKERR pcomm->resolve_shared_ents(0, 3, 0);
 
  226      CHKERR pcomm->resolve_shared_ents(0, 3, 1);
 
  227      CHKERR pcomm->resolve_shared_ents(0, 3, 2);
 
  237    Range CubitSIDESETs_meshsets;
 
  239        SIDESET, CubitSIDESETs_meshsets);
 
  245    CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
 
  260    bool check_if_spatial_field_exist = m_field.
check_field(
"SPATIAL_POSITION");
 
  273    boost::shared_ptr<Hooke<double>> mat_double =
 
  274        boost::make_shared<Hooke<double>>();
 
  275    boost::shared_ptr<MyMat<adouble>> mat_adouble =
 
  276        boost::make_shared<MyMat<adouble>>();
 
  279    CHKERR elastic.setBlocks(mat_double, mat_adouble);
 
  280    CHKERR elastic.addElement(
"ELASTIC", 
"SPATIAL_POSITION");
 
  285    elastic.feRhs.getOpPtrVector().push_back(
 
  287            "D0", elastic.commonData));
 
  288    elastic.feLhs.getOpPtrVector().push_back(
 
  290            "D0", elastic.commonData));
 
  291    CHKERR elastic.setOperators(
"SPATIAL_POSITION");
 
  307    if (flg != PETSC_TRUE) {
 
  332                                                        "MESH_NODE_POSITIONS");
 
  338      CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris, 
true);
 
  345      CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris, 
true);
 
  357    if (!check_if_spatial_field_exist) {
 
  359                                                        "MESH_NODE_POSITIONS");
 
  377    if (is_partitioned) {
 
  392        "ELASTIC_MECHANICS", 
ROW, &
F);
 
  397        ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>(
"ELASTIC_MECHANICS",
 
  416    CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
 
  417    CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
 
  418    CHKERR MatZeroEntries(Aij);
 
  421        "ELASTIC_MECHANICS", 
COL, 
D, INSERT_VALUES, SCATTER_FORWARD);
 
  422    CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
 
  423    CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
 
  432    CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
 
  433    CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
 
  435        "ELASTIC_MECHANICS", 
COL, 
D, INSERT_VALUES, SCATTER_REVERSE);
 
  438    boost::ptr_map<std::string, NodalForce> nodal_forces;
 
  439    string fe_name_str = 
"FORCE_FE";
 
  440    nodal_forces.insert(fe_name_str, 
new NodalForce(m_field));
 
  443    boost::ptr_map<std::string, NodalForce>::iterator fit =
 
  444        nodal_forces.begin();
 
  445    for (; fit != nodal_forces.end(); fit++) {
 
  447                                          fit->second->getLoopFe());
 
  456    elastic.getLoopFeRhs().snes_x = 
D;
 
  457    elastic.getLoopFeRhs().snes_f = 
F;
 
  459                                        elastic.getLoopFeRhs());
 
  467    my_Dirichlet_bc.
snes_B = Aij;
 
  477    elastic.getLoopFeLhs().snes_B = Aij;
 
  479                                        elastic.getLoopFeLhs());
 
  486    CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
 
  487    CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
 
  489    CHKERR MatAssemblyBegin(Aij, MAT_FINAL_ASSEMBLY);
 
  490    CHKERR MatAssemblyEnd(Aij, MAT_FINAL_ASSEMBLY);
 
  500    CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
 
  501    CHKERR KSPSetOperators(solver, Aij, Aij);
 
  502    CHKERR KSPSetFromOptions(solver);
 
  507    CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
 
  508    CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
 
  511    CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
 
  512    CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
 
  515        "ELASTIC_MECHANICS", 
"SPATIAL_POSITION", 
"D0", 
COL, 
D, INSERT_VALUES,
 
  519    CHKERR MatDuplicate(Aij, MAT_SHARE_NONZERO_PATTERN, &Bij);
 
  521    CHKERR MatZeroEntries(Bij);
 
  540    mat_adouble->doAotherwiseB = 
false;
 
  543    my_Dirichlet_bc.
snes_B = Bij;
 
  549    PetscBool is_conservative = PETSC_FALSE;
 
  551                               &is_conservative, &flg);
 
  552    if (is_conservative) {
 
  558    elastic.getLoopFeLhs().snes_B = Bij;
 
  560                                        elastic.getLoopFeLhs());
 
  565    CHKERR MatSetOption(Bij, MAT_SPD, PETSC_TRUE);
 
  566    CHKERR MatAssemblyBegin(Bij, MAT_FINAL_ASSEMBLY);
 
  567    CHKERR MatAssemblyEnd(Bij, MAT_FINAL_ASSEMBLY);
 
  579    PetscInt nev, maxit, its;
 
  584    CHKERR EPSCreate(PETSC_COMM_WORLD, &
eps);
 
  602    PetscPrintf(PETSC_COMM_WORLD, 
" Number of iterations of the method: %D\n",
 
  609    PetscPrintf(PETSC_COMM_WORLD, 
" Solution method: %s\n", type);
 
  610    CHKERR EPSGetDimensions(
eps, &nev, NULL, NULL);
 
  611    PetscPrintf(PETSC_COMM_WORLD, 
" Number of requested eigenvalues: %D\n",
 
  614    PetscPrintf(PETSC_COMM_WORLD, 
" Stopping condition: tol=%.4g, maxit=%D\n",
 
  631    PetscScalar eigr, eigi, nrm2r;
 
  632    for (
int nn = 0; nn < nev; nn++) {
 
  633      CHKERR EPSGetEigenpair(
eps, nn, &eigr, &eigi, 
D, PETSC_NULLPTR);
 
  634      CHKERR VecNorm(
D, NORM_2, &nrm2r);
 
  637          " ncov = %D eigr = %.4g eigi = %.4g (inv eigr = %.4g) nrm2r = %.4g\n",
 
  638          nn, eigr, eigi, 1. / eigr, nrm2r);
 
  639      std::ostringstream o1;
 
  640      o1 << 
"eig_" << nn << 
".h5m";
 
  642          "ELASTIC_MECHANICS", 
"SPATIAL_POSITION", 
"EIGEN_VECTOR", 
COL, 
D,
 
  643          INSERT_VALUES, SCATTER_REVERSE);
 
  649    CHKERR KSPDestroy(&solver);
 
 
Implementation of linear elastic material.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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 add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string name, const bool recursive=true)=0
add entities to finite element
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 modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual bool check_field(const std::string &name) const =0
check if field is in database
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULLPTR)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULLPTR)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode writeFile(const std::string file_name, const char *file_type="MOAB", const char *file_options="PARALLEL=WRITE_PART")
wrote results in (MOAB) format, use "file_name.h5m"
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
MoFEMErrorCode buildProblemOnDistributedMesh(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures, assuming that mesh is distributed (collective)
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::vector< T, ublas::bounded_array< T, N > > VectorBoundedArray
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Set Dirichlet boundary conditions on spatial displacements.
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Matrix manager is used to build and partition problems.
Interface for managing meshsets containing materials and boundary conditions.
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
@ CTX_SNESSETJACOBIAN
Setting up Jacobian matrix computation.
@ CTX_SNESSETFUNCTION
Setting up nonlinear function evaluation.
Vec & snes_f
Reference to residual vector.
Vec & snes_x
Reference to current solution state vector.
Mat & snes_B
Reference to preconditioner of Jacobian matrix.
SNESContext snes_ctx
Current SNES computation context.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
VectorBoundedArray< TYPE, 6 > sTrain
FTensor::Index< 'k', 3 > k
virtual MoFEMErrorCode calculateP_PiolaKirchhoffI(const NonlinearElasticElement::BlockData block_data, boost::shared_ptr< const NumeredEntFiniteElement > fe_ptr)
VectorBoundedArray< TYPE, 6 > sTrain0
FTensor::Index< 'i', 3 > i
FTensor::Index< 'j', 3 > j
VectorBoundedArray< TYPE, 6 > sTress
virtual MoFEMErrorCode setUserActiveVariables(VectorDouble &active_variables)
virtual MoFEMErrorCode setUserActiveVariables(int &nb_active_variables)
add additional active variables
MoFEMErrorCode addPressure(int ms_id)
MoFEMErrorCode addForce(int ms_id)
NonLinear surface pressure element (obsolete implementation)
MyTriangleSpatialFE & getLoopSpatialFe()
data for calculation heat conductivity and heat capacity elements
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
Implementation of elastic (non-linear) St. Kirchhoff equation.
int gG
Gauss point number.
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_sigmaCauchy
MatrixBoundedArray< TYPE, 9 > sigmaCauchy
CommonData * commonDataPtr
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_P
FTensor::Tensor2< FTensor::PackPtr< TYPE *, 0 >, 3, 3 > t_invF
MatrixBoundedArray< TYPE, 9 > invF
MatrixBoundedArray< TYPE, 9 > F
MatrixBoundedArray< TYPE, 9 > P
structure grouping operators and data used for calculation of nonlinear elastic element
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.