46 static Tensor2<double, 3, 3>
Omega;
47 static Tensor4<double, 3, 3, 3, 3>
Is;
61 template <
typename T1,
typename T2>
64 std::vector<double> gaps;
67 auto t_normal = roller.getNormal(t_coords, t_disp);
68 double dist = roller.getGap(t_coords);
72 auto it = std::min_element(gaps.begin(), gaps.end());
73 auto nb = distance(gaps.begin(), it);
105 CHKERR PetscOptionsInsertString(
107 "-mat_mumps_icntl_14 800 -mat_mumps_icntl_24 1 -mat_mumps_icntl_13 1 ");
109 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"",
"none");
111 CHKERR PetscOptionsScalar(
"-gravity_y",
"Gravity load in Y direction",
"",
117 CHKERR PetscOptionsScalar(
"-poisson",
"Poisson ratio",
"",
poisson,
122 double rot_vec[3] = {0, 0, 1};
124 PetscBool is_rotating;
128 if (nb_ent < 3 && is_rotating) {
130 "provide an appropriate number of entries for omega vector (3)");
134 Omega(
i,
k) = levi_civita<double>(
i,
j,
k) * t1_omega(
j);
139 "poisson ratio value not supported");
141 CHKERR PetscOptionsScalar(
"-sigmaY",
"SigmaY - limit stress",
"",
sigmaY,
143 CHKERR PetscOptionsScalar(
"-H",
"Hardening modulus",
"",
H, &
H, PETSC_NULL);
144 CHKERR PetscOptionsScalar(
"-cn_pl",
145 "cn_pl parameter for active set constrains",
"",
148 CHKERR PetscOptionsScalar(
"-scale_constraint",
"constraint scale",
"",
151 CHKERR PetscOptionsScalar(
"-C1_k",
"1st Backstress",
"",
C1_k, &
C1_k,
155 CHKERR PetscOptionsScalar(
"-b_iso",
"Isotropic exponent",
"",
b_iso, &
b_iso,
158 CHKERR PetscOptionsScalar(
"-spring",
159 "Springs stiffness on the boundary (contact)",
"",
163 CHKERR PetscOptionsScalar(
"-cn_cont",
164 "cn_cont parameter for active set constrains",
"",
166 CHKERR PetscOptionsScalar(
167 "-rollers_stop",
"Time at which rollers will stop moving",
"",
170 constexpr
int MAX_NB_OF_ROLLERS = 10;
173 const char *body_list[] = {
"plane",
"sphere",
"cylinder",
"torus",
174 "superquadric",
"stl",
"nurbs"};
175 for (
int ii = 0; ii < MAX_NB_OF_ROLLERS; ++ii) {
182 PetscBool flg = PETSC_FALSE;
184 string param =
"-body" + to_string(ii);
192 param =
"-coords" + to_string(ii);
194 ¢er_coords(0), &nb_coord, &rflg);
195 param =
"-move" + to_string(ii);
200 switch (body_type_id) {
241 ierr = PetscOptionsEnd();
262 extern std::map<int, BlockParamData>
mat_blocks;
279 SmartPetscObj<Vec> m_diag_vec)
296 boost::ptr_vector<RigidBodyData> &roller_data,
297 boost::ptr_vector<DataFromMove> &bc_data,
bool scale_mat)
305 CHKERR VecSetOption(ts_F, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
309 case CTX_TSSETIJACOBIAN: {
310 snes_ctx = CTX_SNESSETJACOBIAN;
316 case CTX_TSSETIFUNCTION: {
317 snes_ctx = CTX_SNESSETFUNCTION;
330 for (
auto &bdata :
bcData) {
331 bdata.scaledValues(0) = bdata.values;
336 if (time <= (*cache).rollers_stop_time) {
337 for (
auto &roller : (*cache).rollerDataVec) {
338 roller.BodyDispScaled = roller.rollerDisp;
340 roller.BodyDispScaled);
358 const Problem *problem_ptr;
359 auto dm =
simple->getDM();
366 auto diag_m_vec = mps_ctx->
mDiagVec;
368 vector<string> fields;
371 for (
auto const &field : (*fields_ptr)) {
372 const std::string field_name = field->getName();
373 fields.push_back(field_name);
379 auto calculate_diagonal_scale = [&]() {
382 const double *f_array;
383 CHKERR VecGetArrayRead(ts_F, &f_array);
385 problem_ptr->getNumeredRowDofsPtr()->get<Unique_mi_tag>();
386 const int nb_local = problem_ptr->getNbLocalDofsRow();
388 size_t nb_fields = fields.size();
389 vector<double> lnorms(nb_fields, 0);
390 vector<double> norms(nb_fields, 0);
392 for (
size_t f = 0;
f != nb_fields; ++
f) {
394 const auto lo_uid = FieldEntity::getLoBitNumberUId(bit_number);
395 const auto hi_uid = FieldEntity::getHiBitNumberUId(bit_number);
396 for (
auto lo = rows_index.lower_bound(lo_uid);
397 lo != rows_index.upper_bound(hi_uid); ++lo) {
398 const int local_idx = (*lo)->getPetscLocalDofIdx();
399 if (local_idx >= 0 && local_idx < nb_local) {
400 const double val = f_array[local_idx];
401 lnorms[
f] += val * val;
406 CHKERR VecRestoreArrayRead(ts_F, &f_array);
407 CHKERR MPIU_Allreduce(lnorms.data(), norms.data(), nb_fields, MPIU_REAL,
408 MPIU_SUM, PetscObjectComm((PetscObject)dm));
410 for (
auto &v : norms)
418 CHKERR VecSet(diag_m_vec, 1.);
420 double *diag_m_array;
421 CHKERR VecGetArray(diag_m_vec, &diag_m_array);
423 for (
size_t f = 0;
f != nb_fields; ++
f) {
425 const auto lo_uid = FieldEntity::getLoBitNumberUId(bit_number);
426 const auto hi_uid = FieldEntity::getHiBitNumberUId(bit_number);
427 for (
auto lo = rows_index.lower_bound(lo_uid);
428 lo != rows_index.upper_bound(hi_uid); ++lo) {
429 const int local_idx = (*lo)->getPetscLocalDofIdx();
430 if (local_idx >= 0 && local_idx < nb_local)
431 diag_m_array[local_idx] = 1. / norms[
f];
435 CHKERR VecRestoreArray(diag_m_vec, &diag_m_array);
438 case CTX_TSSETIJACOBIAN: {
439 CHKERR MatAssemblyBegin(ts_B, MAT_FINAL_ASSEMBLY);
440 CHKERR MatAssemblyEnd(ts_B, MAT_FINAL_ASSEMBLY);
441 CHKERR MatDiagonalScale(ts_B, diag_m_vec, PETSC_NULL);
444 case CTX_TSSETIFUNCTION: {
445 CHKERR VecPointwiseMult(ts_F, ts_F, diag_m_vec);
454 CHKERR calculate_diagonal_scale();
490 Range &essential_boundary_entsX, Range &essential_boundary_entsY,
491 Range &essential_boundary_entsZ)
492 :
mField(m_field),
entX(essential_boundary_entsX),
493 entY(essential_boundary_entsY),
entZ(essential_boundary_entsZ),
504 auto pre_proc_lhs = [&](Range &boundary_ents,
int comp) {
514 CHKERR MatZeroRowsIS(ts_B, bc_is, 0.0, PETSC_NULL, PETSC_NULL);
520 auto pre_proc_rhs = [&](Range &boundary_ents,
int comp) {
531 CHKERR VecISSet(ts_F, bc_is, 0.0);
540 case CTX_TSSETIJACOBIAN: {
541 CHKERR MatSetOption(ts_B, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
542 CHKERR MatSetOption(ts_B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
544 CHKERR MatAssemblyBegin(ts_B, MAT_FINAL_ASSEMBLY);
545 CHKERR MatAssemblyEnd(ts_B, MAT_FINAL_ASSEMBLY);
552 case CTX_TSSETIFUNCTION: {
553 CHKERR VecSetOption(ts_F, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
554 CHKERR VecAssemblyBegin(ts_F);
555 CHKERR VecAssemblyEnd(ts_F);
560 CHKERR VecAssemblyBegin(ts_F);
561 CHKERR VecAssemblyEnd(ts_F);
563 CHKERR VecGhostUpdateBegin(ts_F, ADD_VALUES, SCATTER_REVERSE);
564 CHKERR VecGhostUpdateEnd(ts_F, ADD_VALUES, SCATTER_REVERSE);
multi_index_container< boost::shared_ptr< Field >, indexed_by< hashed_unique< tag< BitFieldId_mi_tag >, const_mem_fun< Field, const BitFieldId &, &Field::getId >, HashBit< BitFieldId >, EqBit< BitFieldId > >, ordered_unique< tag< Meshset_mi_tag >, member< Field, EntityHandle, &Field::meshSet > >, ordered_unique< tag< FieldName_mi_tag >, const_mem_fun< Field, boost::string_ref, &Field::getNameRef > >, ordered_non_unique< tag< BitFieldId_space_mi_tag >, const_mem_fun< Field, FieldSpace, &Field::getSpace > > > > Field_multiIndex
Field_multiIndex for Field.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#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 ...
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::vector< double, DoubleAllocator > VectorDouble
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 SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
MoFEMErrorCode scaleParameters()
MoFEMErrorCode getOptionsFromCommandLine()
static Tensor4< double, 3, 3, 3, 3 > Is
virtual MoFEMErrorCode scaleParameters()=0
Ddg< double, 3, 3 > scaling_mat
static double spring_stiffness
static int closest_roller_id
static Tensor2< double, 3, 3 > Omega
static RigidBodyData * closest_roller
int getClosestRollerID(Tensor1< T1, 3 > &t_coords, Tensor1< T2, 3 > &t_disp)
static boost::ptr_vector< RigidBodyData > rollerDataVec
DataFromMove(size_t c, double val, Range ents_bc)
VectorDouble scaledValues
EraseRows(MoFEM::Interface &m_field, std::string problem_name, Range &essential_boundary_entsX, Range &essential_boundary_entsY, Range &essential_boundary_entsZ)
MoFEMErrorCode operator()()
MoFEMErrorCode preProcess()
MoFEMErrorCode postProcess()
MoFEM::Interface & mField
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
MoFEM::Interface & mField
MoFEMErrorCode preProcess()
MoFEMErrorCode postProcess()
FePrePostProcess(MoFEM::Interface &m_field, boost::ptr_vector< RigidBodyData > &roller_data, boost::ptr_vector< DataFromMove > &bc_data, bool scale_mat)
boost::ptr_vector< RigidBodyData > & rollerDataVec
boost::ptr_vector< DataFromMove > & bcData
boost::ptr_vector< MethodForForceScaling > methodsOp
LoadScale(double &my_scale)
MoFEMErrorCode scaleNf(const FEMethod *fe, VectorDouble &nf)
Class used to scale loads, f.e. in arc-length control.
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
Deprecated interface functions.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MpSnesCtx(MoFEM::Interface &m_field, const std::string &problem_name, SmartPetscObj< Vec > m_diag_vec)
SmartPetscObj< Vec > mDiagVec