v0.13.1
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
19
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;
32 double visH;
34
35 MatrixDouble mtD; // for stiffness matrix
36
37 // for contact
38 double cn_cont;
40
41 // for rotation
43 static double spring_stiffness;
44
46 static double delta;
48 static Tensor2<double, 3, 3> Omega;
49 static Tensor4<double, 3, 3, 3, 3> Is; // symmetric
50
51 static boost::ptr_vector<RigidBodyData> rollerDataVec;
53
56 sigmaY(1e9), cn_pl(1), H(0), C1_k(0), Q_inf(0), b_iso(0), visH(0),
57 cn_cont(1. / young_modulus), /*spring_stiffness(0),
58 closest_roller_id(0),*/
60 mtD.resize(36, 1);
61 }
62
64
65 template <typename T1, typename T2>
67
68 std::vector<double> gaps;
70 for (auto &roller : rollerDataVec) {
71 auto t_normal = roller.getNormal(t_coords, t_disp);
72 double dist = roller.getGap(t_coords);
73 gaps.push_back(dist);
74 }
75
76 auto it = std::min_element(gaps.begin(), gaps.end());
77 auto nb = distance(gaps.begin(), it);
78
79 return nb;
80 };
81};
82
91
98
99static Number<0> N0;
100static Number<1> N1;
101static Number<2> N2;
102
104
107
108 // CHKERR PetscOptionsInsertString(NULL, default_options);
109 CHKERR PetscOptionsInsertString(
110 NULL, "-mat_mumps_icntl_14 1200 -mat_mumps_icntl_24 1 "
111 "-mat_mumps_icntl_13 1 -mat_mumps_icntl_20 0");
112
113 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "none");
114
115 int nb_ent_grav = 3;
116 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-gravity", gravity.data(),
117 &nb_ent_grav, PETSC_NULL);
118 // for elastic
119 CHKERR PetscOptionsScalar("-young", "Young's modulus", "", young_modulus,
120 &young_modulus, PETSC_NULL);
121
122 CHKERR PetscOptionsScalar("-poisson", "Poisson ratio", "", poisson,
123 &poisson, PETSC_NULL);
124 CHKERR PetscOptionsScalar("-density", "Density", "", density, &density,
125 PETSC_NULL);
126
127 PetscBool with_contact = PETSC_FALSE;
128 CHKERR PetscOptionsBool("-with_contact", "run calculations with contact",
129 "", with_contact, &with_contact, PETSC_NULL);
130
131 double rot_vec[3] = {0, 0, 1};
132 int nb_ent = 3;
133 PetscBool is_rotating;
134 // get rotation velocity vector
135 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-rot_omega", rot_vec, &nb_ent,
136 &is_rotating);
137 if (nb_ent < 3 && is_rotating) {
138 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
139 "provide an appropriate number of entries for omega vector (3)");
140 }
141 {
142 Tensor1<double, 3> t1_omega(rot_vec[0], rot_vec[1], rot_vec[2]);
143 omega(i) = t1_omega(i);
144 Omega(i, k) = levi_civita<double>(i, j, k) * t1_omega(j);
145 }
146
147 if (poisson >= 0.5)
148 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
149 "poisson ratio value not supported");
150
151 CHKERR PetscOptionsScalar("-sigmaY", "SigmaY - limit stress", "", sigmaY,
152 &sigmaY, PETSC_NULL);
153 CHKERR PetscOptionsScalar("-H", "Hardening modulus", "", H, &H, PETSC_NULL);
154 CHKERR PetscOptionsScalar("-cn_pl",
155 "cn_pl parameter for active set constrains", "",
156 cn_pl, &cn_pl, PETSC_NULL);
157 CHKERR PetscOptionsScalar("-scale_constraint", "constraint scale", "",
158 scale_constraint, &scale_constraint, PETSC_NULL);
159 CHKERR PetscOptionsScalar("-C1_k", "1st Backstress", "", C1_k, &C1_k,
160 PETSC_NULL);
161 CHKERR PetscOptionsScalar("-Q_inf", "Saturation stress", "", Q_inf, &Q_inf,
162 PETSC_NULL);
163 CHKERR PetscOptionsScalar("-b_iso", "Isotropic exponent", "", b_iso, &b_iso,
164 PETSC_NULL);
165 CHKERR PetscOptionsScalar("-visH", "viscous hardening", "", visH, &visH,
166 PETSC_NULL);
167 // for contact
168 CHKERR PetscOptionsScalar("-spring",
169 "Springs stiffness on the boundary (contact)", "",
170 spring_stiffness, &spring_stiffness, PETSC_NULL);
171
172 CHKERR PetscOptionsScalar("-cn_cont",
173 "cn_cont parameter for active set constrains", "",
174 cn_cont, &cn_cont, PETSC_NULL);
175 CHKERR PetscOptionsScalar(
176 "-rollers_stop", "Time at which rollers will stop moving", "",
178
179 constexpr int MAX_NB_OF_ROLLERS = 10;
180 // if (with_contact) {
181
182 const char *body_list[] = {"plane", "sphere", "cylinder", "cone",
183 "torus", "roller", "stl", "nurbs"};
184 for (int ii = 0; ii < MAX_NB_OF_ROLLERS; ++ii) {
185 int nb_coord = 3;
186 int nb_disp = 3;
187 int nb_dirs = 3;
188 VectorDouble center_coords({0, 0, 0});
189 VectorDouble disp({0, 0, 0});
190 VectorDouble direction({0, 0, 1});
191 PetscBool flg = PETSC_FALSE;
192 PetscBool rflg;
193 string param = "-body" + to_string(ii);
194 // char mesh_file_name[255];
195 // CHKERR PetscOptionsString(param.c_str(), "file name", "", "mesh.vtk",
196 // mesh_file_name, 255, &flg);
198 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, param.c_str(), body_list,
199 LASTBODY, &body_type_id, &flg);
200
201 param = "-coords" + to_string(ii);
202 CHKERR PetscOptionsGetRealArray(NULL, NULL, param.c_str(),
203 &center_coords(0), &nb_coord, &rflg);
204 param = "-move" + to_string(ii);
205 CHKERR PetscOptionsGetRealArray(NULL, NULL, param.c_str(), &disp(0),
206 &nb_disp, &rflg);
207 if (flg) {
208
209 switch (body_type_id) {
210 case PLANE:
211 rollerDataVec.push_back(new PlaneRigidBody(center_coords, disp, ii));
212 break;
213 case SPHERE:
214 rollerDataVec.push_back(new SphereRigidBody(center_coords, disp, ii));
215 break;
216 case CYLINDER:
217 rollerDataVec.push_back(
218 new CylinderRigidBody(center_coords, disp, ii));
219 break;
220 case CONE:
221 rollerDataVec.push_back(new ConeRigidBody(center_coords, disp, ii));
222 break;
223 case TORUS:
224 rollerDataVec.push_back(new TorusRigidBody(center_coords, disp, ii));
225 break;
226 case ROLLER:
227 rollerDataVec.push_back(
228 new RollerRigidBody(center_coords, disp, ii));
229 break;
230 case STL:
231 rollerDataVec.push_back(new STLRigidBody(center_coords, disp, ii));
232 break;
233 case NURBS:
234 // code to be executed if
235 break;
236 default:
237 // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
238 // "this type of rigid body is not supported, check spelling "
239 // "(plane, sphere, cylinder, torus, superquadric, stl, nurbs)");
240 break;
241 }
242 CHKERR rollerDataVec.back().getBodyOptions();
243
244 char pos_data_file[255];
245 PetscBool ctg_flag = PETSC_FALSE;
246 param = "-position_data" + to_string(ii);
247 CHKERR PetscOptionsString(param.c_str(), "", "", "", pos_data_file, 255,
248 &ctg_flag);
249 if (ctg_flag) {
250 rollerDataVec.back().positionDataParamName = string(param);
251 rollerDataVec.back().originCoords.clear();
252 rollerDataVec.back().rollerDisp.clear();
253 }
254 }
255 }
256
257 if (with_contact && rollerDataVec.empty()) {
258 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
259 "provide at least one body data for contact (-body1 "
260 "sphere -radius r and opional -move x,y) ");
261 }
262 // }
263
264 ierr = PetscOptionsEnd();
265 CHKERRQ(ierr);
266
268 }
269
272
274 young_modulus = 1;
280
284 // cn_cont = 1; // perhaps these should be always assigned by user
285 // cn_pl = 1;
286
288 }
289};
290
291extern BlockParamData *cache;
292extern std::map<int, BlockParamData> mat_blocks; // for use in future
293
295 size_t comp;
298 Range ents;
299
300 DataFromMove(size_t c, VectorDouble val, Range ents_bc)
301 : comp(c), values(val), scaledValues(val), ents(ents_bc) {}
302 DataFromMove(VectorDouble val, Range ents_bc)
303 : comp(0), values(val), scaledValues(val), ents(ents_bc) {}
304
305 DataFromMove() = delete;
306};
307
309 double &scale;
310 LoadScale(double &my_scale) : scale(my_scale){};
313 nf *= scale;
315 }
316};
317
318
319// }; // 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 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 Index< 'i', 3 > i
static Number< 0 > N0
@ TORUS
Definition: RigidBodies.hpp:22
@ ROLLER
Definition: RigidBodies.hpp:23
@ STL
Definition: RigidBodies.hpp:24
@ PLANE
Definition: RigidBodies.hpp:18
@ CYLINDER
Definition: RigidBodies.hpp:20
@ NURBS
Definition: RigidBodies.hpp:25
@ CONE
Definition: RigidBodies.hpp:21
@ SPHERE
Definition: RigidBodies.hpp:19
@ LASTBODY
Definition: RigidBodies.hpp:26
Definition: single.cpp:4
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
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
static vector< double > gravity
double cn_cont
double poisson
double young_modulus
static double spring_stiffness
static int closest_roller_id
static Tensor1< double, 3 > omega
double scale_constraint
double density
static Tensor2< double, 3, 3 > Omega
static RigidBodyData * closest_roller
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
MatrixDouble mtD
static double delta
DataFromMove(size_t c, VectorDouble val, Range ents_bc)
DataFromMove()=delete
VectorDouble scaledValues
DataFromMove(VectorDouble val, Range ents_bc)
VectorDouble values
LoadScale(double &my_scale)
double & scale
MoFEMErrorCode scaleNf(const FEMethod *fe, VectorDouble &nf)
Class used to scale loads, f.e. in arc-length control.