v0.14.0
Loading...
Searching...
No Matches
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
42 static vector<double> gravity;
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
63 virtual MoFEMErrorCode scaleParameters() = 0;
64
65 template <typename T1, typename T2>
67
68 std::vector<double> gaps;
69 Index<'i', 3> i;
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
83static Index<'i', 3> i;
84static Index<'j', 3> j;
85static Index<'k', 3> k;
86static Index<'l', 3> l;
87static Index<'m', 3> m;
88static Index<'n', 3> n;
89static Index<'o', 3> o;
90static Index<'p', 3> p;
91
92static Index<'I', 3> I;
93static Index<'J', 3> J;
94static Index<'K', 3> K;
95static Index<'L', 3> L;
96static Index<'M', 3> M;
97static Index<'N', 3> N;
98
99static Number<0> N0;
100static Number<1> N1;
101static Number<2> N2;
102
104
105 MoFEMErrorCode getOptionsFromCommandLine() {
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
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);
197 int body_type_id = LASTBODY;
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 switch (body_type_id) {
209 case PLANE:
210 rollerDataVec.push_back(new PlaneRigidBody(center_coords, disp, ii));
211 break;
212 case SPHERE:
213 rollerDataVec.push_back(new SphereRigidBody(center_coords, disp, ii));
214 break;
215 case CYLINDER:
216 rollerDataVec.push_back(
217 new CylinderRigidBody(center_coords, disp, ii));
218 break;
219 case CONE:
220 rollerDataVec.push_back(new ConeRigidBody(center_coords, disp, ii));
221 break;
222 case TORUS:
223 rollerDataVec.push_back(new TorusRigidBody(center_coords, disp, ii));
224 break;
225 case ROLLER:
226 rollerDataVec.push_back(new RollerRigidBody(center_coords, disp, ii));
227 break;
228 case STL:
229 rollerDataVec.push_back(new STLRigidBody(center_coords, disp, ii));
230 break;
231 case NURBS:
232 // code to be executed if
233 break;
234 default:
235 // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
236 // "this type of rigid body is not supported, check spelling "
237 // "(plane, sphere, cylinder, torus, superquadric, stl, nurbs)");
238 break;
239 }
240 auto &rol = rollerDataVec.back();
241 CHKERR rol.getBodyOptions();
242
243 char pos_data_file[255];
244 PetscBool ctg_flag = PETSC_FALSE;
245 param = "-position_data" + to_string(ii);
246 CHKERR PetscOptionsString(param.c_str(), "", "", "", pos_data_file, 255,
247 &ctg_flag);
248 if (ctg_flag) {
249 rol.methodOpForRollerPosition =
250 boost::make_shared<TimeAccelerogram>(string(param));
251 rol.originCoords.clear();
252 rol.rollerDisp.clear();
253 }
254
255 param = "-direction_data" + to_string(ii);
256 CHKERR PetscOptionsString(param.c_str(), "", "", "", pos_data_file, 255,
257 &ctg_flag);
258 if (ctg_flag) {
259 rol.methodOpForRollerDirection =
260 boost::make_shared<TimeAccelerogram>(string(param));
261 rol.BodyDirectionScaled = VectorDouble({0, 0, 0});
262 }
263
264 }
265 }
266
267 if (with_contact && rollerDataVec.empty()) {
268 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
269 "provide at least one body data for contact (-body1 "
270 "sphere -radius r and opional -move x,y) ");
271 }
272 // }
273
274 ierr = PetscOptionsEnd();
275 CHKERRQ(ierr);
276
278 }
279
280 MoFEMErrorCode scaleParameters() {
282
284 young_modulus = 1;
290
294 // cn_cont = 1; // perhaps these should be always assigned by user
295 // cn_pl = 1;
296
298 }
299};
300
301extern BlockParamData *cache;
302extern std::map<int, BlockParamData> mat_blocks; // for use in future
303
305 size_t comp;
306 VectorDouble values;
307 VectorDouble scaledValues;
309
310 DataFromMove(size_t c, VectorDouble val, Range ents_bc)
311 : comp(c), values(val), scaledValues(val), ents(ents_bc) {}
312 DataFromMove(VectorDouble val, Range ents_bc)
313 : comp(0), values(val), scaledValues(val), ents(ents_bc) {}
314
315 DataFromMove() = delete;
316};
317
319 double &scale;
320 LoadScale(double &my_scale) : scale(my_scale){};
321 MoFEMErrorCode scaleNf(const FEMethod *fe, VectorDouble &nf) {
323 nf *= scale;
325 }
326};
327
328
329// }; // 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
@ ROLLER
@ STL
@ PLANE
@ CYLINDER
@ NURBS
@ CONE
@ SPHERE
@ LASTBODY
static PetscErrorCode ierr
#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
Definition definitions.h:31
#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 ...
const double c
speed of light (cm/ns)
MoFEMErrorCode scaleParameters()
MoFEMErrorCode getOptionsFromCommandLine()
static Tensor4< double, 3, 3, 3, 3 > Is
virtual MoFEMErrorCode scaleParameters()=0
static vector< double > gravity
double young_modulus
static double spring_stiffness
static int closest_roller_id
static Tensor1< double, 3 > omega
double scale_constraint
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.