v0.10.0
BasicFeTools.hpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2  * MoFEM is free software: you can redistribute it and/or modify it under
3  * the terms of the GNU Lesser General Public License as published by the
4  * Free Software Foundation, either version 3 of the License, or (at your
5  * option) any later version.
6  *
7  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10  * License for more details.
11  *
12  * You should have received a copy of the GNU Lesser General Public
13  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14 // #include <BasicFiniteElements.hpp>
15 // namespace ContactPlasticity {
16 #pragma once
17 
18 struct CacheParams {
19 
20  double young_modulus;
22  double poisson;
23  double density;
24  // for plastic
25  double sigmaY;
26  double cn_pl;
27  // hardening
28  double H;
29  double C1_k;
30  double Q_inf;
31  double b_iso;
33 
34  Ddg<double, 3, 3> tD;
35  Ddg<double, 3, 3> scaling_mat; // for flow rule
36 
37  // for contact
38  double cn_cont;
40 
41  // for rotation
42  static double gravity_y;
43  static double spring_stiffness;
44 
45  static int closest_roller_id;
46  static Tensor2<double, 3, 3> Omega;
47  static Tensor4<double, 3, 3, 3, 3> Is; // symmetric
48 
49  static boost::ptr_vector<RigidBodyData> rollerDataVec;
51 
54  sigmaY(1e9), cn_pl(1), H(0), C1_k(0), Q_inf(0), b_iso(0),
55  cn_cont(1. / young_modulus), /*spring_stiffness(0),
56  closest_roller_id(0), gravity_y(0),*/
57  rollers_stop_time(INFINITY), scale_constraint(1) {}
58 
60 
61  template <typename T1, typename T2>
63 
64  std::vector<double> gaps;
66  for (auto &roller : rollerDataVec) {
67  auto t_normal = roller.getNormal(t_coords, t_disp);
68  double dist = roller.getGap(t_coords);
69  gaps.push_back(dist);
70  }
71 
72  auto it = std::min_element(gaps.begin(), gaps.end());
73  auto nb = distance(gaps.begin(), it);
74 
75  return nb;
76  };
77 };
78 
87 
94 
95 static Number<0> N0;
96 static Number<1> N1;
97 static Number<2> N2;
98 
99 struct BlockParamData : public CacheParams {
100 
103 
104  // CHKERR PetscOptionsInsertString(NULL, default_options);
105  CHKERR PetscOptionsInsertString(
106  NULL,
107  "-mat_mumps_icntl_14 800 -mat_mumps_icntl_24 1 -mat_mumps_icntl_13 1 ");
108 
109  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "none");
110 
111  CHKERR PetscOptionsScalar("-gravity_y", "Gravity load in Y direction", "",
112  gravity_y, &gravity_y, PETSC_NULL);
113 
114  // for elastic
115  CHKERR PetscOptionsScalar("-young", "Young's modulus", "", young_modulus,
116  &young_modulus, PETSC_NULL);
117  CHKERR PetscOptionsScalar("-poisson", "Poisson ratio", "", poisson,
118  &poisson, PETSC_NULL);
119  CHKERR PetscOptionsScalar("-density", "Density", "", density, &density,
120  PETSC_NULL);
121 
122  double rot_vec[3] = {0, 0, 1};
123  int nb_ent = 3;
124  PetscBool is_rotating;
125  // get rotation velocity vector
126  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-rot_omega", rot_vec, &nb_ent,
127  &is_rotating);
128  if (nb_ent < 3 && is_rotating) {
129  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
130  "provide an appropriate number of entries for omega vector (3)");
131  }
132  {
133  Tensor1<double, 3> t1_omega(rot_vec[0], rot_vec[1], rot_vec[2]);
134  Omega(i, k) = levi_civita<double>(i, j, k) * t1_omega(j);
135  }
136 
137  if (poisson >= 0.5)
138  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
139  "poisson ratio value not supported");
140 
141  CHKERR PetscOptionsScalar("-sigmaY", "SigmaY - limit stress", "", sigmaY,
142  &sigmaY, PETSC_NULL);
143  CHKERR PetscOptionsScalar("-H", "Hardening modulus", "", H, &H, PETSC_NULL);
144  CHKERR PetscOptionsScalar("-cn_pl",
145  "cn_pl parameter for active set constrains", "",
146  cn_pl, &cn_pl, PETSC_NULL);
147 
148  CHKERR PetscOptionsScalar("-scale_constraint", "constraint scale", "",
149  scale_constraint, &scale_constraint, PETSC_NULL);
150 
151  CHKERR PetscOptionsScalar("-C1_k", "1st Backstress", "", C1_k, &C1_k,
152  PETSC_NULL);
153  CHKERR PetscOptionsScalar("-Q_inf", "Saturation stress", "", Q_inf, &Q_inf,
154  PETSC_NULL);
155  CHKERR PetscOptionsScalar("-b_iso", "Isotropic exponent", "", b_iso, &b_iso,
156  PETSC_NULL);
157  // for contact
158  CHKERR PetscOptionsScalar("-spring",
159  "Springs stiffness on the boundary (contact)", "",
160  spring_stiffness, &spring_stiffness, PETSC_NULL);
161  cn_cont = 1. / young_modulus;
162 
163  CHKERR PetscOptionsScalar("-cn_cont",
164  "cn_cont parameter for active set constrains", "",
165  cn_cont, &cn_cont, PETSC_NULL);
166  CHKERR PetscOptionsScalar(
167  "-rollers_stop", "Time at which rollers will stop moving", "",
168  rollers_stop_time, &rollers_stop_time, PETSC_NULL);
169 
170  constexpr int MAX_NB_OF_ROLLERS = 10;
171  // if (with_contact) {
172 
173  const char *body_list[] = {"plane", "sphere", "cylinder", "torus",
174  "superquadric", "stl", "nurbs"};
175  for (int ii = 0; ii < MAX_NB_OF_ROLLERS; ++ii) {
176  int nb_coord = 3;
177  int nb_disp = 3;
178  int nb_dirs = 3;
179  VectorDouble center_coords({0, 0, 0});
180  VectorDouble disp({0, 0, 0});
181  VectorDouble direction({0, 0, 1});
182  PetscBool flg = PETSC_FALSE;
183  PetscBool rflg;
184  string param = "-body" + to_string(ii);
185  // char mesh_file_name[255];
186  // CHKERR PetscOptionsString(param.c_str(), "file name", "", "mesh.vtk",
187  // mesh_file_name, 255, &flg);
188  int body_type_id = LASTBODY;
189  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, param.c_str(), body_list,
190  LASTBODY, &body_type_id, &flg);
191 
192  param = "-coords" + to_string(ii);
193  CHKERR PetscOptionsGetRealArray(NULL, NULL, param.c_str(),
194  &center_coords(0), &nb_coord, &rflg);
195  param = "-move" + to_string(ii);
196  CHKERR PetscOptionsGetRealArray(NULL, NULL, param.c_str(), &disp(0),
197  &nb_disp, &rflg);
198  if (flg) {
199 
200  switch (body_type_id) {
201  case PLANE:
202  rollerDataVec.push_back(new PlaneRigidBody(center_coords, disp, ii));
203  break;
204  case SPHERE:
205  rollerDataVec.push_back(new SphereRigidBody(center_coords, disp, ii));
206  break;
207  case CYLINDER:
208  rollerDataVec.push_back(
209  new CylinderRigidBody(center_coords, disp, ii));
210  break;
211  case TORUS:
212  rollerDataVec.push_back(new TorusRigidBody(center_coords, disp, ii));
213  break;
214  case SUPERQUADRIC:
215  rollerDataVec.push_back(
216  new SuperquadricRigidBody(center_coords, disp, ii));
217  break;
218  case STL:
219  rollerDataVec.push_back(new STLRigidBody(center_coords, disp, ii));
220  break;
221  case NURBS:
222  // code to be executed if
223  break;
224  default:
225  // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
226  // "this type of rigid body is not supported, check spelling "
227  // "(plane, sphere, cylinder, torus, superquadric, stl, nurbs)");
228  break;
229  }
230  CHKERR rollerDataVec.back().getBodyOptions();
231  }
232  }
233 
234  // if (rollerDataVec.empty()) {
235  // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
236  // "provide at least one roller data for contact (-roller "
237  // "x,y -radius r and opional -move x,y) ");
238  // }
239  // }
240 
241  ierr = PetscOptionsEnd();
242  CHKERRQ(ierr);
243 
245  }
246 
249 
251  young_modulus = 1;
253  H *= young_modulus_inv;
256 
258  }
259 };
260 
261 extern BlockParamData *cache;
262 extern std::map<int, BlockParamData> mat_blocks; // for use in future
263 
264 struct DataFromMove {
265  size_t comp;
266  double values;
268  Range ents;
269  DataFromMove(size_t c, double val, Range ents_bc)
270  : comp(c), values(val), scaledValues(1), ents(ents_bc) {
271  scaledValues(0) = val;
272  }
273  DataFromMove() = delete;
274 };
275 
276 struct MpSnesCtx : public SnesCtx {
277  SmartPetscObj<Vec> mDiagVec;
278  MpSnesCtx(MoFEM::Interface &m_field, const std::string &problem_name,
279  SmartPetscObj<Vec> m_diag_vec)
280  : SnesCtx(m_field, problem_name), mDiagVec(m_diag_vec) {}
281 };
282 
283 /**
284  * \brief Not used at this stage. Could be used to do some calculations,
285  * before assembly of local elements.
286  */
287 struct FePrePostProcess : public FEMethod {
288 
290  boost::ptr_vector<RigidBodyData> &rollerDataVec;
291  boost::ptr_vector<DataFromMove> &bcData;
292  boost::ptr_vector<MethodForForceScaling> methodsOp;
293  bool scaleMat;
296  boost::ptr_vector<RigidBodyData> &roller_data,
297  boost::ptr_vector<DataFromMove> &bc_data, bool scale_mat)
298  : mField(m_field), rollerDataVec(roller_data), bcData(bc_data),
299  scaleMat(scale_mat) {}
300 
302  //
304  runOnceFlag = true;
305  CHKERR VecSetOption(ts_F, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
306 
307 
308  switch (ts_ctx) {
309  case CTX_TSSETIJACOBIAN: {
310  snes_ctx = CTX_SNESSETJACOBIAN;
311  snes_B = ts_B;
312  // CHKERR MatSetOption(ts_B, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
313  // CHKERR MatSetOption(ts_B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
314  break;
315  }
316  case CTX_TSSETIFUNCTION: {
317  snes_ctx = CTX_SNESSETFUNCTION;
318  snes_f = ts_F;
319  break;
320  }
321  default:
322  break;
323  }
324 
325  auto time = ts_t;
326 
327  // forceDataScaled = forceData;
328  // pressureDataScaled = pressureData;
329 
330  for (auto &bdata : bcData) {
331  bdata.scaledValues(0) = bdata.values;
333  bdata.scaledValues);
334  }
335 
336  if (time <= (*cache).rollers_stop_time) {
337  for (auto &roller : (*cache).rollerDataVec) {
338  roller.BodyDispScaled = roller.rollerDisp;
340  roller.BodyDispScaled);
341  }
342  }
344  }
345 
347 
349  if (runOnceFlag) {
350  runOnceFlag = false;
351  } else
353 
354  if (!scaleMat)
356 
357  Simple *simple = mField.getInterface<Simple>();
358  const Problem *problem_ptr;
359  auto dm = simple->getDM();
360  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
361 
362  SnesCtx *snes_ctx;
363  CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
364  MpSnesCtx *mps_ctx = (MpSnesCtx *)snes_ctx;
365 
366  auto diag_m_vec = mps_ctx->mDiagVec;
367 
368  vector<string> fields;
369  const Field_multiIndex *fields_ptr;
370  CHKERR mField.get_fields(&fields_ptr);
371  for (auto const &field : (*fields_ptr)) {
372  const std::string field_name = field->getName();
373  fields.push_back(field_name);
374  }
375 
376  // if (fields.size() < 2)
377  // MoFEMFunctionReturnHot(0);
378 
379  auto calculate_diagonal_scale = [&]() {
381 
382  const double *f_array;
383  CHKERR VecGetArrayRead(ts_F, &f_array);
384  auto &rows_index =
385  problem_ptr->getNumeredRowDofsPtr()->get<Unique_mi_tag>();
386  const int nb_local = problem_ptr->getNbLocalDofsRow();
387 
388  size_t nb_fields = fields.size();
389  vector<double> lnorms(nb_fields, 0);
390  vector<double> norms(nb_fields, 0);
391 
392  for (size_t f = 0; f != nb_fields; ++f) {
393  const auto bit_number = mField.get_field_bit_number(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;
402  }
403  }
404  }
405 
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));
409 
410  for (auto &v : norms)
411  v = sqrt(v);
412 
413  // for (size_t f = 0; f != nb_fields; ++f)
414  // MOFEM_LOG_C("WORLD", Sev::verbose,
415  // "Scaling diagonal for field[ %s ] by %9.8e",
416  // fields[f].c_str(), norms[f]);
417 
418  CHKERR VecSet(diag_m_vec, 1.);
419 
420  double *diag_m_array;
421  CHKERR VecGetArray(diag_m_vec, &diag_m_array);
422 
423  for (size_t f = 0; f != nb_fields; ++f) {
424  const auto bit_number = mField.get_field_bit_number(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];
432  }
433  }
434 
435  CHKERR VecRestoreArray(diag_m_vec, &diag_m_array);
436 
437  switch (ts_ctx) {
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);
442  break;
443  }
444  case CTX_TSSETIFUNCTION: {
445  CHKERR VecPointwiseMult(ts_F, ts_F, diag_m_vec);
446  }
447  default:
448  break;
449  }
450 
452  };
453 
454  CHKERR calculate_diagonal_scale();
455 
457  }
458 };
459 
460 static MoFEMErrorCode custom_snes_rhs(SNES snes, Vec x, Vec f, void *ctx) {
461  MpSnesCtx *snes_ctx = (MpSnesCtx *)ctx;
463  CHKERR SnesRhs(snes, x, f, ctx);
464  CHKERR VecPointwiseMult(f, f, snes_ctx->mDiagVec);
466 }
467 
468 static MoFEMErrorCode custom_snes_mat(SNES snes, Vec x, Mat A, Mat B,
469  void *ctx) {
470  MpSnesCtx *snes_ctx = (MpSnesCtx *)ctx;
472  CHKERR SnesMat(snes, x, A, B, ctx);
473  CHKERR MatDiagonalScale(B, snes_ctx->mDiagVec, PETSC_NULL);
474 
476 };
477 
479  double &scale;
480  LoadScale(double &my_scale) : scale(my_scale){};
483  nf *= scale;
485  }
486 };
487 
488 struct EraseRows : public FEMethod {
489  EraseRows(MoFEM::Interface &m_field, std::string problem_name,
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),
494  problemName(problem_name) {}
495 
497  eraseFlag = true;
498  return 0;
499  }
500 
503 
504  auto pre_proc_lhs = [&](Range &boundary_ents, int comp) {
506 
507  IS bc_is;
508 
509  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
510  problemName, ROW, "U", comp, comp, &bc_is, &boundary_ents);
511 
512  // CHKERR ISView(bc_is, PETSC_VIEWER_STDOUT_SELF);
513 
514  CHKERR MatZeroRowsIS(ts_B, bc_is, 0.0, PETSC_NULL, PETSC_NULL);
515 
516  CHKERR ISDestroy(&bc_is);
518  };
519 
520  auto pre_proc_rhs = [&](Range &boundary_ents, int comp) {
522 
523  IS bc_is;
524 
525  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
526  problemName, ROW, "U", comp, comp, &bc_is, &boundary_ents);
527 
528  // CHKERR ISView(bc_is, PETSC_VIEWER_STDOUT_SELF);
529  // CHKERR VecView(ts_F,PETSC_VIEWER_STDOUT_WORLD);
530 
531  CHKERR VecISSet(ts_F, bc_is, 0.0);
532 
533  CHKERR ISDestroy(&bc_is);
535  };
536 
537  if (eraseFlag) {
538 
539  switch (ts_ctx) {
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);
543 
544  CHKERR MatAssemblyBegin(ts_B, MAT_FINAL_ASSEMBLY);
545  CHKERR MatAssemblyEnd(ts_B, MAT_FINAL_ASSEMBLY);
546  CHKERR pre_proc_lhs(entX, 0);
547  CHKERR pre_proc_lhs(entY, 1);
548  CHKERR pre_proc_lhs(entZ, 2);
549 
550  break;
551  }
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);
556  CHKERR pre_proc_rhs(entX, 0);
557  CHKERR pre_proc_rhs(entY, 1);
558  CHKERR pre_proc_rhs(entZ, 2);
559 
560  CHKERR VecAssemblyBegin(ts_F);
561  CHKERR VecAssemblyEnd(ts_F);
562 
563  CHKERR VecGhostUpdateBegin(ts_F, ADD_VALUES, SCATTER_REVERSE);
564  CHKERR VecGhostUpdateEnd(ts_F, ADD_VALUES, SCATTER_REVERSE);
565 
566  break;
567  }
568  default:
569  break;
570  }
571 
572  eraseFlag = false;
573  }
575  }
576 
577  MoFEMErrorCode postProcess() { return 0; }
578 
579 private:
580  Range entX;
581  Range entY;
582  Range entZ;
583  bool eraseFlag;
585  std::string problemName;
586 };
587 
588 // }; // namespace OpContactTools
std::map< int, BlockParamData > mat_blocks
static Index< 'l', 3 > l
static Index< 'k', 3 > k
static Index< 'm', 3 > m
static Index< 'M', 3 > M
static Index< 'n', 3 > n
static Index< 'j', 3 > j
static Index< 'o', 3 > o
static Index< 'L', 3 > L
static Index< 'J', 3 > J
static MoFEMErrorCode custom_snes_rhs(SNES snes, Vec x, Vec f, void *ctx)
static Number< 2 > N2
static Index< 'p', 3 > p
BlockParamData * cache
static Index< 'N', 3 > N
static Index< 'I', 3 > I
static Number< 1 > N1
static Index< 'K', 3 > K
static MoFEMErrorCode custom_snes_mat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
static Index< 'i', 3 > i
static Number< 0 > N0
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.
@ TORUS
Definition: RigidBodies.hpp:21
@ STL
Definition: RigidBodies.hpp:23
@ SUPERQUADRIC
Definition: RigidBodies.hpp:22
@ PLANE
Definition: RigidBodies.hpp:18
@ CYLINDER
Definition: RigidBodies.hpp:20
@ NURBS
Definition: RigidBodies.hpp:24
@ SPHERE
Definition: RigidBodies.hpp:19
@ LASTBODY
Definition: RigidBodies.hpp:25
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
Definition: single.cpp:4
@ ROW
Definition: definitions.h:192
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:348
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMMoFEM.cpp:953
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
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.
Definition: SnesCtx.cpp:126
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.
Definition: SnesCtx.cpp:17
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
double cn_cont
double poisson
double young_modulus
Ddg< double, 3, 3 > scaling_mat
static double spring_stiffness
static int closest_roller_id
Ddg< double, 3, 3 > tD
double scale_constraint
double density
static Tensor2< double, 3, 3 > Omega
static RigidBodyData * closest_roller
static double gravity_y
double young_modulus_inv
int getClosestRollerID(Tensor1< T1, 3 > &t_coords, Tensor1< T2, 3 > &t_disp)
double rollers_stop_time
static boost::ptr_vector< RigidBodyData > rollerDataVec
DataFromMove()=delete
DataFromMove(size_t c, double val, Range ents_bc)
VectorDouble scaledValues
std::string problemName
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
bool scaleMat
MoFEMErrorCode preProcess()
MoFEMErrorCode postProcess()
FePrePostProcess(MoFEM::Interface &m_field, boost::ptr_vector< RigidBodyData > &roller_data, boost::ptr_vector< DataFromMove > &bc_data, bool scale_mat)
bool runOnceFlag
boost::ptr_vector< RigidBodyData > & rollerDataVec
boost::ptr_vector< DataFromMove > & bcData
boost::ptr_vector< MethodForForceScaling > methodsOp
LoadScale(double &my_scale)
double & 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