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_NULL);
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 nme:nme847.
#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 addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
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.
Mat & snes_B
preconditioner of jacobian matrix
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.