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 
25  // for plastic
26  double sigmaY;
27  double cn_pl;
28  // hardening
29  double H;
30  double C1_k;
31  double Q_inf;
32  double b_iso;
33 
34  // for contact
35  double cn_cont;
37  double gravity_y;
40 
41  // for rotation
42  Tensor2<double, 3, 3> Omega;
43  Tensor4<double, 3, 3, 3, 3> Is; // symmetric
44  Tensor4<double, 3, 3, 3, 3> II; // unit
45 
48  sigmaY(1e9), cn_pl(1), H(0), C1_k(0), Q_inf(0), b_iso(0),
50  rollers_stop_time(INFINITY), closest_roller_id(0) {}
51 };
52 
61 
62 static Number<0> N0;
63 static Number<1> N1;
64 static Number<2> N2;
65 
66 struct BlockParamData : public CacheParams {
67 
70 
71  // CHKERR PetscOptionsInsertString(NULL, default_options);
72 
73  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "none");
74 
75  CHKERR PetscOptionsScalar("-gravity_y", "Gravity load in Y direction", "",
76  gravity_y, &gravity_y, PETSC_NULL);
77 
78  // for elastic
79  CHKERR PetscOptionsScalar("-young", "Young's modulus", "", young_modulus,
80  &young_modulus, PETSC_NULL);
81  CHKERR PetscOptionsScalar("-poisson", "Poisson ratio", "", poisson,
82  &poisson, PETSC_NULL);
83  CHKERR PetscOptionsScalar("-density", "Density", "", density, &density,
84  PETSC_NULL);
85 
86  double rot_vec[3] = {0, 0, 1};
87  int nb_ent = 3;
88  PetscBool is_rotating;
89  // get rotation velocity vector
90  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-rot_omega", rot_vec, &nb_ent,
91  &is_rotating);
92  if (nb_ent < 3 && is_rotating) {
93  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
94  "provide an appropriate number of entries for omega vector (3)");
95  }
96  {
97  Tensor1<double, 3> t1_omega;
98  for (int dd : {0, 1, 2})
99  t1_omega(dd) = rot_vec[dd];
100  // Tensor1<double, 3> t1_omega(rot_vec[0], rot_vec[1], rot_vec[2]);
101  Omega(i, k) = levi_civita<double>(i, j, k) * t1_omega(j);
102  }
103 
104  if (poisson >= 0.5)
105  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
106  "poisson ratio value not supported");
107 
108  CHKERR PetscOptionsScalar("-sigmaY", "SigmaY - limit stress", "", sigmaY,
109  &sigmaY, PETSC_NULL);
110  CHKERR PetscOptionsScalar("-H", "Hardening modulus", "", H, &H, PETSC_NULL);
111  CHKERR PetscOptionsScalar("-cn_pl",
112  "cn_pl parameter for active set constrains", "",
113  cn_pl, &cn_pl, PETSC_NULL);
114 
115  CHKERR PetscOptionsScalar("-C1_k", "1st Backstress", "", C1_k, &C1_k,
116  PETSC_NULL);
117  CHKERR PetscOptionsScalar("-Q_inf", "Saturation stress", "", Q_inf, &Q_inf,
118  PETSC_NULL);
119  CHKERR PetscOptionsScalar("-b_iso", "Isotropic exponent", "", b_iso, &b_iso,
120  PETSC_NULL);
121 
122  // for contact
123  CHKERR PetscOptionsScalar("-spring",
124  "Springs stiffness on the boundary (contact)", "",
125  spring_stiffness, &spring_stiffness, PETSC_NULL);
126  cn_cont = 1. / young_modulus;
127 
128  CHKERR PetscOptionsScalar("-cn_cont",
129  "cn_cont parameter for active set constrains", "",
130  cn_cont, &cn_cont, PETSC_NULL);
131  CHKERR PetscOptionsScalar(
132  "-rollers_stop", "Time at which rollers will stop moving", "",
133  rollers_stop_time, &rollers_stop_time, PETSC_NULL);
134 
135  constexpr int MAX_NB_OF_ROLLERS = 6;
136  // if (with_contact) {
137  // for (int i = 0; i < MAX_NB_OF_ROLLERS; ++i) {
138  // int nb_coord = 2;
139  // int nb_disp = 2;
140  // double center_coords[2] = {0, -1};
141  // double disp[2] = {0, 0};
142  // PetscBool flg = PETSC_FALSE;
143  // PetscBool rflg;
144  // string param = "-roller" + to_string(i);
145  // CHKERR PetscOptionsGetRealArray(NULL, NULL, param.c_str(),
146  // center_coords,
147  // &nb_coord, &rflg);
148  // if (rflg) {
149  // param = "-move" + to_string(i);
150  // CHKERR PetscOptionsGetRealArray(NULL, NULL, param.c_str(), disp,
151  // &nb_disp, &flg);
152  // double radi;
153  // param = "-radius" + to_string(i);
154  // CHKERR PetscOptionsScalar(param.c_str(), "set roller radius", "",
155  // radi,
156  // &radi, &flg);
157  // if (nb_coord < 2 || !flg) {
158  // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
159  // "provide an appropriate number of coords for roller
160  // center " "and corresponding radius");
161  // }
162  // VectorDouble center;
163  // VectorDouble displ;
164  // center.data().assign(center_coords, center_coords + 2);
165  // displ.data().assign(disp, disp + 2);
166  // rollerData.push_back(RollerData(center, displ, radi));
167  // }
168  // }
169 
170  // if (rollerData.empty()) {
171  // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
172  // "provide at least one roller data for contact (-roller "
173  // "x,y -radius r and opional -move x,y) ");
174  // }
175  // }
176 
177  ierr = PetscOptionsEnd();
178  CHKERRQ(ierr);
179 
181  }
182 
185 
187  young_modulus = 1;
189  H *= young_modulus_inv;
192 
194  }
195 };
196 
197 extern BlockParamData *cache;
198 // static std::map<int, BlockParamData> cache_blocks; //for use in future
199 
200 struct DataFromMove {
201  size_t comp;
202  double values;
204  Range ents;
205  DataFromMove(size_t c, double val, Range ents_bc)
206  : comp(c), values(val), scaled_values(1), ents(ents_bc) {
207  scaled_values(0) = val;
208  }
209  DataFromMove() = delete;
210 };
211 
212 struct RollerData {
216  const double rAdius;
217  double gAp0;
220 
221  RollerData(VectorDouble c_coords, VectorDouble roller_disp, double radius)
222  : cCoords(c_coords), rollerDisp(roller_disp), rAdius(radius) {
224  rollerDispScaled.clear();
225  }
226 
229  return Tensor1<double, 2>(pos(0), pos(1));
230  };
231 
232  template <typename T1, typename T2>
234  Tensor1<T2, 2> &t_disp) {
235  auto t_center = getRollerCenter();
236  Tensor1<double, 2> t_d;
237  t_d(i) = t_coords(i) - t_center(i) + t_disp(i);
238  double dist = sqrt(t_d(i) * t_d(i));
239  this->t_norm(i) = t_d(i) / dist;
240 
241  // calculate the gap0 (from reference)
242  t_d(i) = t_d(i) - t_disp(i);
243  dist = sqrt(t_d(i) * t_d(i));
244  double gap0 = -t_norm(i) * (t_d(i) - t_norm(i) * rAdius);
245  this->gAp0 = gap0;
246  return t_norm;
247  };
248 
249  inline double getGap0() {
250  // auto norm = getNormal(t_coords);
251  return gAp0;
252  };
253 };
254 
255 template <typename T1, typename T2>
257  vector<RollerData> &rollerData) {
258 
259  std::vector<double> dists;
261 
262  for (auto &roller : rollerData) {
263  // auto t_center = roller.getRollerCenter();
264  // Tensor1<double, 2> t_d;
265  // t_d(i) = t_coords(i) - t_center(i) + t_disp(i);
266  auto t_normal = roller.getNormal(t_coords, t_disp);
267  double dist = t_disp(i) * t_normal(i) - roller.getGap0();
268  // Tensor1<double, 2> t_normal;
269  // t_normal(i) = t_d(i) / dist;
270  // dist = -t_normal(i) * (t_d(i) - t_normal(i) * roller.rAdius);
271  // dist = abs(dist);
272  dists.push_back(dist);
273  }
274  auto it = std::min_element(dists.begin(), dists.end());
275  auto nb = distance(dists.begin(), it);
276 
277  return nb;
278 };
279 
280 /**
281  * \brief Not used at this stage. Could be used to do some calculations,
282  * before assembly of local elements.
283  */
284 struct FePrePostProcess : public FEMethod {
285 
286  vector<RollerData> &rollerData;
287  vector<shared_ptr<DataFromMove>> bcData;
288  ;
289  boost::ptr_vector<MethodForForceScaling> methodsOp;
290 
291  FePrePostProcess(vector<RollerData> &roller_data,
292  vector<shared_ptr<DataFromMove>> bc_data)
293  : rollerData(roller_data), bcData(bc_data) {}
294 
296  //
298 
299  CHKERR VecSetOption(ts_F, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
300  // CHKERR MatSetOption(ts_B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
301 
302  switch (ts_ctx) {
303  case CTX_TSSETIJACOBIAN: {
304  snes_ctx = CTX_SNESSETJACOBIAN;
305  snes_B = ts_B;
306  break;
307  }
308  case CTX_TSSETIFUNCTION: {
309  snes_ctx = CTX_SNESSETFUNCTION;
310  snes_f = ts_F;
311  break;
312  }
313  default:
314  break;
315  }
316 
317  auto time = ts_t;
318 
319  // forceDataScaled = forceData;
320  // pressureDataScaled = pressureData;
321 
322  for (auto &bdata : bcData) {
323  bdata->scaled_values(0) = bdata->values;
325  bdata->scaled_values);
326  }
327 
328  // if (time <= rollers_stop_time) {
329  // if (with_contact)
330  // for (auto &roller : rollerData) {
331  // roller.rollerDispScaled = roller.rollerDisp;
332  // CHKERR MethodForForceScaling::applyScale(this, methodsOp,
333  // roller.rollerDispScaled);
334  // }
335  // }
337  }
338 
339  MoFEMErrorCode postProcess() { return 0; }
340 };
341 
342 struct LoadScale : public MethodForForceScaling {
343  double &scale;
344  LoadScale(double &my_scale) : scale(my_scale){};
347  nf *= scale;
349  }
350 };
351 
352 struct EraseLHSRows : public FEMethod {
353  EraseLHSRows(MoFEM::Interface &m_field, std::string problem_name,
354  Range &essential_boundary_entsX, Range &essential_boundary_entsY,
355  Range &essential_boundary_entsZ)
356  : mField(m_field), entX(essential_boundary_entsX),
357  entY(essential_boundary_entsY), entZ(essential_boundary_entsZ),
358  problemName(problem_name) {}
359 
361  eraseFlag = true;
362  return 0;
363  }
364 
367 
368  auto pre_proc = [&](Range &boundary_ents, int comp) {
370  // if (boundary_ents.empty())
371  // MoFEMFunctionReturnHot(0);
372  IS bc_is;
373 
374  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
375  problemName, ROW, "U", comp, comp, &bc_is, &boundary_ents);
376 
377  // CHKERR ISView(bc_is, PETSC_VIEWER_STDOUT_SELF);
378 
379  CHKERR MatZeroRowsIS(ts_B, bc_is, 0.0, PETSC_NULL, PETSC_NULL);
380 
381  CHKERR ISDestroy(&bc_is);
383  };
384  if (eraseFlag) {
385 
386  CHKERR MatSetOption(ts_B, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
387  CHKERR MatSetOption(ts_B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
388  // CHKERR MatAssemblyBegin(ts_B, MAT_FLUSH_ASSEMBLY);
389  // CHKERR MatAssemblyEnd(ts_B, MAT_FLUSH_ASSEMBLY);
390 
391  CHKERR MatAssemblyBegin(ts_B, MAT_FINAL_ASSEMBLY);
392  CHKERR MatAssemblyEnd(ts_B, MAT_FINAL_ASSEMBLY);
393  CHKERR pre_proc(entX, 0);
394  CHKERR pre_proc(entY, 1);
395  CHKERR pre_proc(entZ, 2);
396  // CHKERR MatAssemblyBegin(ts_B, MAT_FLUSH_ASSEMBLY);
397  // CHKERR MatAssemblyEnd(ts_B, MAT_FLUSH_ASSEMBLY);
398  // CHKERR MatAssemblyBegin(ts_B, MAT_FINAL_ASSEMBLY);
399  // CHKERR MatAssemblyEnd(ts_B, MAT_FINAL_ASSEMBLY);
400  eraseFlag = false;
401  }
403  }
404 
405  MoFEMErrorCode postProcess() { return 0; }
406 
407 private:
408  Range entX;
409  Range entY;
410  Range entZ;
411  bool eraseFlag;
413  std::string problemName;
414 };
415 
416 struct EraseRHSRows : public FEMethod {
417  EraseRHSRows(MoFEM::Interface &m_field, std::string problem_name,
418  Range &essential_boundary_entsX, Range &essential_boundary_entsY,
419  Range &essential_boundary_entsZ)
420  : mField(m_field), entX(essential_boundary_entsX),
421  entY(essential_boundary_entsY), entZ(essential_boundary_entsZ),
422  problemName(problem_name) {}
423 
424  MoFEMErrorCode postProcess() { return 0; }
425 
427  eraseFlag = true;
428  return 0;
429  }
430 
433 
434  auto pre_proc = [&](Range &boundary_ents, int comp) {
436 
437  IS bc_is;
438 
439  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
440  problemName, ROW, "U", comp, comp, &bc_is, &boundary_ents);
441 
442  // CHKERR ISView(bc_is, PETSC_VIEWER_STDOUT_SELF);
443  // CHKERR VecView(ts_F,PETSC_VIEWER_STDOUT_WORLD);
444 
445  CHKERR VecISSet(ts_F, bc_is, 0.0);
446 
447  CHKERR ISDestroy(&bc_is);
449  };
450 
451  if (eraseFlag) {
452  // CHKERR VecGhostUpdateBegin(ts_F, ADD_VALUES, SCATTER_REVERSE);
453  // CHKERR VecGhostUpdateEnd(ts_F, ADD_VALUES, SCATTER_REVERSE);
454  CHKERR VecAssemblyBegin(ts_F);
455  CHKERR VecAssemblyEnd(ts_F);
456  CHKERR pre_proc(entX, 0);
457  CHKERR pre_proc(entY, 1);
458  CHKERR pre_proc(entZ, 2);
459 
460  // CHKERR VecGhostUpdateBegin(d0, INSERT_VALUES, SCATTER_FORWARD);
461  // CHKERR VecGhostUpdateEnd(d0, INSERT_VALUES, SCATTER_FORWARD);
462  // CHKERR DMoFEMPreProcessFiniteElements(sub_dm_rho,
463  // density_bc_ptr.get()); CHKERR VecGhostUpdateBegin(d0, ADD_VALUES,
464  // SCATTER_REVERSE); CHKERR VecGhostUpdateEnd(d0, ADD_VALUES,
465  // SCATTER_REVERSE);
466 
467  CHKERR VecAssemblyBegin(ts_F);
468  CHKERR VecAssemblyEnd(ts_F);
469  CHKERR VecGhostUpdateBegin(ts_F, ADD_VALUES, SCATTER_REVERSE);
470  CHKERR VecGhostUpdateEnd(ts_F, ADD_VALUES, SCATTER_REVERSE);
471 
472  // CHKERR VecAssemblyBegin(ts_F);
473  // CHKERR VecAssemblyEnd(ts_F);
474  // CHKERR VecGhostUpdateBegin(ts_F, ADD_VALUES, SCATTER_FORWARD);
475  // CHKERR VecGhostUpdateEnd(ts_F, ADD_VALUES, SCATTER_FORWARD);
476 
477  eraseFlag = false;
478  }
480  }
481 
482 private:
483  Range entX;
484  Range entY;
485  Range entZ;
486  bool eraseFlag;
488  std::string problemName;
489 };
490 
491 // }; // namespace OpContactTools
double young_modulus_inv
EraseRHSRows(MoFEM::Interface &m_field, std::string problem_name, Range &essential_boundary_entsX, Range &essential_boundary_entsY, Range &essential_boundary_entsZ)
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
double gap0(FTensor::Tensor1< T, 3 > &t_coords, FTensor::Tensor1< double, 2 > &t_normal)
Definition: ContactOps.hpp:143
Deprecated interface functions.
MoFEMErrorCode postProcess()
MoFEM::Interface & mField
std::string problemName
RollerData(VectorDouble c_coords, VectorDouble roller_disp, double radius)
Tensor4< double, 3, 3, 3, 3 > II
boost::ptr_vector< MethodForForceScaling > methodsOp
Tensor2< double, 3, 3 > Omega
Class used to scale loads, f.e. in arc-length control.
LoadScale(double &my_scale)
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
FePrePostProcess(vector< RollerData > &roller_data, vector< shared_ptr< DataFromMove >> bc_data)
static Index< 'p', 3 > p
MoFEMErrorCode scaleParameters()
Index< 'i', 2 > i
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
DataFromMove()=delete
double getGap0()
MoFEMErrorCode preProcess()
double poisson
Tensor1< double, 2 > & getNormal(Tensor1< T1, 3 > &t_coords, Tensor1< T2, 2 > &t_disp)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
static Index< 'l', 3 > l
MoFEM::Interface & mField
MoFEMErrorCode postProcess()
DataFromMove(size_t c, double val, Range ents_bc)
static Index< 'n', 3 > n
vector< shared_ptr< DataFromMove > > bcData
MoFEMErrorCode scaleNf(const FEMethod *fe, VectorDouble &nf)
VectorDouble scaled_values
auto pre_proc
int closest_roller_id
Tensor1< double, 2 > t_norm
static Index< 'm', 3 > m
const VectorDouble rollerDisp
MoFEMErrorCode operator()()
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEMErrorCode preProcess()
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
VectorDouble rollerDispScaled
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
static Index< 'i', 3 > i
Definition: single.cpp:3
const VectorDouble cCoords
static Number< 2 > N2
double density
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
static Index< 'j', 3 > j
double cn_cont
Tensor1< double, 2 > getRollerCenter()
double gravity_y
static Index< 'o', 3 > o
double rollers_stop_time
std::string problemName
int getClosestRollerID(Tensor1< T1, 3 > &t_coords, Tensor1< T2, 2 > &t_disp, vector< RollerData > &rollerData)
#define CHKERR
Inline error check.
Definition: definitions.h:604
static Number< 1 > N1
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
vector< RollerData > & rollerData
double & scale
MoFEMErrorCode postProcess()
MoFEMErrorCode operator()()
static Index< 'k', 3 > k
const double rAdius
EraseLHSRows(MoFEM::Interface &m_field, std::string problem_name, Range &essential_boundary_entsX, Range &essential_boundary_entsY, Range &essential_boundary_entsZ)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
double young_modulus
static Number< 0 > N0
BlockParamData * cache
Implement the scale factor for Newton-Raphson method.
Tensor4< double, 3, 3, 3, 3 > Is
double spring_stiffness
MoFEMErrorCode getOptionsFromCommandLine()
MoFEMErrorCode preProcess()
MatZeroRowsIS(B, vel, 0.0, NULL, NULL)