|
| v0.14.0
|
Go to the documentation of this file.
22 using namespace MoFEM;
24 #include <boost/program_options.hpp>
26 namespace po = boost::program_options;
28 #include <boost/numeric/ublas/vector_proxy.hpp>
29 #include <boost/numeric/ublas/matrix.hpp>
30 #include <boost/numeric/ublas/matrix_proxy.hpp>
31 #include <boost/numeric/ublas/vector.hpp>
32 #include <boost/numeric/ublas/symmetric.hpp>
34 #include <adolc/adolc.h>
50 using namespace boost::numeric;
52 static char help[] =
"...\n"
55 int main(
int argc,
char *argv[]) {
63 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
65 pcomm =
new ParallelComm(&moab, PETSC_COMM_WORLD);
67 PetscBool flg = PETSC_TRUE;
71 if (flg != PETSC_TRUE) {
72 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
77 if (flg != PETSC_TRUE) {
80 PetscInt bubble_order;
83 if (flg != PETSC_TRUE) {
89 PetscBool is_partitioned = PETSC_FALSE;
91 &is_partitioned, &flg);
94 if (is_partitioned == PETSC_TRUE) {
96 option =
"PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION="
97 "PARALLEL_PARTITION;";
132 PetscPrintf(PETSC_COMM_WORLD,
"young_modulus = %4.2e \n",
young_modulus);
133 PetscPrintf(PETSC_COMM_WORLD,
"poisson_ratio = %4.2e \n",
poisson_ratio);
134 PetscPrintf(PETSC_COMM_WORLD,
"sIgma_y = %4.2e \n", cp.
sIgma_y);
135 PetscPrintf(PETSC_COMM_WORLD,
"H = %4.2e \n", cp.
H);
136 PetscPrintf(PETSC_COMM_WORLD,
"K = %4.2e \n", cp.
K);
137 PetscPrintf(PETSC_COMM_WORLD,
"phi = %4.2e \n", cp.
phi);
139 cp.
sTrain.resize(6,
false);
155 vector<BitRefLevel> bit_levels;
161 Range CubitSideSets_meshsets;
162 CHKERR m_field.get_cubit_meshsets(
SIDESET, CubitSideSets_meshsets);
191 m_field,
"MESH_NODE_POSITIONS");
192 CHKERR m_field.
loop_dofs(
"MESH_NODE_POSITIONS", ent_method_material, 0);
205 "PLASTIC",
"MESH_NODE_POSITIONS");
209 CHKERR MetaNeummanForces::addNeumannBCElements(m_field,
"DISPLACEMENT");
218 DMType dm_name =
"PLASTIC_PROB";
222 CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
223 CHKERR DMSetType(dm, dm_name);
227 CHKERR DMSetFromOptions(dm);
244 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
245 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
247 CHKERR MatZeroEntries(Aij);
257 PetscBool bbar = PETSC_TRUE;
264 "DISPLACEMENT", small_strain_plasticity.
commonData));
267 m_field,
"DISPLACEMENT", small_strain_plasticity.
commonData,
271 "DISPLACEMENT", small_strain_plasticity.
commonData));
274 "DISPLACEMENT", small_strain_plasticity.
commonData));
277 m_field,
"DISPLACEMENT", small_strain_plasticity.
commonData,
281 "DISPLACEMENT", small_strain_plasticity.
commonData));
284 "DISPLACEMENT", small_strain_plasticity.
commonData));
287 m_field,
"DISPLACEMENT", small_strain_plasticity.
commonData,
291 "DISPLACEMENT", small_strain_plasticity.
commonData));
296 boost::ptr_map<string, NeummanForcesSurface> neumann_forces;
297 boost::ptr_map<string, NodalForce> nodal_forces;
298 boost::ptr_map<string, EdgeForce> edge_forces;
302 CHKERR MetaNeummanForces::setMomentumFluxOperators(
303 m_field, neumann_forces, PETSC_NULL,
"DISPLACEMENT");
310 for (boost::ptr_map<string, NeummanForcesSurface>::iterator mit =
311 neumann_forces.begin();
312 mit != neumann_forces.end(); mit++) {
313 mit->second->methodsOp.push_back(
316 for (boost::ptr_map<string, NodalForce>::iterator mit =
317 nodal_forces.begin();
318 mit != nodal_forces.end(); mit++) {
319 mit->second->methodsOp.push_back(
322 for (boost::ptr_map<string, EdgeForce>::iterator mit =
324 mit != edge_forces.end(); mit++) {
325 mit->second->methodsOp.push_back(
336 boost::ptr_map<string, NeummanForcesSurface>::iterator fit;
337 fit = neumann_forces.begin();
338 for (; fit != neumann_forces.end(); fit++) {
340 dm, fit->first.c_str(), &fit->second->getLoopFe(), NULL, NULL);
342 fit = nodal_forces.begin();
343 for (; fit != nodal_forces.end(); fit++) {
345 dm, fit->first.c_str(), &fit->second->getLoopFe(), NULL, NULL);
347 fit = edge_forces.begin();
348 for (; fit != edge_forces.end(); fit++) {
350 dm, fit->first.c_str(), &fit->second->getLoopFe(), NULL, NULL);
354 &small_strain_plasticity.
feRhs,
355 PETSC_NULL, PETSC_NULL);
363 dm,
"PLASTIC", &small_strain_plasticity.
feLhs, NULL, NULL);
372 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
377 CHKERR SNESSetFromOptions(snes);
389 double final_time = 1, delta_time = 0.1;
392 double delta_time0 = delta_time;
407 rval = moab.get_entities_by_type(meshset, MBVERTEX, nodes,
true);
409 node_set.merge(nodes);
411 PetscPrintf(PETSC_COMM_WORLD,
"Nb. nodes in load path: %u\n",
415 boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_rows =
416 problem_ptr->getNumeredDofsRows();
417 Range::iterator nit = node_set.begin();
418 for (; nit != node_set.end(); nit++) {
419 NumeredDofEntityByEnt::iterator dit, hi_dit;
420 dit = numered_dofs_rows->get<
Ent_mi_tag>().lower_bound(*nit);
421 hi_dit = numered_dofs_rows->get<
Ent_mi_tag>().upper_bound(*nit);
422 for (; dit != hi_dit; dit++) {
423 load_path_dofs_view.insert(*dit);
430 SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
431 for (;
t < final_time; step++) {
434 PetscPrintf(PETSC_COMM_WORLD,
"Step %d Time %6.4g final time %3.2g\n",
435 step,
t, final_time);
439 dirichlet_bc.
ts_t =
t;
441 boost::ptr_map<string, NeummanForcesSurface>::iterator fit;
442 fit = neumann_forces.begin();
443 for (; fit != neumann_forces.end(); fit++) {
444 fit->second->getLoopFe().ts_t =
t;
446 fit = nodal_forces.begin();
447 for (; fit != nodal_forces.end(); fit++) {
448 fit->second->getLoopFe().ts_t =
t;
450 fit = edge_forces.begin();
451 for (; fit != edge_forces.end(); fit++) {
452 fit->second->getLoopFe().ts_t =
t;
464 if (step == 0 || reason < 0) {
465 CHKERR SNESSetLagJacobian(snes, -2);
467 CHKERR SNESSetLagJacobian(snes, -1);
469 CHKERR SNESSolve(snes, PETSC_NULL,
D);
472 CHKERR SNESGetIterationNumber(snes, &its);
478 CHKERR PetscPrintf(PETSC_COMM_WORLD,
479 "number of Newton iterations = %D\n", its);
481 CHKERR SNESGetConvergedReason(snes, &reason);
489 const double gamma = 0.5;
490 const double max_reudction = 1;
491 const double min_reduction = 1e-1;
493 reduction = pow((
double)its_d / (
double)(its + 1), gamma);
494 if (delta_time >= max_reudction * delta_time0 && reduction > 1) {
496 }
else if (delta_time <= min_reduction * delta_time0 &&
503 "reduction delta_time = %6.4e delta_time = %6.4e\n",
504 reduction, delta_time);
505 delta_time *= reduction;
506 if (reduction > 1 && delta_time < min_reduction * delta_time0) {
507 delta_time = min_reduction * delta_time0;
513 dm,
"PLASTIC", &small_strain_plasticity.
feUpdate);
518 NumeredDofEntity_multiIndex_uid_view_ordered::iterator it,
520 it = load_path_dofs_view.begin();
521 hi_it = load_path_dofs_view.end();
522 for (; it != hi_it; it++) {
523 PetscPrintf(PETSC_COMM_WORLD,
524 "load_path %s [ %d ] %6.4e %6.4e %6.4e\n",
525 (*it)->getName().c_str(),
527 (*it)->getDofCoeffIdx(), (*it)->getFieldData(),
t,
537 std::ostringstream stm;
538 stm <<
"out_" << step <<
".h5m";
541 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"out file %s\n",
#define MYPCOMM_INDEX
default communicator number PCOMM
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.
double young_modulus
Young modulus.
moab::Interface & postProcMesh
#define CHKERRQ_MOAB(a)
check error code of MoAB function
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
VectorDouble plasticStrain
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Projection of edge entities with one mid-node on hierarchical basis.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
PetscErrorCode snesCreate()
Set Dirichlet boundary conditions on displacements.
virtual PetscErrorCode createMatAVecR()
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.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Deprecated interface functions.
DeprecatedCoreInterface Interface
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
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
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
#define CHKERR
Inline error check.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
implementation of Data Operators for Forces and Sources
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
int main(int argc, char *argv[])
Force scale operator for reading two columns.
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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.
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< const_mem_fun< NumeredDofEntity::interface_type_DofEntity, UId, &NumeredDofEntity::getLocalUniqueId > > > > NumeredDofEntity_multiIndex_uid_view_ordered
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Small Strain plasticity material models.
double poisson_ratio
Poisson ratio.
J2 plasticity (Kinematic Isotropic (Linear) Hardening)
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
Post-process fields on refined mesh.
constexpr double t
plate stiffness
Small strain finite element implementation.
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Operators and data structures for small strain plasticity.
VectorDouble internalVariables
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
MoFEMErrorCode getForceScale(const double ts_t, double &scale)
#define CATCH_ERRORS
Catch errors.
PetscErrorCode createInternalVariablesTag()
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
boost::ptr_vector< MethodForForceScaling > methodsOp
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
#define _IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(MESHSET_MANAGER, NAME, IT)
Iterator that loops over Cubit BlockSet having a particular name.
const double D
diffusivity
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
Interface for nonlinear (SNES) solver.
keeps basic data about problem
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.
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.