v0.14.0
Loading...
Searching...
No Matches
multifield_plasticity.cpp
Go to the documentation of this file.
1/**
2 * \file multifield_plasticity.cpp
3 *
4 * Multifield plasticity with contact
5 *
6 */
7
8/* This file is part of MoFEM.
9 * MoFEM is free software: you can redistribute it and/or modify it under
10 * the terms of the GNU Lesser General Public License as published by the
11 * Free Software Foundation, either version 3 of the License, or (at your
12 * option) any later version.
13 *
14 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
15 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 * License for more details.
18 *
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
21
22#include <MoFEM.hpp>
23
24#include <chrono> //for debugging
25#include <boost/tokenizer.hpp>
26#include <IntegrationRules.hpp>
27#include <BlockMatData.hpp>
28#include <MultifieldPlasticity.hpp>
29
30struct MeasureTime {
31 int time;
32 chrono::high_resolution_clock::time_point start;
33 chrono::high_resolution_clock::time_point stop;
34 MeasureTime() { start = chrono::high_resolution_clock::now(); }
35 ~MeasureTime() {
36 stop = chrono::high_resolution_clock::now();
37 auto duration = chrono::duration_cast<chrono::microseconds>(stop - start);
38 cout << "Time taken by function: " << duration.count() << " us." << endl;
39 }
40};
41
42using namespace MoFEM;
43using namespace FTensor;
44
45using OpLogStrainMatrixLhs = FormsIntegrators<DomainEleOp>::Assembly<
46 USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpGradTensorGrad<1, 3, 3, 1>;
47using OpStiffnessMatrixLhs = FormsIntegrators<DomainEleOp>::Assembly<
48 USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpGradSymTensorGrad<1, 3, 3, 0>;
49using OpBodyForce = FormsIntegrators<DomainEleOp>::Assembly<
50 USER_ASSEMBLE>::LinearForm<GAUSS>::OpSource<1, 3>;
51
52using OpMixDivULhs = FormsIntegrators<DomainEleOp>::Assembly<
53 USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMixDivTimesVec<3>;
54using OpLambdaGraULhs = FormsIntegrators<DomainEleOp>::Assembly<
55 USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMixTensorTimesGrad<3>;
56using OpMixDivURhs = FormsIntegrators<DomainEleOp>::Assembly<
57 USER_ASSEMBLE>::LinearForm<GAUSS>::OpMixDivTimesU<3, 3, 3>;
58using OpMiXLambdaGradURhs = FormsIntegrators<DomainEleOp>::Assembly<
59 USER_ASSEMBLE>::LinearForm<GAUSS>::OpMixTensorTimesGradU<3>;
60using OpMixUTimesDivLambdaRhs = FormsIntegrators<DomainEleOp>::Assembly<
61 USER_ASSEMBLE>::LinearForm<GAUSS>::OpMixVecTimesDivLambda<3>;
62using OpMixUTimesLambdaRhs = FormsIntegrators<DomainEleOp>::Assembly<
63 USER_ASSEMBLE>::LinearForm<GAUSS>::OpGradTimesTensor<1, 3, 3>;
64
65using OpSpringLhs = FormsIntegrators<BoundaryEleOp>::Assembly<
66 USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<1, 3>;
67using OpSpringLhsNoFS = FormsIntegrators<BoundaryEleOp>::Assembly<
68 USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<1, 3>;
69using OpSpringRhs = FormsIntegrators<BoundaryEleOp>::Assembly<
70 USER_ASSEMBLE>::LinearForm<GAUSS>::OpBaseTimesVector<1, 3, 1>;
71
72using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
73 USER_ASSEMBLE>::LinearForm<GAUSS>::OpBaseTimesVector<1, 3, 1>;
74
75using OpCentrifugalForce2 = FormsIntegrators<DomainEleOp>::Assembly<
76 USER_ASSEMBLE>::LinearForm<GAUSS>::OpSource<1, 3>;
77
78// using OpScaleL2 = MoFEM::OpScaleBaseBySpaceInverseOfMeasure<DomainEleOp>;
79
80#include <BasicFiniteElements.hpp>
81#include <AnalyticalSurfaces.hpp>
82#include <RigidBodies.hpp>
83#include <BasicFeTools.hpp>
84#include <MatrixFunction.hpp>
85
86// global variables
87BlockParamData *cache;
88std::map<int, BlockParamData> mat_blocks;
89boost::weak_ptr<MoFEM::BlockMatContainer> MoFEM::block_mat_container_ptr;
90
91#include <ElasticOperators.hpp>
92#include <PlasticOperators.hpp>
93#include <ContactOperators.hpp>
94#include <RotatingFrameOperators.hpp>
95#include <DualBase.hpp>
96
97using namespace OpElasticTools;
98using namespace OpPlasticTools;
99using namespace OpContactTools;
100using namespace OpRotatingFrameTools;
101
102// for testing
103constexpr bool TEST_H1_SPACE = false;
104EntityType zero_type = TEST_H1_SPACE ? MBVERTEX : MBTET;
105FieldSpace space_test = (TEST_H1_SPACE ? H1 : L2);
106
107struct ContactPlasticity {
108
109 ContactPlasticity(MoFEM::Interface &m_field) : mField(m_field) {}
110
111 MoFEMErrorCode runProblem();
112 MoFEMErrorCode loadMeshBlockSets();
113 static MoFEMErrorCode getAnalysisParameters();
114 static MoFEMErrorCode myMeshPartition(MoFEM::Interface &m_field);
115
116private:
117 MoFEM::Interface &mField;
118 MoFEMErrorCode setUP();
119 MoFEMErrorCode createCommonData();
120 MoFEMErrorCode OPs();
121 MoFEMErrorCode bC();
122 MoFEMErrorCode tsSolve();
123 MoFEMErrorCode postProcess();
124 MoFEMErrorCode checkResults();
125
126 MatrixDouble invJac, jAc;
127 boost::shared_ptr<OpContactTools::CommonData> commonDataPtr;
128 boost::shared_ptr<PostProcVolumeOnRefinedMesh> postProcFe;
129 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
130 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
131 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
132
133 // mofem boundary conditions
134 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
135 boost::ptr_map<std::string, NodalForce> nodal_forces;
136 boost::ptr_map<std::string, EdgeForce> edge_forces;
137 boost::shared_ptr<DirichletDisplacementRemoveDofsBc> dirichlet_bc_ptr;
138
139 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
140 // for forces scaling
141 boost::ptr_vector<DataFromMove> bcData;
142 boost::shared_ptr<DomainSideEle> volSideFe;
143
144 // boost::shared_ptr<DomainEle> contactVolFeLhs;
145 // boost::shared_ptr<DomainEle> contactVolFeRhs;
146 boost::shared_ptr<DomainEle> feLambdaRhs;
147 boost::shared_ptr<DomainEle> feLambdaLhs;
148
149 Range getEntsOnMeshSkin();
150
151 // SmartPetscObj<DM> dmElastic;
152 SmartPetscObj<DM> dmEP;
153 SmartPetscObj<DM> dmTAU;
154 SmartPetscObj<DM> dmContact;
155
156public:
157 FieldApproximationBase base;
158 // global params
159 struct ProblemData {
160 int order;
161 PetscBool no_output;
162 PetscBool debug_info;
163 int output_freq;
164 int atom_test_nb;
165 int reaction_id;
166 PetscBool scale_params;
167
168 PetscBool is_large_strain;
169 PetscBool is_fieldsplit_bc_only;
170
171 PetscBool is_neohooke;
172 PetscBool with_plasticity;
173 PetscBool calculate_reactions;
174 PetscBool is_ho_geometry;
175 PetscBool save_restarts;
176 PetscBool is_restart;
177
178 PetscBool is_reduced_integration;
179 PetscBool test_user_base_tau;
180 PetscBool is_alm;
181
182 PetscBool with_contact;
183 PetscBool print_rollers; //
184 PetscBool is_rotating;
185
186 // PetscBool use_fieldsplit_for_bc;
187
188 string contact_history;
189 string move_history;
190 string dirichlet_history;
191 string mesh_file_str;
192 string restart_file_str;
193
194 boost::shared_ptr<BoundaryEle> update_roller;
195 boost::shared_ptr<DomainEle> calc_reactions;
196 boost::shared_ptr<DomainEle> contact_pipeline;
197
198 moab::Core mb_post_debug;
199 moab::Interface *moab_debug;
200
201 int restart_step;
202
203 Range skinEnts;
204 std::map<EntityHandle, int> mapBlockTets;
205 ProblemData() {
206 order = 2;
207 no_output = PETSC_FALSE;
208 debug_info = PETSC_FALSE;
209 output_freq = 1;
210 atom_test_nb = 0;
211 scale_params = PETSC_FALSE;
212
213 is_large_strain = PETSC_FALSE;
214 is_fieldsplit_bc_only = PETSC_FALSE;
215
216 is_reduced_integration = PETSC_FALSE;
217 test_user_base_tau = PETSC_FALSE;
218 is_alm = PETSC_FALSE;
219
220 is_neohooke = PETSC_FALSE;
221 with_plasticity = PETSC_FALSE;
222
223 calculate_reactions = PETSC_FALSE;
224 reaction_id = -1;
225
226 with_contact = PETSC_FALSE;
227 print_rollers = PETSC_FALSE; //
228 is_rotating = PETSC_FALSE;
229
230 // use_fieldsplit_for_bc = PETSC_TRUE;
231
232 is_ho_geometry = PETSC_FALSE;
233
234 save_restarts = PETSC_FALSE;
235 is_restart = PETSC_FALSE;
236
237 restart_step = -1;
238
239 contact_history = move_history = dirichlet_history = "-load_history";
240 moab_debug = &mb_post_debug;
241 };
242 };
243
244 static ProblemData data;
245
246 struct MMonitor : public FEMethod {
247
248 MMonitor(
249 SmartPetscObj<DM> &dm, MoFEM::Interface &m_field,
250 boost::shared_ptr<OpContactTools::CommonData> common_data_ptr,
251 boost::shared_ptr<PostProcVolumeOnRefinedMesh> post_proc_fe,
252 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
253 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter,
254 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uz_scatter)
255 : dM(dm), mField(m_field), commonDataPtr(common_data_ptr),
256 postProcFe(post_proc_fe), uXScatter(ux_scatter),
257 uYScatter(uy_scatter), uZScatter(uz_scatter) {
258 volSideFe = boost::make_shared<MoFEM::DomainSideEle>(mField);
259 volSideFe->getOpPtrVector().push_back(
260 new OpVolumeSideCalculateTAU("TAU", commonDataPtr));
261 volSideFe->getOpPtrVector().push_back(
262 new OpVolumeSideCalculateEP("EP", commonDataPtr));
263 postProcFeSkeleton =
264 boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
265 postProcFeSkeleton->generateReferenceElementMesh();
266 postProcFeSkeleton->getOpPtrVector().push_back(
267 new OpCalculateJumpOnSkeleton("SKELETON_LAMBDA", commonDataPtr,
268 volSideFe));
269 postProcFeSkeleton->getOpPtrVector().push_back(
270 new OpPostProcPlasticSkeleton(
271 "SKELETON_LAMBDA", postProcFeSkeleton->postProcMesh,
272 postProcFeSkeleton->mapGaussPts, commonDataPtr));
273
274 postProcFeSkin = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
275 CHKERR postProcFeSkin->generateReferenceElementMesh();
276 CHKERR postProcFeSkin->addFieldValuesPostProc("U", "DISPLACEMENT");
277 if (data.test_user_base_tau)
278 postProcFeSkin->getUserPolynomialBase() =
279 boost::shared_ptr<BaseFunction>(new BubbleTauPolynomialBase());
280
281 if (data.is_ho_geometry)
282 postProcFeSkin->addFieldValuesPostProc("MESH_NODE_POSITIONS");
283
284 if (data.with_plasticity) {
285 CHKERR postProcFeSkin->addFieldValuesPostProcOnSkin("TAU");
286 CHKERR postProcFeSkin->addFieldValuesPostProcOnSkin("EP");
287 }
288
289 if (data.with_contact) {
290 CHKERR postProcFeSkin->addFieldValuesPostProcOnSkin("SIGMA");
291 }
292 };
293
294 MoFEMErrorCode preProcess() {
295 MoFEMFunctionBegin;
296 CHKERR TSGetTimeStep(ts, &(commonDataPtr->timeStepSize));
297 if (!ts_step || ts_step == data.restart_step) {
298 for (auto &roller : (*cache).rollerDataVec) {
299 if (roller.methodOpForRollerPosition) {
300 roller.BodyDispScaled.clear();
301 CHKERR MethodForForceScaling::applyScale(
302 this, roller.methodOpForRollerPosition, roller.BodyDispScaled);
303 }
304 if (roller.methodOpForRollerDirection) {
305 roller.BodyDirectionScaled.clear();
306 CHKERR MethodForForceScaling::applyScale(
307 this, roller.methodOpForRollerDirection,
308 roller.BodyDirectionScaled);
309 }
310 }
311 }
312
313 MoFEMFunctionReturn(0);
314 }
315 MoFEMErrorCode operator()() { return 0; }
316
317 MoFEMErrorCode postProcess() {
318 MoFEMFunctionBegin;
319
320 auto make_vtk = [&]() {
321 MoFEMFunctionBegin;
322 if (data.no_output || ts_step % data.output_freq != 0)
323 MoFEMFunctionReturnHot(0);
324 CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProcFe);
325 CHKERR postProcFe->writeFile("out_plastic_" +
326 boost::lexical_cast<std::string>(ts_step) +
327 ".h5m");
328 MoFEMFunctionReturn(0);
329 };
330
331 auto save_skeleton = [&]() {
332 MoFEMFunctionBegin;
333
334 CHKERR DMoFEMLoopFiniteElements(dM, "sFE", postProcFeSkeleton);
335 CHKERR postProcFeSkeleton->writeFile(
336 "out_skeleton_" + boost::lexical_cast<std::string>(ts_step) +
337 ".h5m");
338 MoFEMFunctionReturn(0);
339 };
340
341 auto save_skin = [&]() {
342 MoFEMFunctionBegin;
343 if (data.no_output || ts_step % data.output_freq != 0)
344 MoFEMFunctionReturnHot(0);
345
346 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcFeSkin);
347 CHKERR postProcFeSkin->writeFile(
348 "out_skin_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
349 MoFEMFunctionReturn(0);
350 };
351
352 double max_disp, min_disp;
353
354 auto print_max_min = [&](auto &tuple, const std::string msg) {
355 MoFEMFunctionBegin;
356 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
357 INSERT_VALUES, SCATTER_FORWARD);
358 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
359 INSERT_VALUES, SCATTER_FORWARD);
360 CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max_disp);
361 CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min_disp);
362 PetscPrintf(PETSC_COMM_WORLD, "%s time %3.4e min %3.4e max %3.4e\n",
363 msg.c_str(), ts_t, min_disp, max_disp);
364 MoFEMFunctionReturn(0);
365 };
366
367 double gauss_area;
368
369 if (data.with_contact) {
370 auto &l_reactions = commonDataPtr->bodyReactions;
371 auto &lgauss_pts = commonDataPtr->stateArrayArea;
372
373 int roller_nb = (*cache).rollerDataVec.size();
374 l_reactions.resize(roller_nb * 3);
375 std::fill(l_reactions.begin(), l_reactions.end(), 0);
376
377 CHKERR DMoFEMLoopFiniteElements(dM, "bFE", data.update_roller);
378
379 if (true) {
380 std::vector<double> gauss_pts(2, 0);
381 auto dm = mField.getInterface<Simple>()->getDM();
382 CHKERR MPIU_Allreduce(lgauss_pts.data(), gauss_pts.data(), 2,
383 MPIU_REAL, MPIU_SUM,
384 PetscObjectComm((PetscObject)dm));
385 MOFEM_LOG_C("WORLD", Sev::inform,
386 "\t TS %d Total contact area (ref): %g / %g", ts_step,
387 gauss_pts[0], gauss_pts[1]);
388 gauss_area = gauss_pts[0];
389 lgauss_pts[0] = lgauss_pts[1] = 0;
390
391 vector<double> g_reactions(l_reactions.size() * 3, 0);
392 CHKERR MPIU_Allreduce(l_reactions.data(), g_reactions.data(),
393 roller_nb * 3, MPIU_REAL, MPIU_SUM,
394 PetscObjectComm((PetscObject)dm));
395 for (int dd = 0; dd != roller_nb; ++dd)
396 MOFEM_LOG_C("WORLD", Sev::inform,
397 "\t Body %d reactions Time %g Force X: %3.4e Y: "
398 "%3.4e Z: %3.4e",
399 (*cache).rollerDataVec[dd].iD, ts_t,
400 g_reactions[3 * dd + 0], g_reactions[3 * dd + 1],
401 g_reactions[3 * dd + 2]);
402 }
403 }
404
405 if (ts_step % data.output_freq == 0 && !data.no_output) {
406 std::ostringstream ostrm, ostrm2;
407 ostrm << "out_debug_" << ts_step << ".vtk";
408 ostrm2 << "out_debug_" << ts_step << ".h5m";
409 if (data.with_contact && data.debug_info) {
410 if (mField.get_comm_size() == 1)
411 CHKERR data.moab_debug->write_file(ostrm.str().c_str(), "VTK", "");
412 else
413 CHKERR data.moab_debug->write_file(ostrm2.str().c_str(), "MOAB",
414 "PARALLEL=WRITE_PART");
415
416 data.moab_debug->delete_mesh();
417 }
418 }
419
420 if (data.calculate_reactions) {
421 auto &vec = commonDataPtr->reactionsVec;
422 CHKERR VecZeroEntries(vec);
423 CHKERR VecGhostUpdateBegin(vec, INSERT_VALUES, SCATTER_FORWARD);
424 CHKERR VecGhostUpdateEnd(vec, INSERT_VALUES, SCATTER_FORWARD);
425 CHKERR DMoFEMLoopFiniteElements(dM, "dFE", data.calc_reactions);
426 Simple *simple = mField.getInterface<Simple>();
427 Range ents;
428 CHKERR mField.get_moab().get_entities_by_dimension(0, 0, ents, true);
429
430 ParallelComm *pcomm =
431 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
432 CHKERR pcomm->reduce_tags(commonDataPtr->reactionTag, MPI_SUM, ents);
433
434 CHKERR make_vtk();
435 CHKERR save_skin();
436
437 // clear reactions tag
438 double def_value = 0.;
439 CHKERR mField.get_moab().tag_clear_data(commonDataPtr->reactionTag,
440 ents, &def_value);
441
442 // CHKERR VecGhostUpdateBegin(vec, ADD_VALUES, SCATTER_REVERSE);
443 // CHKERR VecGhostUpdateEnd(vec, ADD_VALUES, SCATTER_REVERSE);
444 CHKERR VecAssemblyBegin(vec);
445 CHKERR VecAssemblyEnd(vec);
446 CHKERR VecGhostUpdateBegin(vec, ADD_VALUES, SCATTER_REVERSE);
447 CHKERR VecGhostUpdateEnd(vec, ADD_VALUES, SCATTER_REVERSE);
448 CHKERR VecGhostUpdateBegin(vec, INSERT_VALUES, SCATTER_FORWARD);
449 CHKERR VecGhostUpdateEnd(vec, INSERT_VALUES, SCATTER_FORWARD);
450
451 const double *react_ptr;
452 CHKERR VecGetArrayRead(vec, &react_ptr);
453
454 PetscPrintf(PETSC_COMM_WORLD,
455 "Reactions time %3.4e X: %3.4e Y: %3.4e Z: %3.4e \n", ts_t,
456 react_ptr[0], react_ptr[1], react_ptr[2]);
457
458 CHKERR VecRestoreArrayRead(vec, &react_ptr);
459 } else {
460 CHKERR make_vtk();
461 CHKERR save_skin();
462 }
463 // FIXME: this does not work properly
464 if (data.debug_info && data.is_rotating && data.with_plasticity)
465 CHKERR save_skeleton();
466
467 if (data.with_contact && data.print_rollers &&
468 ts_step % data.output_freq == 0 && !data.no_output) {
469 auto &m_field = postProcFe->mField;
470 data.moab_debug->delete_mesh();
471 if (m_field.get_comm_rank() == 0) {
472 std::ostringstream ost;
473 for (auto &rol : (*cache).rollerDataVec) {
474 EntityHandle vertex;
475 VectorDouble roller_coords = rol.BodyDispScaled;
476 roller_coords += rol.originCoords;
477
478 CHKERR data.moab_debug->create_vertex(&roller_coords(0), vertex);
479 CHKERR rol.saveBasicDataOnTag(*data.moab_debug, vertex);
480 }
481
482 ost << "out_roller_" << ts_step << ".vtk";
483 CHKERR data.moab_debug->write_file(ost.str().c_str(), "VTK", "");
484 data.moab_debug->delete_mesh();
485 }
486 }
487
488 if (data.save_restarts && ts_step % data.output_freq == 0) {
489
490 CHKERR VecGhostUpdateBegin(ts_u, INSERT_VALUES, SCATTER_FORWARD);
491 CHKERR VecGhostUpdateEnd(ts_u, INSERT_VALUES, SCATTER_FORWARD);
492
493 string rest_name =
494 "restart_" + to_string(ts_step) + "_" + to_string(ts_t) + ".dat";
495 PetscViewer viewer;
496 CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, rest_name.c_str(),
497 FILE_MODE_WRITE, &viewer);
498 CHKERR PetscViewerBinarySetSkipInfo(viewer, PETSC_TRUE);
499 VecView(ts_u, viewer);
500 PetscViewerDestroy(&viewer);
501 }
502
503 CHKERR print_max_min(uXScatter, "Ux");
504 CHKERR print_max_min(uYScatter, "Uy");
505 CHKERR print_max_min(uZScatter, "Uz");
506
507 switch (data.atom_test_nb) {
508 case 1: {
509 if (ts_step == 10) {
510 const double *react_ptr;
511 CHKERR VecGetArrayRead(commonDataPtr->reactionsVec, &react_ptr);
512 if (fabs(react_ptr[0] + 0.12519621572) > 1e-6)
513 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
514 "atom test diverged!");
515 CHKERR VecRestoreArrayRead(commonDataPtr->reactionsVec, &react_ptr);
516 }
517 break;
518 }
519 case 2: {
520 if (ts_step == 5) {
521 const double *react_ptr;
522 CHKERR VecGetArrayRead(commonDataPtr->reactionsVec, &react_ptr);
523 // hertz solution 0.25*(4/3)*(1/(1-0.3*0.3))*(10^0.5)*0.1^(3/2)
524 if ((0.03663003663 - react_ptr[2]) / 0.03663003663 > 1e-2)
525 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
526 "atom test diverged!");
527 CHKERR VecRestoreArrayRead(commonDataPtr->reactionsVec, &react_ptr);
528 }
529 break;
530 }
531 case 3: {
532 if (ts_t == 1) {
533 double min_disp;
534 CHKERR VecMin(std::get<0>(uXScatter), PETSC_NULL, &min_disp);
535 if (fabs((min_disp + 2.70) / 2.70) > 5e-2)
536 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
537 "atom test diverged!");
538 }
539 break;
540 }
541 case 4: {
542 if (ts_step == 10) {
543 const double *react_ptr;
544 CHKERR VecGetArrayRead(commonDataPtr->reactionsVec, &react_ptr);
545 // Reaction Force at first yield Johnson = Contact Mechanics = Eq6.10
546 // Fy = (yield_stress/Em)^2 * yield_stress *(Pi^3 * 1.6^3 * R^2)/6
547 // Em = E /(1-poisson_ratio^2)
548 if (abs(0.056090728 - (react_ptr[2] * 4.0)) / 0.056090728 > 5.0e-2)
549 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
550 "atom test diverged!");
551 }
552 case 5: {
553 if (ts_step == 14) {
554 const double *react_ptr;
555 //comparison with data from miehe.dat file
556 CHKERR VecGetArrayRead(commonDataPtr->reactionsVec, &react_ptr);
557 if (abs(19030.7 + react_ptr[2]) / 19030.7 > 1e-3)
558 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
559 "atom test diverged!");
560 }
561 break;
562 }
563 break;
564 }
565
566 default:
567 break;
568 }
569
570 MoFEMFunctionReturn(0);
571 }
572
573 private:
574 SmartPetscObj<DM> dM;
575 MoFEM::Interface &mField;
576 boost::shared_ptr<OpContactTools::CommonData> commonDataPtr;
577 boost::shared_ptr<PostProcVolumeOnRefinedMesh> postProcFe;
578
579 boost::shared_ptr<DomainSideEle> volSideFe;
580 boost::shared_ptr<PostProcFaceOnRefinedMesh> postProcFeSkeleton;
581 boost::shared_ptr<PostProcFaceOnRefinedMesh> postProcFeSkin;
582
583 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
584 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
585 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
586 };
587};
588
589ContactPlasticity::ProblemData ContactPlasticity::data;
590
591MoFEMErrorCode ContactPlasticity::runProblem() {
592 MoFEMFunctionBegin;
593 CHKERR setUP();
594 CHKERR createCommonData();
595 CHKERR bC();
596 CHKERR OPs();
597 CHKERR tsSolve();
598 CHKERR postProcess();
599 CHKERR checkResults();
600 MoFEMFunctionReturn(0);
601}
602MoFEMErrorCode ContactPlasticity::getAnalysisParameters() {
603 MoFEMFunctionBegin;
604
605 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "none");
606 // Set approximation order
607 CHKERR PetscOptionsInt("-order", "approximation order", "", data.order,
608 &data.order, PETSC_NULL);
609 CHKERR PetscOptionsBool("-no_output", "save no output files", "",
610 data.no_output, &data.no_output, PETSC_NULL);
611 CHKERR PetscOptionsBool("-debug_info", "print debug info", "",
612 data.debug_info, &data.debug_info, PETSC_NULL);
613 CHKERR PetscOptionsBool("-scale_params",
614 "scale parameters (young modulus = 1)", "",
615 data.scale_params, &data.scale_params, PETSC_NULL);
616 CHKERR PetscOptionsInt("-output_every",
617 "frequncy how often results are dumped on hard disk",
618 "", 1, &data.output_freq, PETSC_NULL);
619 CHKERR PetscOptionsBool("-calculate_reactions", "calculate reactions", "",
620 data.calculate_reactions, &data.calculate_reactions,
621 PETSC_NULL);
622 CHKERR PetscOptionsInt("-reaction_id", "meshset id for computing reactions",
623 "", data.reaction_id, &data.reaction_id, PETSC_NULL);
624 CHKERR PetscOptionsInt("-atom_test", "number of atom test", "",
625 data.atom_test_nb, &data.atom_test_nb, PETSC_NULL);
626 CHKERR PetscOptionsBool("-is_large_strain", "is large strains regime", "",
627 data.is_large_strain, &data.is_large_strain,
628 PETSC_NULL);
629 CHKERR PetscOptionsBool(
630 "-is_fieldsplit_bc_only", "use fieldsplit for boundary only", "",
631 data.is_fieldsplit_bc_only, &data.is_fieldsplit_bc_only, PETSC_NULL);
632 CHKERR PetscOptionsBool(
633 "-is_reduced_integration",
634 "use axiator split approach, with reduced integration", "",
635 data.is_reduced_integration, &data.is_reduced_integration, PETSC_NULL);
636 CHKERR PetscOptionsBool("-test_user_base_tau", "test_user_base_for_tau", "",
637 data.test_user_base_tau, &data.test_user_base_tau,
638 PETSC_NULL);
639 CHKERR PetscOptionsBool("-is_alm",
640 "use Augmented Lagrangian method for enforcing the "
641 "KKT condtions for plasticity",
642 "", data.is_alm, &data.is_alm, PETSC_NULL);
643
644 CHKERR PetscOptionsBool("-is_neohooke", "is neohooke material", "",
645 data.is_neohooke, &data.is_neohooke, PETSC_NULL);
646 // CHKERR PetscOptionsBool(
647 // "-use_fieldsplit_for_bc", "use fieldsplit for boundary conditions", "",
648 // data.use_fieldsplit_for_bc, &data.use_fieldsplit_for_bc, PETSC_NULL);
649
650 CHKERR PetscOptionsBool(
651 "-with_plasticity", "run calculations with plasticity", "",
652 data.with_plasticity, &data.with_plasticity, PETSC_NULL);
653
654 CHKERR PetscOptionsBool("-is_rotating", "is rotating frame used", "",
655 data.is_rotating, &data.is_rotating, PETSC_NULL);
656 CHKERR PetscOptionsBool("-with_contact", "run calculations with contact", "",
657 data.with_contact, &data.with_contact, PETSC_NULL);
658 CHKERR PetscOptionsBool("-print_rollers",
659 "output file with rollers positions", "",
660 data.print_rollers, &data.print_rollers, PETSC_NULL);
661
662 char load_hist_file[255];
663 PetscBool ctg_flag = PETSC_FALSE;
664 CHKERR PetscOptionsString("-contact_history", "load_history.in", "",
665 "load_history.in", load_hist_file, 255, &ctg_flag);
666 if (ctg_flag)
667 data.contact_history = "-contact_history";
668 CHKERR PetscOptionsString("-move_history", "load_history.in", "",
669 "load_history.in", load_hist_file, 255, &ctg_flag);
670
671 if (ctg_flag)
672 data.move_history = "-move_history";
673
674 CHKERR PetscOptionsString("-dirichlet_history", "load_history.in", "",
675 "load_history.in", load_hist_file, 255, &ctg_flag);
676 if (ctg_flag)
677 data.dirichlet_history = "-dirichlet_history";
678
679 char mesh_file_name[255];
680 CHKERR PetscOptionsString("-file_name", "file name", "", "mesh.h5m",
681 mesh_file_name, 255, PETSC_NULL);
682 data.mesh_file_str = string(mesh_file_name);
683
684 char restart_file_name[255];
685 CHKERR PetscOptionsString("-restart_file", "restart file name", "",
686 "restart.dat", restart_file_name, 255,
687 &data.is_restart);
688 data.restart_file_str = string(restart_file_name);
689
690 CHKERR PetscOptionsBool("-save_restarts",
691 "save restart files at each iteration", "",
692 data.save_restarts, &data.save_restarts, PETSC_NULL);
693
694 ierr = PetscOptionsEnd();
695 CHKERRQ(ierr);
696
697 if (data.is_large_strain && data.with_plasticity && data.is_neohooke)
698 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
699 "Neohooke material is not supported with plasticity.");
700
701 MoFEMFunctionReturn(0);
702}
703
704MoFEMErrorCode ContactPlasticity::myMeshPartition(MoFEM::Interface &m_field) {
705 MoFEMFunctionBegin;
706
707 auto &moab = m_field.get_moab();
708 Range tets;
709 CHKERR moab.get_entities_by_dimension(0, 3, tets, false);
710 Range faces;
711 CHKERR moab.get_adjacencies(tets, 2, true, faces, moab::Interface::UNION);
712 Range edges;
713 CHKERR moab.get_adjacencies(tets, 1, true, edges, moab::Interface::UNION);
714
715 CHKERR m_field.getInterface<ProblemsManager>()->partitionMesh(
716 tets, 3, 0, m_field.get_comm_size());
717
718 EntityHandle part_set;
719 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
720 if (pcomm == NULL)
721 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
722 CHKERR moab.create_meshset(MESHSET_SET, part_set);
723 Tag part_tag = pcomm->part_tag();
724 Range proc_ents;
725 Range all_proc_ents;
726 Range off_proc_ents;
727 Range tagged_sets;
728 CHKERR moab.get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag, NULL, 1,
729 tagged_sets, moab::Interface::UNION);
730
731 for (auto &mit : tagged_sets) {
732 int part;
733 CHKERR moab.tag_get_data(part_tag, &mit, 1, &part);
734 if (part == m_field.get_comm_rank()) {
735 CHKERR moab.get_entities_by_dimension(mit, 3, proc_ents, true);
736 CHKERR moab.get_entities_by_handle(mit, all_proc_ents, true);
737 } else
738 CHKERR moab.get_entities_by_handle(mit, off_proc_ents, true);
739 }
740
741 Skinner skin(&moab);
742 Range proc_ents_skin[4];
743 proc_ents_skin[3] = proc_ents;
744
745 Range all_tets, all_skin;
746 CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, all_tets, false);
747 CHKERR skin.find_skin(0, all_tets, false, all_skin);
748 CHKERR skin.find_skin(0, proc_ents, false, proc_ents_skin[2]);
749
750 proc_ents_skin[2] = subtract(proc_ents_skin[2], all_skin);
751 CHKERR moab.get_adjacencies(proc_ents_skin[2], 1, false, proc_ents_skin[1],
752 moab::Interface::UNION);
753 CHKERR moab.get_connectivity(proc_ents_skin[1], proc_ents_skin[0], true);
754 for (int dd = 0; dd != 3; dd++)
755 CHKERR moab.add_entities(part_set, proc_ents_skin[dd]);
756 CHKERR moab.add_entities(part_set, all_proc_ents);
757
758 // not sure if necessary...
759 {
760 Range all_ents;
761 CHKERR moab.get_entities_by_handle(0, all_ents);
762 std::vector<unsigned char> pstat_tag(all_ents.size(), 0);
763 CHKERR moab.tag_set_data(pcomm->pstatus_tag(), all_ents,
764 &*pstat_tag.begin());
765 }
766
767 Range ents_to_keep;
768 CHKERR moab.get_entities_by_handle(part_set, ents_to_keep, false);
769 off_proc_ents = subtract(off_proc_ents, ents_to_keep);
770 Range meshsets;
771 CHKERR moab.get_entities_by_type(0, MBENTITYSET, meshsets, false);
772 for (auto m : meshsets) {
773 CHKERR moab.remove_entities(m, off_proc_ents);
774 // CHKERR moab.add_entities(part_set, &m, 1);
775 }
776 CHKERR moab.delete_entities(off_proc_ents);
777
778 CHKERR pcomm->resolve_shared_ents(0, proc_ents, 3, -1, proc_ents_skin);
779
780 BitRefLevel bit_level0;
781 bit_level0.set(0);
782 Simple *simple = m_field.getInterface<Simple>();
783
784 Range ents;
785 CHKERR m_field.get_moab().get_entities_by_dimension(part_set, 3, ents, true);
786 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
787 ents, simple->getBitRefLevel(), false);
788 simple->getMeshSet() = part_set;
789
790 MoFEMFunctionReturn(0);
791}
792
793MoFEMErrorCode ContactPlasticity::loadMeshBlockSets() {
794 MoFEMFunctionBegin;
795 string block_name = "MAT_ELASTIC";
796 // TODO: implement blocks for plasticity
797 if (!data.with_plasticity)
798 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
799 if (it->getName().compare(0, block_name.length(), block_name) == 0) {
800 std::vector<double> block_data;
801 CHKERR it->getAttributes(block_data);
802 const int id = it->getMeshsetId();
803 EntityHandle meshset = it->getMeshset();
804 Range tets;
805 // on a particular part, the meshset can be empty
806 rval =
807 mField.get_moab().get_entities_by_dimension(meshset, 3, tets, true);
808 rval = MB_SUCCESS;
809 for (auto ent : tets)
810 data.mapBlockTets[ent] = id;
811
812 mat_blocks[id] = BlockParamData();
813 mat_blocks.at(id).young_modulus = block_data[0];
814 mat_blocks.at(id).poisson = block_data[1];
815 if (block_data.size() > 2)
816 mat_blocks.at(id).density = block_data[2];
817 mat_blocks.at(id).cn_cont = 1. / block_data[0];
818
819 if (block_data.size() < 2)
820 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
821 "provide an appropriate number of entries (2) parameters for "
822 "MAT_ELASTIC block");
823 }
824 }
825
826 if (mat_blocks.size() < 2)
827 mat_blocks[0] = BlockParamData();
828
829 cache = &mat_blocks.begin()->second;
830
831 MoFEMFunctionReturn(0);
832}
833
834MoFEMErrorCode ContactPlasticity::setUP() {
835 MoFEMFunctionBegin;
836 Simple *simple = mField.getInterface<Simple>();
837
838 auto get_base = [this](bool has_hexes = false) {
839 MoFEMFunctionBeginHot;
840
841 // Select base
842 // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
843 // Demkowicz base. We need to implement Demkowicz H1 base on tet.
844 enum bases { AINSWORTH, DEMKOWICZ, BERNSTEIN, LASBASETOPT };
845 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz",
846 "bernstein"};
847 PetscInt choice_base_value = has_hexes ? DEMKOWICZ : AINSWORTH;
848 if (has_hexes)
849 choice_base_value = DEMKOWICZ;
850 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
851 LASBASETOPT, &choice_base_value, PETSC_NULL);
852
853 switch (choice_base_value) {
854 case AINSWORTH:
855 base = AINSWORTH_LEGENDRE_BASE;
856 MOFEM_LOG("WORLD", Sev::inform)
857 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
858 break;
859 case DEMKOWICZ:
860 base = DEMKOWICZ_JACOBI_BASE;
861 MOFEM_LOG("WORLD", Sev::inform)
862 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
863 break;
864 case BERNSTEIN:
865 base = AINSWORTH_BERNSTEIN_BEZIER_BASE;
866 MOFEM_LOG("WORLD", Sev::inform)
867 << "Set AINSWORTH_BERNSTEIN_BEZIER_BASE for displacements";
868 break;
869 default:
870 base = LASTBASE;
871 break;
872 }
873 MoFEMFunctionReturnHot(0);
874 };
875
876 auto skin_ents = getEntsOnMeshSkin();
877 data.skinEnts = skin_ents;
878
879 Range tets, nodes;
880 int nb_tets, nb_hexes;
881 CHKERR mField.get_moab().get_entities_by_dimension(0, 3, tets);
882 CHKERR mField.get_moab().get_connectivity(&*tets.begin(), 1, nodes);
883 CHKERR mField.get_moab().get_number_entities_by_type(0, MBTET, nb_tets, true);
884 CHKERR mField.get_moab().get_number_entities_by_type(0, MBHEX, nb_hexes,
885 true);
886
887 if ((nb_tets && nodes.size() > 4) || (nb_hexes && nodes.size() > 8)) {
888 data.is_ho_geometry = PETSC_TRUE;
889 CHKERR simple->addDataField("MESH_NODE_POSITIONS", H1, base, 3);
890 CHKERR simple->setFieldOrder("MESH_NODE_POSITIONS", 2);
891 // auto all_edges_ptr = boost::make_shared<Range>();
892 // CHKERR mField.get_moab().get_entities_by_type(0, MBEDGE, *all_edges_ptr,
893 // true);
894 // CHKERR simple->setFieldOrder("MESH_NODE_POSITIONS", 2,
895 // all_edges_ptr.get());
896 }
897
898 CHKERR get_base(nb_hexes);
899 // Add field
900 CHKERR simple->addDomainField("U", H1, base, 3);
901 CHKERR simple->addBoundaryField("U", H1, base, 3);
902 CHKERR simple->setFieldOrder("U", data.order);
903
904 if (data.with_contact) {
905 CHKERR simple->addDomainField("SIGMA", HDIV, DEMKOWICZ_JACOBI_BASE, 3);
906 CHKERR simple->addBoundaryField("SIGMA", HDIV, DEMKOWICZ_JACOBI_BASE, 3);
907 CHKERR simple->setFieldOrder("SIGMA", 0);
908 CHKERR simple->setFieldOrder("SIGMA", data.order - 1, &skin_ents);
909 }
910 if (data.with_plasticity) {
911 if (data.test_user_base_tau) {
912 CHKERR simple->addDomainField("TAU", space_test, USER_BASE, 1);
913 CHKERR simple->setFieldOrder("TAU", std::max(0, data.order - 1));
914 CHKERR simple->addDomainField("EP", space_test, USER_BASE, 6);
915 CHKERR simple->setFieldOrder("EP", std::max(0, data.order - 1));
916
917 auto *pipeline_mng = mField.getInterface<PipelineManager>();
918
919 // this creates domain_elements
920 pipeline_mng->setDomainLhsIntegrationRule(
921 [](int, int, int approx_order) { return 1; });
922 pipeline_mng->setDomainRhsIntegrationRule(
923 [](int, int, int approx_order) { return 1; });
924
925 // get access to "TAU" data structure
926 for (string name : {"TAU", "EP"}) {
927
928 auto field_ptr = mField.get_field_structure(name);
929 // get table associating number of dofs to entities depending on
930 // approximation order set on those entities.
931 auto field_order_table =
932 const_cast<Field *>(field_ptr)->getFieldOrderTable();
933
934 // function set zero number of dofs
935 auto get_tau_bubble_order_zero = [](int p) { return 0; };
936 // function set non-zero number of dofs on tetrahedrons
937 auto get_tau_bubble_order_tet = [](int p) -> int {
938 return std::min(p + 1, 2);
939 };
940
941 field_order_table[MBVERTEX] = get_tau_bubble_order_zero;
942 field_order_table[MBEDGE] = get_tau_bubble_order_zero;
943 field_order_table[MBTRI] = get_tau_bubble_order_zero;
944 field_order_table[MBTET] = get_tau_bubble_order_tet;
945 const_cast<Field *>(field_ptr)->rebuildDofsOrderMap();
946 }
947
948 } else {
949 CHKERR simple->addDomainField("TAU", space_test, base, 1);
950 CHKERR simple->setFieldOrder("TAU", std::max(0, data.order - 1));
951 CHKERR simple->addDomainField("EP", space_test, base, 6);
952 CHKERR simple->setFieldOrder("EP", std::max(0, data.order - 1));
953 }
954 }
955
956 if (data.is_rotating && data.with_plasticity && false) {
957 CHKERR simple->addSkeletonField("SKELETON_LAMBDA", L2, base, 7);
958 CHKERR simple->setFieldOrder("SKELETON_LAMBDA",
959 std::max(0, data.order - 1));
960 }
961
962 CHKERR simple->defineFiniteElements();
963
964 // Add Neumann forces elements
965 CHKERR MetaNeumannForces::addNeumannBCElements(mField, "U");
966 CHKERR MetaNodalForces::addElement(mField, "U");
967 CHKERR MetaEdgeForces::addElement(mField, "U");
968
969 simple->getOtherFiniteElements().push_back("FORCE_FE");
970 simple->getOtherFiniteElements().push_back("PRESSURE_FE");
971 if (data.is_rotating && data.with_plasticity && false)
972 CHKERR simple->setSkeletonAdjacency();
973
974 // PetscBool is_partitioned = PETSC_TRUE;
975 CHKERR simple->defineProblem(PETSC_TRUE);
976 CHKERR simple->buildFields();
977
978 if (nb_hexes) {
979 CHKERR mField.set_field_order(0, MBQUAD, "U", 0);
980 CHKERR mField.set_field_order(0, MBHEX, "U", 0);
981 CHKERR mField.build_fields();
982 }
983 if (data.is_ho_geometry) {
984 Projection10NodeCoordsOnField ent_method_material(mField,
985 "MESH_NODE_POSITIONS");
986 CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
987 }
988
989 CHKERR simple->buildFiniteElements();
990 if (data.is_rotating && data.with_plasticity && false) {
991 CHKERR mField.remove_ents_from_finite_element("sFE", skin_ents);
992 CHKERR mField.build_finite_elements("sFE");
993 }
994 CHKERR simple->buildProblem();
995 if (data.with_plasticity)
996 CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
997
998 MoFEMFunctionReturn(0);
999}
1000
1001MoFEMErrorCode ContactPlasticity::createCommonData() {
1002 MoFEMFunctionBegin;
1003
1004 auto get_material_stiffness = [&]() {
1005 MoFEMFunctionBegin;
1006
1007 constexpr auto t_kd = Kronecker_Delta<int>();
1008 constexpr auto t_kds = Kronecker_Delta_symmetric<int>();
1009
1010 // auto &t_D = commonDataPtr->tD;
1011 (*cache).Is(i, j, k, l) =
1012 (0.5 * t_kd(i, k) * t_kd(j, l)) + (0.5 * t_kd(i, l) * t_kd(j, k));
1013
1014 // Tensor4<double, 3, 3, 3, 3> II;
1015 // II(i, j, k, l) = t_kd(i, k) * t_kd(j, l);
1016
1017 auto Is_ddg = tensor4_to_ddg((*cache).Is);
1018
1019 for (auto &mit : mat_blocks) {
1020
1021 cache = &mit.second;
1022 const double G = (*cache).young_modulus / (2. * (1. + (*cache).poisson));
1023 const double K =
1024 (*cache).young_modulus / (3. * (1. - 2. * (*cache).poisson));
1025
1026 // for plane stress
1027 double A = 6. * G / (3. * K + 4. * G);
1028 // is_plane_strain and 3D
1029 A = 1.;
1030
1031 auto t_D_tmp = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
1032 auto t_D_axiator =
1033 getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Axiator);
1034 auto t_D_deviator =
1035 getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Deviator);
1036
1037 t_D_axiator(i, j, k, l) =
1038 A * (K - (2. / 3.) * G) * t_kds(i, j) * t_kds(k, l);
1039 t_D_deviator(i, j, k, l) = 2 * G * ((t_kds(i, k) ^ t_kds(j, l)) / 4.);
1040
1041 t_D_tmp(i, j, k, l) = t_D_deviator(i, j, k, l) + t_D_axiator(i, j, k, l);
1042
1043 (*cache).mtD = *commonDataPtr->mtD;
1044 (*cache).scale_constraint = 1.;
1045 }
1046
1047 // commonDataPtr->isDualBase = data.is_dual_base;
1048
1049 MoFEMFunctionReturn(0);
1050 };
1051
1052 commonDataPtr = boost::make_shared<OpContactTools::CommonData>();
1053
1054 commonDataPtr->mtD = boost::make_shared<MatrixDouble>();
1055 commonDataPtr->mtD->resize(36, 1);
1056 commonDataPtr->mtD_Axiator = boost::make_shared<MatrixDouble>();
1057 commonDataPtr->mtD_Axiator->resize(36, 1);
1058 commonDataPtr->mtD_Deviator = boost::make_shared<MatrixDouble>();
1059 commonDataPtr->mtD_Deviator->resize(36, 1);
1060
1061 commonDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
1062 commonDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
1063 commonDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
1064 commonDataPtr->mPiolaStressPtr = boost::make_shared<MatrixDouble>();
1065 commonDataPtr->mPiolaStressPtr->resize(9, 1, false);
1066
1067 commonDataPtr->mDeformationGradient = boost::make_shared<MatrixDouble>();
1068 commonDataPtr->matC = boost::make_shared<MatrixDouble>();
1069 commonDataPtr->matEigenVal = boost::make_shared<MatrixDouble>();
1070 commonDataPtr->matEigenVector = boost::make_shared<MatrixDouble>();
1071 commonDataPtr->detF = boost::make_shared<VectorDouble>();
1072 commonDataPtr->materialTangent = boost::make_shared<MatrixDouble>();
1073 commonDataPtr->dE_dF_mat = boost::make_shared<MatrixDouble>();
1074 commonDataPtr->dE_dF_mat->resize(81, 1, false);
1075
1076 commonDataPtr->contactStressPtr = boost::make_shared<MatrixDouble>();
1077 commonDataPtr->contactStressDivergencePtr =
1078 boost::make_shared<MatrixDouble>();
1079 commonDataPtr->contactNormalPtr = boost::make_shared<MatrixDouble>();
1080 commonDataPtr->contactGapPtr = boost::make_shared<VectorDouble>();
1081 commonDataPtr->contactGapDiffPtr = boost::make_shared<MatrixDouble>();
1082 commonDataPtr->contactNormalDiffPtr = boost::make_shared<MatrixDouble>();
1083
1084 commonDataPtr->contactTractionPtr = boost::make_shared<MatrixDouble>();
1085 commonDataPtr->mDispPtr = boost::make_shared<MatrixDouble>();
1086
1087 commonDataPtr->plasticSurfacePtr = boost::make_shared<VectorDouble>();
1088 commonDataPtr->plasticFlowPtr = boost::make_shared<MatrixDouble>();
1089 commonDataPtr->plasticTauPtr = boost::make_shared<VectorDouble>();
1090 commonDataPtr->plasticStrainPtr = boost::make_shared<MatrixDouble>();
1091
1092 commonDataPtr->plasticTauDotPtr = boost::make_shared<VectorDouble>();
1093 commonDataPtr->plasticStrainDotPtr = boost::make_shared<MatrixDouble>();
1094
1095 commonDataPtr->plasticTauJumpPtr = boost::make_shared<VectorDouble>();
1096 commonDataPtr->plasticStrainJumpPtr = boost::make_shared<MatrixDouble>();
1097 commonDataPtr->guidingVelocityPtr = boost::make_shared<MatrixDouble>();
1098 commonDataPtr->guidingVelocityPtr->resize(3, 1, false);
1099
1100 commonDataPtr->plasticGradTauPtr = boost::make_shared<MatrixDouble>();
1101 commonDataPtr->plasticGradStrainPtr = boost::make_shared<MatrixDouble>();
1102
1103 jAc.resize(3, 3, false);
1104 invJac.resize(3, 3, false);
1105
1106 if (data.with_plasticity)
1107 CHKERR commonDataPtr->calculateDiffDeviator();
1108
1109 CHKERR get_material_stiffness();
1110 if (data.is_neohooke)
1111 commonDataPtr->isNeoHook = true;
1112
1113 MoFEMFunctionReturn(0);
1114}
1115
1116MoFEMErrorCode ContactPlasticity::bC() {
1117 MoFEMFunctionBegin;
1118 auto bc_mng = mField.getInterface<BcManager>();
1119 auto prb_mng = mField.getInterface<ProblemsManager>();
1120 auto simple = mField.getInterface<Simple>();
1121
1122 auto get_ents = [&](const std::string blockset_name) {
1123 Range fix_ents;
1124 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
1125 if (it->getName().compare(0, blockset_name.length(), blockset_name) ==
1126 0) {
1127 CHKERR mField.get_moab().get_entities_by_handle(it->meshset, fix_ents,
1128 true);
1129 }
1130 }
1131 return fix_ents;
1132 };
1133
1134 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
1135 "U", 0, 0);
1136 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
1137 "U", 1, 1);
1138 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
1139 "U", 2, 2);
1140 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
1141 "U", 0, 2);
1142
1143 if (data.with_contact) {
1144 Range boundary_ents;
1145 for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName())
1146 if (bc_mng->checkBlock(bc, "FIX_"))
1147 boundary_ents.merge(*bc.second->getBcEntsPtr());
1148 boundary_ents.merge(get_ents("NO_CONTACT"));
1149 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
1150 boundary_ents);
1151 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
1152 mField.getInterface<Simple>()->getProblemName(), "SIGMA", boundary_ents,
1153 0, 2);
1154 }
1155
1156 std::vector<DataFromBc> bcData;
1157 dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementRemoveDofsBc>(
1158 mField, "U", simple->getProblemName(), "DISPLACEMENT", true);
1159 dirichlet_bc_ptr->iNitialize();
1160
1161 auto dm = simple->getDM();
1162 // FIXME: contact with fieldsplit does not work in parallel, requires further
1163 // investigation
1164 // TODO: will come back into this in the future
1165 if (data.with_contact && data.with_plasticity && false) {
1166 CHKERR set_contact_dm(mField, dm, dmContact);
1167 CHKERR set_plastic_ep_dm(mField, dmContact, dmEP);
1168 CHKERR set_plastic_tau_dm(mField, dmEP, dmTAU);
1169 } else if (data.with_plasticity) {
1170 CHKERR set_plastic_ep_dm(mField, dm, dmEP);
1171 CHKERR set_plastic_tau_dm(mField, dmEP, dmTAU);
1172 }
1173 // CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1174
1175 // forces and pressures on surface
1176 CHKERR MetaNeumannForces::setMomentumFluxOperators(mField, neumann_forces,
1177 PETSC_NULL, "U");
1178 // nodal forces
1179 CHKERR MetaNodalForces::setOperators(mField, nodal_forces, PETSC_NULL, "U");
1180 // edge forces
1181 CHKERR MetaEdgeForces::setOperators(mField, edge_forces, PETSC_NULL, "U");
1182
1183 // FIXME: add higher order operators !!
1184 // TODO: add higher order operators (check nonlinear_dynamics)
1185 for (auto mit = neumann_forces.begin(); mit != neumann_forces.end(); mit++) {
1186 if (data.is_ho_geometry)
1187 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", mit->second->getLoopFe(),
1188 false, false);
1189 mit->second->methodsOp.push_back(
1190 new TimeForceScale("-load_history", false));
1191 mit->second->methodsOp.push_back(new LoadScale((*cache).young_modulus_inv));
1192 }
1193 for (auto mit = nodal_forces.begin(); mit != nodal_forces.end(); mit++) {
1194 if (data.is_ho_geometry)
1195 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", mit->second->getLoopFe(),
1196 false, false);
1197 mit->second->methodsOp.push_back(
1198 new TimeForceScale("-load_history", false));
1199 mit->second->methodsOp.push_back(new LoadScale((*cache).young_modulus_inv));
1200 }
1201 for (auto mit = edge_forces.begin(); mit != edge_forces.end(); mit++) {
1202 if (data.is_ho_geometry)
1203 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", mit->second->getLoopFe(),
1204 false, false);
1205 mit->second->methodsOp.push_back(
1206 new TimeForceScale("-load_history", false));
1207 mit->second->methodsOp.push_back(new LoadScale((*cache).young_modulus_inv));
1208 }
1209
1210 // dirichlet_bc_ptr->dIag = 1;
1211 dirichlet_bc_ptr->methodsOp.push_back(
1212 new TimeForceScale(data.dirichlet_history, false));
1213
1214 MoFEMFunctionReturn(0);
1215}
1216
1217MoFEMErrorCode ContactPlasticity::OPs() {
1218 MoFEMFunctionBegin;
1219 PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
1220 Simple *simple = mField.getInterface<Simple>();
1221 auto problem_manager = mField.getInterface<ProblemsManager>();
1222 auto bc_mng = mField.getInterface<BcManager>();
1223 const Field_multiIndex *fields_ptr;
1224 CHKERR mField.get_fields(&fields_ptr);
1225
1226 auto get_boundary_markers = [&](string prob, Range bound_ents, int lo,
1227 int hi) {
1228 auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
1229
1230 CHKERR problem_manager->modifyMarkDofs(
1231 prob, ROW, "U", lo, hi, ProblemsManager::MarkOP::OR, 1, *marker_ptr);
1232
1233 CHKERR problem_manager->markDofs(prob, ROW, ProblemsManager::AND,
1234 bound_ents, *marker_ptr);
1235
1236 return marker_ptr;
1237 };
1238
1239 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "MOVE_X", "U",
1240 0, 0);
1241 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "MOVE_Y", "U",
1242 1, 1);
1243 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "MOVE_Z", "U",
1244 2, 2);
1245
1246 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "ROTATE_X",
1247 "U", 1, 2);
1248 // TODO: implement this
1249 // CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(),
1250 // "ROTATE_Y",
1251 // "U", 0, 0);
1252 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "ROTATE_Z",
1253 "U", 0, 1);
1254 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "ROTATE_ALL",
1255 "U", 0, 3);
1256
1257 auto &bc_map = bc_mng->getBcMapByBlockName();
1258 boundaryMarker =
1259 bc_mng->getMergedBlocksMarker(vector<string>{"MOVE_", "ROTATE_"});
1260
1261 auto add_disp_bc_ops = [&]() {
1262 MoFEMFunctionBeginHot;
1263 std::map<char, size_t> xyz_map = {{'X', 0}, {'Y', 1}, {'Z', 2}};
1264
1265 if (data.with_plasticity)
1266 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
1267 new OpSchurBeginBoundary("U", commonDataPtr));
1268
1269 for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
1270 if (bc_mng->checkBlock(bc, "MOVE_")) {
1271
1272 int idx = xyz_map.at(bc.first[(bc.first.find("MOVE_") + 5)]);
1273
1274 auto disp_vec = bc.second->bcAttributes;
1275 if (disp_vec.empty())
1276 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1277 "provide an appropriate number of params (1) for MOVE block");
1278 VectorDouble disp({0, 0, 0});
1279 disp(idx) = disp_vec[0];
1280
1281 bcData.push_back(new DataFromMove(disp, *bc.second->getBcEntsPtr()));
1282 // rhs
1283 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
1284 new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
1285 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
1286 new OpEssentialRHS_U("U", commonDataPtr, bcData.back().scaledValues,
1287 bcData.back().ents));
1288 pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpUnSetBc("U"));
1289 // lhs
1290 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
1291 new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
1292 pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpSpringLhsNoFS(
1293 "U", "U", [](double, double, double) { return 1.; },
1294 boost::make_shared<Range>(bcData.back().ents)));
1295 pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpUnSetBc("U"));
1296 }
1297 }
1298
1299 for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
1300 if (std::regex_match(bc.first, std::regex("(.*)_ROTATE_(.*)"))) {
1301 MOFEM_LOG("WORLD", Sev::inform)
1302 << "Set boundary (rotation) residual for " << bc.first;
1303 auto angles = boost::make_shared<VectorDouble>(3);
1304 auto c_coords = boost::make_shared<VectorDouble>(3);
1305 if (bc.second->bcAttributes.size() < 6)
1306 SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1307 "Wrong size of boundary attributes vector. Set the correct "
1308 "block size attributes (3) angles and (3) coordinates for "
1309 "center of rotation. Size of attributes %d",
1310 bc.second->bcAttributes.size());
1311 std::copy(&bc.second->bcAttributes[0], &bc.second->bcAttributes[3],
1312 angles->data().begin());
1313 std::copy(&bc.second->bcAttributes[3], &bc.second->bcAttributes[6],
1314 c_coords->data().begin());
1315 bcData.push_back(new DataFromMove(*angles, *bc.second->getBcEntsPtr()));
1316 // rhs
1317 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
1318 new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
1319 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
1320 new OpRotate("U", bcData.back().scaledValues, c_coords,
1321 bc.second->getBcEntsPtr()));
1322 pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpSpringRhs(
1323 "U", commonDataPtr->mDispPtr,
1324 [](double, double, double) { return 1.; },
1325 bc.second->getBcEntsPtr()));
1326 pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpUnSetBc("U"));
1327 // lhs
1328 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
1329 new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
1330 pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpSpringLhsNoFS(
1331 "U", "U", [](double, double, double) { return 1.; },
1332 bc.second->getBcEntsPtr()));
1333 pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpUnSetBc("U"));
1334 }
1335 }
1336
1337 MoFEMFunctionReturnHot(0);
1338 };
1339
1340 auto add_ho_ops = [&](auto &pipeline) {
1341 MoFEMFunctionBeginHot;
1342 auto jac_ptr = boost::make_shared<MatrixDouble>();
1343 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1344 auto det_ptr = boost::make_shared<VectorDouble>();
1345 pipeline.push_back(new OpCalculateVectorFieldGradient<3, 3>(
1346 "MESH_NODE_POSITIONS", jac_ptr));
1347 // pipeline.push_back(new OpCalculateHOJacVolume(jac_ptr));
1348 pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
1349 pipeline.push_back(
1350 new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
1351 pipeline.push_back(new OpSetHOWeights(det_ptr));
1352 pipeline.push_back(new OpSetHOInvJacVectorBase(HDIV, inv_jac_ptr));
1353 pipeline.push_back(new OpSetHOInvJacToScalarBases<3>(H1, inv_jac_ptr));
1354 pipeline.push_back(new OpSetHOWeights(det_ptr));
1355 // pipeline.push_back(new OpMakeHighOrderGeometryWeightsOnVolume());
1356 // pipeline.push_back(new OpSetHOGeometryCoordsOnVolume());
1357 MoFEMFunctionReturnHot(0);
1358 };
1359
1360 auto add_domain_base_ops = [&](auto &pipeline) {
1361 MoFEMFunctionBeginHot;
1362 pipeline.push_back(
1363 new OpSetMaterialBlock("U", commonDataPtr, data.mapBlockTets));
1364 if (data.is_ho_geometry)
1365 CHKERR add_ho_ops(pipeline);
1366
1367 pipeline.push_back(
1368 new OpCalculateVectorFieldValues<3>("U", commonDataPtr->mDispPtr));
1369 pipeline.push_back(
1370 new OpCalculateVectorFieldGradient<3, 3>("U", commonDataPtr->mGradPtr));
1371
1372 if (data.with_plasticity) {
1373 pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
1374 "EP", commonDataPtr->plasticStrainPtr));
1375 pipeline.push_back(new OpCalculateScalarFieldValues(
1376 "TAU", commonDataPtr->plasticTauPtr));
1377 pipeline.push_back(new OpCalculateScalarFieldValuesDot(
1378 "TAU", commonDataPtr->plasticTauDotPtr));
1379 pipeline.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<3>(
1380 "EP", commonDataPtr->plasticStrainDotPtr));
1381
1382 if (data.is_rotating) {
1383 pipeline.push_back(new OpCalculateTensor2SymmetricFieldGradient<3, 3>(
1384 "EP", commonDataPtr->plasticGradStrainPtr));
1385 pipeline.push_back(new OpCalculateScalarFieldGradient<3>(
1386 "TAU", commonDataPtr->plasticGradTauPtr));
1387 pipeline.push_back(
1388 new OpCalculatePlasticConvRotatingFrame("TAU", commonDataPtr));
1389 }
1390 }
1391
1392 if (data.is_large_strain) {
1393 pipeline.push_back(new OpLogStrain("U", commonDataPtr));
1394 } else
1395 pipeline.push_back(new OpStrain("U", commonDataPtr));
1396
1397 MoFEMFunctionReturnHot(0);
1398 };
1399
1400 auto add_domain_stress_ops = [&](auto &pipeline, auto mat_D_ptr) {
1401 MoFEMFunctionBeginHot;
1402
1403 if (data.with_plasticity) {
1404 pipeline.push_back(new OpPlasticStress("U", commonDataPtr, mat_D_ptr));
1405
1406 pipeline.push_back(
1407 new OpCalculatePlasticSurfaceAndFlow("U", commonDataPtr));
1408
1409 } else
1410 pipeline.push_back(new OpStress("U", commonDataPtr, mat_D_ptr));
1411
1412 MoFEMFunctionReturnHot(0);
1413 };
1414
1415 auto add_skeleton_base_ops = [&](auto &pipeline) {
1416 volSideFe = boost::make_shared<MoFEM::DomainSideEle>(mField);
1417
1418 volSideFe->getOpPtrVector().push_back(
1419 new OpVolumeSideCalculateTAU("TAU", commonDataPtr));
1420 volSideFe->getOpPtrVector().push_back(
1421 new OpVolumeSideCalculateEP("EP", commonDataPtr));
1422 // pipeline.push_back(new OpCalculateJumpOnSkeleton("SKELETON_LAMBDA",
1423 // commonDataPtr,
1424 // volSideFe));
1425 };
1426
1427 auto add_domain_ops_lhs = [&](auto &pipeline, auto m_D_ptr) {
1428 MoFEMFunctionBeginHot;
1429
1430 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1431 if (data.with_plasticity)
1432 pipeline.push_back(new OpSchurBegin("U", commonDataPtr));
1433
1434 if (data.is_large_strain) {
1435 if (data.is_neohooke)
1436 pipeline.push_back(
1437 new OpCalculateNeoHookeStressTangent("U", commonDataPtr));
1438 else
1439 pipeline.push_back(
1440 new OpCalculateLogStressTangent<true>("U", commonDataPtr, m_D_ptr));
1441 pipeline.push_back(
1442 new OpLogStrainMatrixLhs("U", "U", commonDataPtr->materialTangent));
1443 } else
1444 pipeline.push_back(new OpStiffnessMatrixLhs("U", "U", m_D_ptr));
1445 if (data.with_plasticity) {
1446 if (data.is_large_strain)
1447 pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dEP<true>(
1448 "U", "EP", commonDataPtr, m_D_ptr));
1449 else
1450 pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dEP<false>(
1451 "U", "EP", commonDataPtr, m_D_ptr));
1452 // FIXME: dummy one for testing
1453 // TODO: this can be deleted
1454 // pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dTAU(
1455 // "U", "TAU", commonDataPtr));
1456 }
1457
1458 if (data.is_rotating)
1459 pipeline.push_back(new OpRotatingFrameLhs("U", "U", commonDataPtr));
1460
1461 pipeline.push_back(new OpUnSetBc("U"));
1462
1463 MoFEMFunctionReturnHot(0);
1464 };
1465
1466 auto add_domain_ops_rhs = [&](auto &pipeline, auto m_D_ptr) {
1467 MoFEMFunctionBeginHot;
1468
1469 auto get_body_force = [this](const double, const double, const double) {
1470 auto *pipeline_mng = mField.getInterface<PipelineManager>();
1471 auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
1472 Tensor1<double, 3> t_source;
1473 for (int i = 0; i != 3; i++)
1474 t_source(i) = (*cache).gravity[i] * (*cache).density;
1475 const auto time = fe_domain_rhs->ts_t;
1476 t_source(i) *= time;
1477 return t_source;
1478 };
1479 auto get_centrifugal_force = [this](const double x, const double y,
1480 const double z) {
1481 Tensor1<double, 3> t_coords(x, y, z);
1482 Tensor1<double, 3> t_source;
1483 t_source(i) = (*cache).density * (*cache).Omega(i, k) *
1484 (*cache).Omega(k, j) * t_coords(j);
1485 // auto *pipeline_mng = mField.getInterface<PipelineManager>();
1486 // auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
1487 // const auto time = fe_domain_rhs->ts_t;
1488 // t_source(i) *= time;
1489 return t_source;
1490 };
1491
1492 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1493
1494 if (data.is_rotating) {
1495 pipeline.push_back(new OpRotatingFrameRhs("U", commonDataPtr));
1496 // pipeline.push_back(new OpCentrifugalForce2("U",
1497 // get_centrifugal_force));
1498 }
1499
1500 pipeline.push_back(new OpBodyForce("U", get_body_force));
1501 if (data.is_large_strain) {
1502 if (data.is_neohooke)
1503 pipeline.push_back(
1504 new OpCalculateNeoHookeStressTangent("U", commonDataPtr));
1505 else
1506 pipeline.push_back(new OpCalculateLogStressTangent<false>(
1507 "U", commonDataPtr, m_D_ptr));
1508 pipeline.push_back(
1509 new OpInternalForceRhs<true, false>("U", commonDataPtr));
1510 } else
1511 pipeline.push_back(
1512 new OpInternalForceRhs<false, false>("U", commonDataPtr));
1513
1514 pipeline.push_back(new OpUnSetBc("U"));
1515 MoFEMFunctionReturnHot(0);
1516 };
1517
1518 auto add_domain_ops_rhs_constraints = [&](auto &pipeline) {
1519 MoFEMFunctionBeginHot;
1520 if (data.with_plasticity) {
1521 pipeline.push_back(new OpCalculatePlasticFlowRhs("EP", commonDataPtr));
1522 pipeline.push_back(new OpCalculateConstraintRhs("TAU", commonDataPtr));
1523 if (data.debug_info)
1524 pipeline.push_back(new OpGetGaussPtsPlasticState("U", commonDataPtr));
1525 }
1526 MoFEMFunctionReturnHot(0);
1527 };
1528
1529 auto add_domain_ops_lhs_constraints = [&](auto &pipeline, auto m_D_ptr) {
1530 MoFEMFunctionBeginHot;
1531
1532 if (data.with_plasticity) {
1533
1534 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1535
1536 if (data.is_large_strain) {
1537 pipeline.push_back(new OpCalculateLogStressTangent<false>(
1538 "U", commonDataPtr, m_D_ptr));
1539 pipeline.push_back(
1540 new OpCalculatePlasticFlowLhs_dU<true>("EP", "U", commonDataPtr));
1541 pipeline.push_back(
1542 new OpCalculateConstraintLhs_dU<true>("TAU", "U", commonDataPtr));
1543 } else {
1544 pipeline.push_back(
1545 new OpCalculatePlasticFlowLhs_dU<false>("EP", "U", commonDataPtr));
1546 pipeline.push_back(
1547 new OpCalculateConstraintLhs_dU<false>("TAU", "U", commonDataPtr));
1548 }
1549
1550 pipeline.push_back(
1551 new OpCalculatePlasticFlowLhs_dEP("EP", "EP", commonDataPtr));
1552 pipeline.push_back(
1553 new OpCalculatePlasticFlowLhs_dTAU("EP", "TAU", commonDataPtr));
1554 pipeline.push_back(
1555 new OpCalculateConstraintLhs_dEP("TAU", "EP", commonDataPtr));
1556 pipeline.push_back(
1557 new OpCalculateConstraintLhs_dTAU("TAU", "TAU", commonDataPtr));
1558
1559 if (data.with_plasticity) {
1560 pipeline.push_back(new OpSchurPreconditionMassTau(
1561 "TAU", "TAU", [](double, double, double) { return 1.; },
1562 commonDataPtr));
1563 pipeline.push_back(new OpSchurEnd("U", commonDataPtr));
1564 }
1565 pipeline.push_back(new OpUnSetBc("U"));
1566 }
1567
1568 MoFEMFunctionReturnHot(0);
1569 };
1570
1571 auto add_domain_ops_contact_lhs = [&](auto &pipeline) {
1572 if (data.with_contact) {
1573 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1574
1575 pipeline.push_back(new OpMixDivULhs(
1576 "SIGMA", "U", []() { return 1.; }, true));
1577 pipeline.push_back(new OpLambdaGraULhs(
1578 "SIGMA", "U", []() { return 1.; }, true));
1579
1580 pipeline.push_back(new OpSchurPreconditionMassContactVol(
1581 "SIGMA", "SIGMA", [](double, double, double) { return 1.; },
1582 commonDataPtr));
1583
1584 pipeline.push_back(new OpUnSetBc("U"));
1585 }
1586 };
1587
1588 auto add_domain_ops_contact_rhs = [&](auto &pipeline) {
1589 if (data.with_contact) {
1590 pipeline.push_back(new OpCalculateHVecTensorField<3, 3>(
1591 "SIGMA", commonDataPtr->contactStressPtr));
1592 pipeline.push_back(new OpCalculateHVecTensorDivergence<3, 3>(
1593 "SIGMA", commonDataPtr->contactStressDivergencePtr));
1594 pipeline.push_back(
1595 new OpMixDivURhs("SIGMA", commonDataPtr->mDispPtr,
1596 [](double, double, double) { return 1; }));
1597 pipeline.push_back(
1598 new OpMiXLambdaGradURhs("SIGMA", commonDataPtr->mGradPtr));
1599
1600 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1601 pipeline.push_back(new OpMixUTimesDivLambdaRhs(
1602 "U", commonDataPtr->contactStressDivergencePtr));
1603 pipeline.push_back(
1604 new OpMixUTimesLambdaRhs("U", commonDataPtr->contactStressPtr));
1605 pipeline.push_back(new OpUnSetBc("U"));
1606 }
1607 };
1608
1609 auto add_boundary_base_ops = [&](auto &pipeline) {
1610 MoFEMFunctionBeginHot;
1611 pipeline.push_back(
1612 new OpSetMaterialBlockBoundary("U", commonDataPtr, data.mapBlockTets));
1613 if (data.is_ho_geometry) {
1614 pipeline.push_back(new OpSetHOWeightsOnFace());
1615 pipeline.push_back(new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
1616 pipeline.push_back(new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
1617 }
1618 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
1619 // pipeline.push_back(new OpHOSetCovariantPiolaTransformOnFace3D(HDIV));
1620
1621 pipeline.push_back(
1622 new OpCalculateVectorFieldValues<3>("U", commonDataPtr->mDispPtr));
1623 if (data.with_contact) {
1624 pipeline.push_back(
1625 new OpGetClosestRigidBodyData(mField, "U", commonDataPtr, false));
1626 pipeline.push_back(new OpCalculateHVecTensorTrace<3, BoundaryEleOp>(
1627 "SIGMA", commonDataPtr->contactTractionPtr));
1628 }
1629 MoFEMFunctionReturnHot(0);
1630 };
1631
1632 auto add_boundary_ops_lhs = [&](auto &pipeline) {
1633 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1634
1635 // if (data.with_plasticity)
1636 // pipeline.push_back(new OpSchurBeginBoundary("U", commonDataPtr));
1637
1638 if (data.with_contact) {
1639 pipeline.push_back(
1640 new OpConstrainBoundaryLhs_dU("SIGMA", "U", commonDataPtr));
1641 pipeline.push_back(new OpConstrainBoundaryLhs_dTraction("SIGMA", "SIGMA",
1642 commonDataPtr));
1643 // FIXME: precondition matrix
1644 pipeline.push_back(new OpSchurPreconditionMassContactBound(
1645 "SIGMA", "SIGMA", [](double, double, double) { return 1.; },
1646 commonDataPtr));
1647 }
1648
1649 // pipeline.push_back(new OpSpringLhs(
1650 // "U", "U",
1651 // [this](double, double, double) { return (*cache).spring_stiffness;
1652 // }
1653 // ));
1654
1655 if (data.is_rotating) {
1656 auto vol_side_fe_ptr = boost::make_shared<MoFEM::DomainSideEle>(mField);
1657
1658 if (data.is_ho_geometry)
1659 CHKERR add_ho_ops(vol_side_fe_ptr->getOpPtrVector());
1660
1661 vol_side_fe_ptr->getOpPtrVector().push_back(
1662 new OpVolumeSideAssembleRowData("U", commonDataPtr));
1663
1664 pipeline.push_back(new OpRotatingFrameBoundaryLhs("U", "U", commonDataPtr,
1665 vol_side_fe_ptr));
1666 }
1667 if (data.with_plasticity)
1668 pipeline.push_back(new OpSchurEndBoundary("U", commonDataPtr));
1669
1670 pipeline.push_back(new OpUnSetBc("U"));
1671 };
1672
1673 auto add_boundary_ops_rhs = [&](auto &pipeline) {
1674 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1675 if (data.debug_info && data.with_contact)
1676 pipeline.push_back(new OpGetGaussPtsContactState("SIGMA", commonDataPtr));
1677
1678 if (data.with_contact)
1679 pipeline.push_back(new OpConstrainBoundaryRhs("SIGMA", commonDataPtr));
1680 // pipeline.push_back(new OpSpringRhs(
1681 // "U", commonDataPtr->mDispPtr,
1682 // [this](double, double, double) { return (*cache).spring_stiffness;
1683 // }));
1684
1685 if (data.is_rotating) {
1686 auto my_vol_side_fe_ptr =
1687 boost::make_shared<MoFEM::DomainSideEle>(mField);
1688 my_vol_side_fe_ptr->getOpPtrVector().push_back(
1689 new OpCalculateVectorFieldGradient<3, 3>("U",
1690 commonDataPtr->mGradPtr));
1691 pipeline.push_back(new OpRotatingFrameBoundaryRhs("U", commonDataPtr,
1692 my_vol_side_fe_ptr));
1693 }
1694 pipeline.push_back(new OpUnSetBc("U"));
1695 };
1696
1697 auto &tD = commonDataPtr->mtD;
1698 auto &tD_axi = commonDataPtr->mtD_Axiator;
1699 auto &tD_dev = commonDataPtr->mtD_Deviator;
1700
1701 feLambdaLhs = boost::make_shared<DomainEle>(mField);
1702 feLambdaRhs = boost::make_shared<DomainEle>(mField);
1703
1704 if (data.test_user_base_tau) {
1705 feLambdaLhs->getUserPolynomialBase() =
1706 boost::shared_ptr<BaseFunction>(new BubbleTauPolynomialBase());
1707 feLambdaRhs->getUserPolynomialBase() =
1708 boost::shared_ptr<BaseFunction>(new BubbleTauPolynomialBase());
1709 }
1710
1711 auto &lam_pipeline_lhs = feLambdaLhs->getOpPtrVector();
1712 auto &lam_pipeline_rhs = feLambdaRhs->getOpPtrVector();
1713
1714 // boundary
1715 CHKERR add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
1716 CHKERR add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
1717 CHKERR add_disp_bc_ops();
1718
1719 // pipeline for mechanical
1720 if (data.is_reduced_integration) {
1721
1722 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
1723 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
1724 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(),
1725 tD_dev);
1726 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(),
1727 tD_dev);
1728 CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline(), tD_dev);
1729 CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline(), tD_dev);
1730 CHKERR add_domain_ops_lhs_constraints(
1731 pipeline_mng->getOpDomainLhsPipeline(), tD_dev);
1732 CHKERR add_domain_ops_rhs_constraints(
1733 pipeline_mng->getOpDomainRhsPipeline());
1734
1735 CHKERR add_domain_base_ops(lam_pipeline_lhs); //
1736 CHKERR add_domain_base_ops(lam_pipeline_rhs); //
1737 CHKERR add_domain_stress_ops(lam_pipeline_lhs, tD_axi);
1738 CHKERR add_domain_stress_ops(lam_pipeline_rhs, tD_axi);
1739 CHKERR add_domain_ops_lhs(lam_pipeline_lhs, tD_axi);
1740 CHKERR add_domain_ops_rhs(lam_pipeline_rhs, tD_axi);
1741
1742 } else {
1743
1744 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
1745 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
1746 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(), tD);
1747 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(), tD);
1748 CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline(), tD);
1749 CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline(), tD);
1750
1751 // pipeline for constraints (only mu-part)
1752 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(),
1753 tD_dev);
1754 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(),
1755 tD_dev);
1756 CHKERR add_domain_ops_lhs_constraints(
1757 pipeline_mng->getOpDomainLhsPipeline(), tD_dev);
1758 CHKERR add_domain_ops_rhs_constraints(
1759 pipeline_mng->getOpDomainRhsPipeline());
1760 }
1761
1762 // contact integration volume
1763 add_domain_ops_contact_lhs(pipeline_mng->getOpDomainLhsPipeline());
1764 add_domain_ops_contact_rhs(pipeline_mng->getOpDomainRhsPipeline());
1765
1766 if (data.is_rotating && data.with_plasticity && false)
1767 add_skeleton_base_ops(pipeline_mng->getOpSkeletonRhsPipeline());
1768
1769 // contact integration on boundary
1770 add_boundary_ops_lhs(pipeline_mng->getOpBoundaryLhsPipeline());
1771 add_boundary_ops_rhs(pipeline_mng->getOpBoundaryRhsPipeline());
1772
1773 // setting integration rules
1774 auto integration_rule_bound = [](int, int, int approx_order) {
1775 return 2 * approx_order + 1;
1776 // return 3 * approx_order;
1777 };
1778
1779 auto integration_rule = [this](int, int, int approx_order) {
1780 if (base == DEMKOWICZ_JACOBI_BASE)
1781 return 2 * approx_order - 1;
1782 // return 3 * approx_order;
1783 return 2 * (approx_order - 1);
1784 };
1785
1786 auto integration_rule_lambda = [](int, int, int approx_order) {
1787 // return 3 * approx_order;
1788 return 1;
1789 };
1790
1791 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
1792 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
1793 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_bound);
1794 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_bound);
1795
1796 feLambdaLhs->getRuleHook = integration_rule_lambda;
1797 feLambdaRhs->getRuleHook = integration_rule_lambda;
1798
1799 if (data.with_contact) {
1800 // store data about the roller
1801 int def_val = 0;
1802 CHKERR mField.get_moab().tag_get_handle(
1803 "_ROLLER_NB", 1, MB_TYPE_INTEGER, commonDataPtr->rollerTag,
1804 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1805
1806 data.update_roller = boost::make_shared<BoundaryEle>(mField);
1807 auto &pipeline = data.update_roller->getOpPtrVector();
1808 if (data.is_ho_geometry) {
1809 pipeline.push_back(new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
1810 pipeline.push_back(new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
1811 }
1812 pipeline.push_back(new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
1813
1814 pipeline.push_back(
1815 new OpCalculateVectorFieldValues<3>("U", commonDataPtr->mDispPtr));
1816 pipeline.push_back(new OpCalculateHVecTensorTrace<3, BoundaryEleOp>(
1817 "SIGMA", commonDataPtr->contactTractionPtr));
1818 pipeline.push_back(
1819 new OpGetClosestRigidBodyData(mField, "U", commonDataPtr, true));
1820 if (data.debug_info)
1821 pipeline.push_back(
1822 new OpPostProcVertex("SIGMA", commonDataPtr, data.moab_debug));
1823 pipeline.push_back(new OpGetContactArea("SIGMA", commonDataPtr));
1824 data.update_roller->getRuleHook = integration_rule_bound;
1825 }
1826
1827 auto create_reactions_element = [&]() {
1828 MoFEMFunctionBeginHot;
1829 Range &reaction_ents = commonDataPtr->reactionEnts;
1830
1831 auto get_reactions_meshset_ents = [&]() {
1832 MoFEMFunctionBeginHot;
1833 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it))
1834 if (it->getMeshsetId() == data.reaction_id) {
1835 Range ents;
1836 CHKERR mField.get_moab().get_entities_by_handle(it->meshset, ents,
1837 true);
1838 CHKERR mField.get_moab().get_connectivity(ents, reaction_ents, true);
1839 CHKERR mField.get_moab().get_adjacencies(
1840 ents, 1, false, reaction_ents, moab::Interface::UNION);
1841 reaction_ents.merge(ents);
1842 break;
1843 }
1844
1845 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
1846 reaction_ents);
1847 int lents_size, ents_size;
1848 lents_size = reaction_ents.size();
1849 auto dm = simple->getDM();
1850 CHKERR MPIU_Allreduce(&lents_size, &ents_size, 1, MPIU_INT, MPIU_SUM,
1851 PetscObjectComm((PetscObject)dm));
1852 if (!ents_size /*&& data.reaction_id != -1 */) {
1853 MOFEM_LOG_C(
1854 "WORLD", Sev::warning,
1855 "Warning: Provided meshset (-reaction_id: %d) for calculating "
1856 "reactions is EMPTY! Calculating reactions for all FIX_ blocks.",
1857 data.reaction_id);
1858 for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName())
1859 if (bc_mng->checkBlock(bc, "FIX_"))
1860 reaction_ents.merge(*bc.second->getBcEntsPtr());
1861 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
1862 reaction_ents);
1863 }
1864
1865 MoFEMFunctionReturnHot(0);
1866 };
1867
1868 CHKERR get_reactions_meshset_ents();
1869
1870 double defs[] = {0, 0, 0};
1871 CHKERR mField.get_moab().tag_get_handle("_REACTION_TAG", 3, MB_TYPE_DOUBLE,
1872 commonDataPtr->reactionTag,
1873 MB_TAG_CREAT | MB_TAG_SPARSE, defs);
1874
1875 // create calculate reactions element
1876 data.calc_reactions = boost::make_shared<DomainEle>(mField);
1877
1878 std::vector<int> ghosts{0, 1, 2};
1879 commonDataPtr->reactionsVec = createSmartGhostVector(
1880 mField.get_comm(), (mField.get_comm_rank() ? 0 : 3), 3,
1881 (mField.get_comm_rank() ? 3 : 0), &*ghosts.begin());
1882
1883 data.calc_reactions->getRuleHook = integration_rule;
1884
1885 auto &pipeline = data.calc_reactions->getOpPtrVector();
1886 pipeline.push_back(
1887 new OpSetMaterialBlock("U", commonDataPtr, data.mapBlockTets));
1888 if (data.is_ho_geometry)
1889 CHKERR add_ho_ops(pipeline);
1890
1891 pipeline.push_back(
1892 new OpCalculateVectorFieldGradient<3, 3>("U", commonDataPtr->mGradPtr));
1893 if (data.is_large_strain) {
1894 pipeline.push_back(new OpLogStrain("U", commonDataPtr));
1895 } else
1896 pipeline.push_back(new OpStrain("U", commonDataPtr));
1897
1898 if (data.with_plasticity) {
1899 pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
1900 "EP", commonDataPtr->plasticStrainPtr));
1901 pipeline.push_back(
1902 new OpPlasticStress("U", commonDataPtr, commonDataPtr->mtD));
1903 } else
1904 pipeline.push_back(new OpStress("U", commonDataPtr, commonDataPtr->mtD));
1905
1906 // FIXME: include body forces
1907 // pipeline.push_back(new OpBodyForce("U", commonDataPtr, gravity));
1908 // if (data.is_rotating) {
1909 // pipeline.push_back(new OpRotatingFrameRhs("U", commonDataPtr));
1910 // pipeline.push_back(new OpCentrifugalForce2("U", commonDataPtr));
1911 // }
1912
1913 if (data.is_large_strain) {
1914 if (data.is_neohooke)
1915 pipeline.push_back(
1916 new OpCalculateNeoHookeStressTangent("U", commonDataPtr));
1917 else
1918 pipeline.push_back(new OpCalculateLogStressTangent<false>(
1919 "U", commonDataPtr, commonDataPtr->mtD));
1920 pipeline.push_back(
1921 new OpInternalForceRhs<true, true>("U", commonDataPtr));
1922 } else
1923 pipeline.push_back(
1924 new OpInternalForceRhs<false, true>("U", commonDataPtr));
1925
1926 MoFEMFunctionReturnHot(0);
1927 };
1928
1929 if (data.calculate_reactions)
1930 CHKERR create_reactions_element();
1931
1932 MoFEMFunctionReturn(0);
1933}
1934
1935MoFEMErrorCode ContactPlasticity::tsSolve() {
1936 MoFEMFunctionBegin;
1937
1938 Simple *simple = mField.getInterface<Simple>();
1939 PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
1940 ISManager *is_manager = mField.getInterface<ISManager>();
1941
1942 auto preProc =
1943 boost::make_shared<FePrePostProcess>(mField, bcData, commonDataPtr);
1944
1945 auto create_custom_ts = [&]() {
1946 auto set_dm_section = [&](auto dm) {
1947 MoFEMFunctionBegin;
1948 PetscSection section;
1949 CHKERR mField.getInterface<ISManager>()->sectionCreate(
1950 simple->getProblemName(), &section);
1951 CHKERR DMSetDefaultSection(dm, section);
1952 CHKERR DMSetDefaultGlobalSection(dm, section);
1953 CHKERR PetscSectionDestroy(&section);
1954 MoFEMFunctionReturn(0);
1955 };
1956 auto dm = simple->getDM();
1957
1958 CHKERR set_dm_section(dm);
1959 boost::shared_ptr<FEMethod> null;
1960 preProc->methodsOpForMove = boost::shared_ptr<MethodForForceScaling>(
1961 new TimeForceScale(data.move_history, false));
1962
1963 preProc->methodsOpForRollersDisp = boost::shared_ptr<MethodForForceScaling>(
1964 new TimeForceScale(data.contact_history, false));
1965 // // here we use accelerogram as position data
1966 // for (auto &roller : (*cache).rollerDataVec) {
1967 // if (!roller.positionDataParamName.empty())
1968 // preProc->methodsOpForRollersPosition.push_back(
1969 // boost::shared_ptr<MethodForForceScaling>(
1970 // new TimeAccelerogram(roller.positionDataParamName)));
1971 // if (!roller.directionDataParamName.empty())
1972 // preProc->methodsOpForRollersDirection.push_back(
1973 // boost::shared_ptr<MethodForForceScaling>(
1974 // new TimeAccelerogram(roller.directionDataParamName)));
1975 // }
1976
1977 // array<Range, 3> bc_ents;
1978 // for (auto &bdata : bcData)
1979 // bc_ents[bdata.comp].merge(bdata.ents);
1980
1981 // auto eraseRows = boost::make_shared<EraseRows>(
1982 // mField, simple->getProblemName(), bc_ents[0], bc_ents[1],
1983 // bc_ents[2]);
1984
1985 if (true) {
1986
1987 // Add element to calculate lhs of stiff part
1988 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), null,
1989 dirichlet_bc_ptr, null);
1990 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), null, preProc,
1991 null);
1992
1993 auto postproc_bound_el = [&]() {
1994 MoFEMFunctionBeginHot;
1995 if (auto bmc_ptr = block_mat_container_ptr.lock()) {
1996 auto &bmc = *bmc_ptr;
1997 bmc.clear();
1998 // print schur complement after assembling the boundary contributions
1999 // if (commonDataPtr->STauTau)
2000 // CHKERR MatView(commonDataPtr->STauTau,
2001 // PETSC_VIEWER_DRAW_WORLD);
2002 }
2003 MoFEMFunctionReturnHot(0);
2004 };
2005
2006 if (data.with_plasticity)
2007 pipeline_mng->getBoundaryLhsFE()->postProcessHook = postproc_bound_el;
2008
2009 if (pipeline_mng->getBoundaryLhsFE())
2010 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getBoundaryFEName(),
2011 pipeline_mng->getBoundaryLhsFE(), null,
2012 pipeline_mng->getBoundaryLhsFE());
2013
2014 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(),
2015 pipeline_mng->getDomainLhsFE(), null, null);
2016 if (data.is_reduced_integration)
2017 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), feLambdaLhs,
2018 null, null);
2019
2020 if (pipeline_mng->getSkeletonLhsFE())
2021 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getSkeletonFEName(),
2022 pipeline_mng->getSkeletonLhsFE(), null,
2023 null);
2024
2025 // CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), null, null,
2026 // dirichlet_bc_ptr);
2027 CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), null, null,
2028 preProc);
2029
2030 if (pipeline_mng->getDomainRhsFE()) {
2031
2032 // add dirichlet boundary conditions
2033 CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), null,
2034 dirichlet_bc_ptr, null);
2035 CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), null,
2036 preProc, null);
2037 if (data.is_reduced_integration)
2038 CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(),
2039 feLambdaRhs, null, null);
2040
2041 CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(),
2042 pipeline_mng->getDomainRhsFE(), null,
2043 null);
2044
2045 if (pipeline_mng->getBoundaryRhsFE())
2046 CHKERR DMMoFEMTSSetIFunction(dm, simple->getBoundaryFEName(),
2047 pipeline_mng->getBoundaryRhsFE(), null,
2048 null);
2049 if (pipeline_mng->getSkeletonRhsFE())
2050 CHKERR DMMoFEMTSSetIFunction(dm, simple->getSkeletonFEName(),
2051 pipeline_mng->getSkeletonRhsFE(), null,
2052 null);
2053 // Add surface forces
2054 for (auto fit = neumann_forces.begin(); fit != neumann_forces.end();
2055 fit++) {
2056 CHKERR DMMoFEMTSSetIFunction(dm, fit->first.c_str(),
2057 &fit->second->getLoopFe(), NULL, NULL);
2058 }
2059
2060 // Add edge forces
2061 for (auto fit = edge_forces.begin(); fit != edge_forces.end(); fit++) {
2062 CHKERR DMMoFEMTSSetIFunction(dm, fit->first.c_str(),
2063 &fit->second->getLoopFe(), NULL, NULL);
2064 }
2065
2066 // Add nodal forces
2067 for (auto fit = nodal_forces.begin(); fit != nodal_forces.end();
2068 fit++) {
2069 CHKERR DMMoFEMTSSetIFunction(dm, fit->first.c_str(),
2070 &fit->second->getLoopFe(), NULL, NULL);
2071 }
2072
2073 // CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), null,
2074 // null,
2075 // dirichlet_bc_ptr);
2076 CHKERR DMMoFEMTSSetIFunction(dm, simple->getDomainFEName(), null, null,
2077 preProc);
2078 }
2079 }
2080
2081 auto print_active_set = [&](std::array<int, 2> &lgauss_pts, string name,
2082 MoFEM::Interface &mField) {
2083 MoFEMFunctionBeginHot;
2084 std::vector<int> gauss_pts(2, 0);
2085 auto dm = mField.getInterface<Simple>()->getDM();
2086 CHKERR MPIU_Allreduce(lgauss_pts.data(), gauss_pts.data(), 2, MPIU_INT,
2087 MPIU_SUM, PetscObjectComm((PetscObject)dm));
2088
2089 MOFEM_LOG_C("WORLD", Sev::inform, "\t \t Active %s gauss pts: %d / %d",
2090 name.c_str(), (int)gauss_pts[0], (int)gauss_pts[1]);
2091 lgauss_pts[0] = lgauss_pts[1] = 0;
2092 MoFEMFunctionReturnHot(0);
2093 };
2094
2095 auto postproc_dom = [&]() {
2096 MoFEMFunctionBeginHot;
2097 CHKERR print_active_set(commonDataPtr->stateArrayPlast, "plastic",
2098 mField);
2099 MoFEMFunctionReturnHot(0);
2100 };
2101 auto postproc_bound = [&]() {
2102 MoFEMFunctionBeginHot;
2103 CHKERR print_active_set(commonDataPtr->stateArrayCont, "contact", mField);
2104 MoFEMFunctionReturnHot(0);
2105 };
2106
2107 if (data.debug_info) {
2108 if (data.with_plasticity)
2109 pipeline_mng->getDomainRhsFE()->postProcessHook = postproc_dom;
2110
2111 if (data.with_contact)
2112 pipeline_mng->getBoundaryRhsFE()->postProcessHook = postproc_bound;
2113 }
2114
2115 auto ts = MoFEM::createTS(mField.get_comm());
2116 CHKERR TSSetDM(ts, dm);
2117 return ts;
2118 };
2119
2120 auto solver = create_custom_ts();
2121
2122 CHKERR TSSetExactFinalTime(solver, TS_EXACTFINALTIME_MATCHSTEP);
2123
2124 auto dm = simple->getDM();
2125 auto D = smartCreateDMVector(dm);
2126
2127 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
2128 CHKERR TSSetFromOptions(solver);
2129 CHKERR TSSetSolution(solver, D);
2130
2131 CHKERR TSSetUp(solver);
2132
2133 // constexpr bool run_fieldsplit_for_boundary = true;
2134 if (true) {
2135 SNES snes;
2136 CHKERR TSGetSNES(solver, &snes);
2137 KSP ksp;
2138 CHKERR SNESGetKSP(snes, &ksp);
2139 PC pC;
2140 CHKERR KSPGetPC(ksp, &pC);
2141 PetscBool is_pcfs = PETSC_FALSE;
2142 PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_pcfs);
2143 SmartPetscObj<AO> sigma_ao; // contact ordering
2144
2145 auto make_is_field_map = [&]() {
2146 PetscSection section;
2147 CHKERR DMGetDefaultSection(dm, &section);
2148
2149 int num_fields;
2150 CHKERR PetscSectionGetNumFields(section, &num_fields);
2151
2152 const MoFEM::Problem *prb_ptr;
2153 CHKERR DMMoFEMGetProblemPtr(dm, &prb_ptr);
2154
2155 map<std::string, SmartPetscObj<IS>> is_map;
2156 for (int ff = 0; ff != num_fields; ff++) {
2157
2158 const char *field_name;
2159 CHKERR PetscSectionGetFieldName(section, ff, &field_name);
2160
2161 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
2162 prb_ptr->getName(), ROW, field_name, 0, MAX_DOFS_ON_ENTITY,
2163 is_map[field_name]);
2164 };
2165 return is_map;
2166 };
2167 auto is_map = make_is_field_map();
2168
2169 auto set_fieldsplit_on_bc = [&](PC &bc_pc, string name_prb) {
2170 MoFEMFunctionBeginHot;
2171
2172 PetscBool is_bc_field_split;
2173 PetscObjectTypeCompare((PetscObject)bc_pc, PCFIELDSPLIT,
2174 &is_bc_field_split);
2175
2176 auto block_prefix = simple->getProblemName();
2177 auto bc_mng = mField.getInterface<BcManager>();
2178
2179 auto is_all_bc =
2180 bc_mng->getBlockIS(block_prefix, "MOVE_Z", "U", name_prb, 2, 2);
2181 is_all_bc = bc_mng->getBlockIS(block_prefix, "MOVE_Y", "U", name_prb, 1,
2182 1, is_all_bc);
2183 is_all_bc = bc_mng->getBlockIS(block_prefix, "MOVE_X", "U", name_prb, 0,
2184 0, is_all_bc);
2185 is_all_bc = bc_mng->getBlockIS(block_prefix, "ROTATE_X", "U", name_prb, 1,
2186 2, is_all_bc);
2187 is_all_bc = bc_mng->getBlockIS(block_prefix, "ROTATE_Z", "U", name_prb, 0,
2188 1, is_all_bc);
2189 is_all_bc = bc_mng->getBlockIS(block_prefix, "ROTATE_ALL", "U", name_prb,
2190 0, 2, is_all_bc);
2191 is_all_bc = bc_mng->getBlockIS(block_prefix, "MOVE_ALL", "U", name_prb, 0,
2192 2, is_all_bc);
2193
2194 int is_all_bc_size;
2195 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
2196
2197 if (is_bc_field_split && is_all_bc_size) {
2198 MOFEM_LOG("WORLD", Sev::inform)
2199 << "Fieldsplit preconditioner for boundary block size: "
2200 << is_all_bc_size;
2201 CHKERR PCFieldSplitSetIS(bc_pc, NULL,
2202 is_all_bc); // boundary block
2203 //FIXME: it could be beneficial to have 'bc' and 'non_bc' in naming
2204
2205 // CHKERR PCFieldSplitSetIS(bc_pc, "bc",
2206 // is_all_bc); // boundary block
2207 // IS is_no_bc;
2208 // CHKERR PCFieldSplitGetISByIndex(bc_pc, 1, &is_no_bc);
2209 // CHKERR PCFieldSplitGetISByIndex(bc_pc, 0, &is_no_bc);
2210 // PetscBool is_equal;
2211 // CHKERR ISEqual(IS is1, IS is2, PetscBool * flg);
2212 // CHKERR PCFieldSplitSetIS(bc_pc, "non_bc",
2213 // is_no_bc); // non-boundary block
2214 // CHKERR ISView(is_no_bc, PETSC_VIEWER_STDOUT_SELF);
2215 // CHKERR ISDestroy(&is_no_bc);
2216 CHKERR PCSetUp(bc_pc);
2217 }
2218 // else {
2219 // CHKERR PCSetFromOptions(bc_pc);
2220 // CHKERR PCSetType(bc_pc, PCLU);
2221 // CHKERR PCSetUp(bc_pc);
2222 // }
2223
2224 MoFEMFunctionReturnHot(0);
2225 };
2226
2227 auto set_global_mat_and_vec = [&]() {
2228 MoFEMFunctionBeginHot;
2229 auto fe = mField.getInterface<PipelineManager>()->getDomainLhsFE();
2230 CHKERR DMCreateMatrix_MoFEM(dm, &fe->ts_B);
2231 CHKERR TSSetIFunction(solver, fe->ts_F, PETSC_NULL, PETSC_NULL);
2232 CHKERR TSSetIJacobian(solver, fe->ts_B, fe->ts_B, PETSC_NULL, PETSC_NULL);
2233 CHKERR KSPSetOperators(ksp, fe->ts_B, fe->ts_B);
2234 MoFEMFunctionReturnHot(0);
2235 };
2236
2237 // Setup fieldsplit (block) solver - optional: yes/no
2238 if (is_pcfs == PETSC_TRUE && !data.is_fieldsplit_bc_only) {
2239
2240 CHKERR set_global_mat_and_vec();
2241
2242 PetscBool is_l0_field_split = PETSC_TRUE; // level0 fieldsplit
2243 PetscBool is_l1_field_split = PETSC_FALSE; // level1 fieldsplit
2244
2245 // PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_l0_field_split);
2246
2247 auto set_fieldsplit_contact = [&]() {
2248 MoFEMFunctionBeginHot;
2249
2250 const MoFEM::Problem *fs_sig_ptr;
2251 CHKERR DMMoFEMGetProblemPtr(dmContact, &fs_sig_ptr);
2252 if (auto sig_data = fs_sig_ptr->subProblemData) {
2253
2254 is_map["E_IS_SIG"] = sig_data->rowIs;
2255 sigma_ao = sig_data->getSmartRowMap();
2256 // CHKERR AOView(sigma_ao, PETSC_VIEWER_STDOUT_SELF);
2257 CHKERR PCFieldSplitSetIS(pC, NULL, is_map["SIGMA"]);
2258 CHKERR PCFieldSplitSetIS(pC, NULL, is_map["E_IS_SIG"]);
2259 // CHKERR PCFieldSplitSetType(pC, PC_COMPOSITE_ADDITIVE);
2260 // CHKERR PCFieldSplitSetType(pC,
2261 // PC_COMPOSITE_MULTIPLICATIVE);
2262 CHKERR PCSetUp(pC);
2263
2264 // SmartPetscObj<Mat> mat_test;
2265 // CHKERR DMCreateMatrix_MoFEM(dmContact, mat_test);
2266
2267 PetscInt n;
2268 KSP *ct_ksp;
2269 CHKERR PCFieldSplitGetSubKSP(pC, &n, &ct_ksp);
2270 PC l1_pc;
2271 CHKERR KSPGetPC(ct_ksp[1], &l1_pc);
2272 PetscObjectTypeCompare((PetscObject)l1_pc, PCFIELDSPLIT,
2273 &is_l1_field_split);
2274 // important!
2275 // in case of contact analysis we swap the first preconditioner
2276 pC = l1_pc;
2277 CHKERR PetscFree(ct_ksp);
2278 }
2279 MoFEMFunctionReturnHot(0);
2280 };
2281
2282 //FIXME: contact with fieldsplit does not work in parallel, requires further investigation
2283 //TODO: will come back into this in the future
2284 if (data.with_contact && is_l0_field_split && false)
2285 CHKERR set_fieldsplit_contact();
2286
2287 PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_l1_field_split);
2288 if (data.with_plasticity && is_l1_field_split) {
2289
2290 const MoFEM::Problem *schur_ep_ptr;
2291 CHKERR DMMoFEMGetProblemPtr(dmEP, &schur_ep_ptr);
2292 if (auto ep_data = schur_ep_ptr->subProblemData) {
2293
2294 is_map["E_IS_EP"] = ep_data->rowIs;
2295 // CHKERR ISView(is_map["EP"], PETSC_VIEWER_STDOUT_SELF);
2296 if (sigma_ao)
2297 CHKERR AOApplicationToPetscIS(sigma_ao, is_map["EP"]);
2298 // CHKERR ISSort(is_map["EP"]);
2299
2300 CHKERR PCFieldSplitSetIS(pC, "ep", is_map["EP"]);
2301 CHKERR PCFieldSplitSetIS(pC, "etau", is_map["E_IS_EP"]);
2302
2303 PCCompositeType pc_type;
2304 CHKERR PCFieldSplitGetType(pC, &pc_type);
2305
2306 if (pc_type == PC_COMPOSITE_SCHUR) {
2307 CHKERR DMCreateMatrix_MoFEM(dmEP, commonDataPtr->SEpEp);
2308 commonDataPtr->aoSEpEp = ep_data->getSmartRowMap();
2309 if (sigma_ao)
2310 commonDataPtr->aoSEpEp = sigma_ao;
2311
2312 commonDataPtr->blockMatContainerPtr =
2313 boost::make_shared<BlockMatContainer>();
2314 MoFEM::block_mat_container_ptr =
2315 commonDataPtr->blockMatContainerPtr;
2316
2317 CHKERR PCFieldSplitSetSchurPre(pC, PC_FIELDSPLIT_SCHUR_PRE_USER,
2318 commonDataPtr->SEpEp);
2319
2320 CHKERR PCSetUp(pC);
2321 }
2322
2323 PetscInt n;
2324 KSP *tt_ksp;
2325 CHKERR PCFieldSplitGetSubKSP(pC, &n, &tt_ksp);
2326 PC tau_pc;
2327 CHKERR KSPGetPC(tt_ksp[1], &tau_pc);
2328 PetscBool is_tau_field_split;
2329 PetscObjectTypeCompare((PetscObject)tau_pc, PCFIELDSPLIT,
2330 &is_tau_field_split);
2331
2332 if (is_tau_field_split) {
2333 const MoFEM::Problem *schur_tau_ptr;
2334 CHKERR DMMoFEMGetProblemPtr(dmTAU, &schur_tau_ptr);
2335 if (auto tau_data = schur_tau_ptr->subProblemData) {
2336
2337 // CHKERR tau_data->getRowIs(is_map["E_IS_TAU"]);
2338 is_map["E_IS_TAU"] = tau_data->rowIs;
2339
2340 AO ep_ao;
2341 CHKERR ep_data->getRowMap(&ep_ao);
2342
2343 if (sigma_ao)
2344 CHKERR AOApplicationToPetscIS(sigma_ao, is_map["TAU"]);
2345 CHKERR AOApplicationToPetscIS(ep_ao, is_map["TAU"]);
2346
2347 CHKERR PCFieldSplitSetIS(tau_pc, "tau", is_map["TAU"]);
2348 CHKERR PCFieldSplitSetIS(tau_pc, "elastic", is_map["E_IS_TAU"]);
2349
2350 CHKERR PCFieldSplitGetType(tau_pc, &pc_type);
2351 if (pc_type == PC_COMPOSITE_SCHUR) {
2352 CHKERR DMCreateMatrix_MoFEM(dmTAU, commonDataPtr->STauTau);
2353 commonDataPtr->aoSTauTau = tau_data->getSmartRowMap();
2354
2355 CHKERR PCFieldSplitSetSchurPre(tau_pc,
2356 PC_FIELDSPLIT_SCHUR_PRE_USER,
2357 commonDataPtr->STauTau);
2358 }
2359 CHKERR PCSetUp(tau_pc);
2360
2361 PetscInt bc_n;
2362 KSP *bc_ksp;
2363 CHKERR PCFieldSplitGetSubKSP(tau_pc, &bc_n, &bc_ksp);
2364 PC bc_pc;
2365 CHKERR KSPGetPC(bc_ksp[1], &bc_pc);
2366
2367 CHKERR set_fieldsplit_on_bc(bc_pc, schur_tau_ptr->getName());
2368 CHKERR PetscFree(bc_ksp);
2369 }
2370 }
2371 CHKERR PetscFree(tt_ksp);
2372 }
2373 }
2374 } else if (is_pcfs && data.is_fieldsplit_bc_only) {
2375
2376 char mumps_options[] = "-fieldsplit_0_mat_mumps_icntl_14 1200 "
2377 "-fieldsplit_0_mat_mumps_icntl_24 1 "
2378 "-fieldsplit_0_mat_mumps_icntl_20 0 "
2379 "-fieldsplit_0_mat_mumps_icntl_13 1 "
2380 "-fieldsplit_1_mat_mumps_icntl_14 1200 "
2381 "-fieldsplit_1_mat_mumps_icntl_24 1 "
2382 "-fieldsplit_1_mat_mumps_icntl_20 0 "
2383 "-fieldsplit_1_mat_mumps_icntl_13 1";
2384 CHKERR PetscOptionsInsertString(NULL, mumps_options);
2385 CHKERR set_global_mat_and_vec();
2386 CHKERR set_fieldsplit_on_bc(pC, simple->getProblemName());
2387
2388 }
2389 // else
2390 // {
2391 // SNES snes;
2392 // CHKERR TSGetSNES(solver, &snes);
2393 // KSP ksp;
2394 // CHKERR SNESGetKSP(snes, &ksp);
2395 // PC pCc;
2396 // CHKERR KSPGetPC(ksp, &pCc);
2397 // CHKERR PCSetFromOptions(pCc);
2398 // CHKERR PCSetType(pCc, PCLU);
2399 // }
2400 }
2401
2402 auto set_section_monitor = [&]() {
2403 MoFEMFunctionBegin;
2404 SNES snes;
2405 CHKERR TSGetSNES(solver, &snes);
2406 PetscViewerAndFormat *vf;
2407 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
2408 PETSC_VIEWER_DEFAULT, &vf);
2409 CHKERR SNESMonitorSet(
2410 snes,
2411 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
2412 vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
2413 MoFEMFunctionReturn(0);
2414 };
2415
2416 auto create_post_process_element = [&]() {
2417 MoFEMFunctionBegin;
2418 postProcFe = boost::make_shared<PostProcVolumeOnRefinedMesh>(mField);
2419 postProcFe->generateReferenceElementMesh();
2420 if (data.test_user_base_tau)
2421 postProcFe->getUserPolynomialBase() =
2422 boost::shared_ptr<BaseFunction>(new BubbleTauPolynomialBase());
2423 postProcFe->getOpPtrVector().push_back(
2424 new OpSetMaterialBlock("U", commonDataPtr, data.mapBlockTets));
2425
2426 postProcFe->getOpPtrVector().push_back(
2427 new OpCalculateVectorFieldGradient<3, 3>("U", commonDataPtr->mGradPtr));
2428 if (data.is_large_strain) {
2429 postProcFe->getOpPtrVector().push_back(
2430 new OpLogStrain("U", commonDataPtr));
2431 } else
2432 postProcFe->getOpPtrVector().push_back(new OpStrain("U", commonDataPtr));
2433 if (data.with_plasticity) {
2434 postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
2435 "TAU", commonDataPtr->plasticTauPtr));
2436 postProcFe->getOpPtrVector().push_back(
2437 new OpCalculateTensor2SymmetricFieldValues<3>(
2438 "EP", commonDataPtr->plasticStrainPtr));
2439 postProcFe->getOpPtrVector().push_back(
2440 new OpPlasticStress("U", commonDataPtr, commonDataPtr->mtD));
2441 postProcFe->getOpPtrVector().push_back(
2442 new OpCalculatePlasticSurfaceAndFlow("U", commonDataPtr));
2443 } else
2444 postProcFe->getOpPtrVector().push_back(
2445 new OpStress("U", commonDataPtr, commonDataPtr->mtD));
2446
2447 if (data.is_large_strain) {
2448 if (data.is_neohooke)
2449 postProcFe->getOpPtrVector().push_back(
2450 new OpCalculateNeoHookeStressTangent("U", commonDataPtr));
2451 else
2452 postProcFe->getOpPtrVector().push_back(
2453 new OpCalculateLogStressTangent<false>("U", commonDataPtr,
2454 commonDataPtr->mtD));
2455 }
2456 if (data.with_contact) {
2457 postProcFe->getOpPtrVector().push_back(
2458 new OpCalculateHVecTensorDivergence<3, 3>(
2459 "SIGMA", commonDataPtr->contactStressDivergencePtr));
2460 postProcFe->getOpPtrVector().push_back(
2461 new OpCalculateHVecTensorField<3, 3>(
2462 "SIGMA", commonDataPtr->contactStressPtr));
2463
2464 postProcFe->getOpPtrVector().push_back(
2465 new OpPostProcContact("SIGMA", postProcFe->postProcMesh,
2466 postProcFe->mapGaussPts, commonDataPtr));
2467 }
2468 if (data.is_large_strain)
2469 postProcFe->getOpPtrVector().push_back(
2470 new OpPostProcElastic<true>("U", postProcFe->postProcMesh,
2471 postProcFe->mapGaussPts, commonDataPtr));
2472 else
2473 postProcFe->getOpPtrVector().push_back(
2474 new OpPostProcElastic<false>("U", postProcFe->postProcMesh,
2475 postProcFe->mapGaussPts, commonDataPtr));
2476
2477 if (data.with_plasticity)
2478 postProcFe->getOpPtrVector().push_back(
2479 new OpPostProcPlastic("U", postProcFe->postProcMesh,
2480 postProcFe->mapGaussPts, commonDataPtr));
2481
2482 if (data.calculate_reactions)
2483 // FIXME: use addFieldValuesPostProc("U", "REACTIONS", vec_rhs); instead
2484 postProcFe->getOpPtrVector().push_back(
2485 new OpSaveReactionForces(mField, "U", postProcFe->postProcMesh,
2486 postProcFe->mapGaussPts, commonDataPtr));
2487 postProcFe->addFieldValuesPostProc("U", "DISPLACEMENT");
2488 if (data.is_ho_geometry)
2489 postProcFe->addFieldValuesPostProc("MESH_NODE_POSITIONS");
2490 MoFEMFunctionReturn(0);
2491 };
2492
2493 auto scatter_create = [&](auto coeff) {
2494 SmartPetscObj<IS> is;
2495 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
2496 ROW, "U", coeff, coeff, is);
2497 int loc_size;
2498 CHKERR ISGetLocalSize(is, &loc_size);
2499 Vec v;
2500 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
2501 VecScatter scatter;
2502 CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
2503 return std::make_tuple(SmartPetscObj<Vec>(v),
2504 SmartPetscObj<VecScatter>(scatter));
2505 };
2506
2507 auto set_time_monitor = [&]() {
2508 MoFEMFunctionBegin;
2509 boost::shared_ptr<MMonitor> monitor_ptr(
2510 new MMonitor(dm, mField, commonDataPtr, postProcFe, uXScatter,
2511 uYScatter, uZScatter));
2512 boost::shared_ptr<ForcesAndSourcesCore> null;
2513 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
2514 monitor_ptr, null, null);
2515 MoFEMFunctionReturn(0);
2516 };
2517
2518 CHKERR set_section_monitor();
2519 CHKERR create_post_process_element();
2520
2521 if (data.is_restart) {
2522
2523 CHKERR PetscPrintf(PETSC_COMM_WORLD,
2524 "Reading vector in binary from %s file...\n",
2525 data.restart_file_str.c_str());
2526 PetscViewer viewer;
2527 CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD,
2528 data.restart_file_str.c_str(), FILE_MODE_READ,
2529 &viewer);
2530 CHKERR VecLoad(D, viewer);
2531
2532 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
2533 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
2534 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
2535
2536 typedef boost::tokenizer<boost::char_separator<char>> Tokenizer;
2537 string s = data.restart_file_str;
2538 Tokenizer tok{s, boost::char_separator<char>("_")};
2539 auto it = tok.begin();
2540 const int restart_step = std::stoi((++it)->c_str());
2541 data.restart_step = restart_step;
2542 string restart_tt = *(++it);
2543 restart_tt.erase(restart_tt.length() - 4);
2544 const double restart_time = std::atof(restart_tt.c_str());
2545 CHKERR TSSetStepNumber(solver, restart_step);
2546 CHKERR TSSetTime(solver, restart_time);
2547 }
2548
2549 uXScatter = scatter_create(0);
2550 uYScatter = scatter_create(1);
2551 uZScatter = scatter_create(2);
2552 CHKERR set_time_monitor();
2553
2554 // Add iteration post-proc for debugging
2555 // boost::shared_ptr<FEMethod> null_fe;
2556 // SNES snes;
2557 // CHKERR TSGetSNES(solver, &snes);
2558 // auto write_fe = boost::make_shared<FEMethod>();
2559 // write_fe->postProcessHook = [&]() {
2560 // MoFEMFunctionBegin;
2561 // int iter;
2562 // CHKERR SNESGetIterationNumber(snes, &iter);
2563 // CHKERR postProcFe->writeFile(
2564 // "it_out_plastic_" + boost::lexical_cast<std::string>(iter) + ".h5m");
2565 // MoFEMFunctionReturn(0);
2566 // };
2567 // CHKERR DMMoFEMTSSetIJacobian(dm, simple->getDomainFEName(), postProcFe,
2568 // null_fe, write_fe);
2569
2570 CHKERR TSSolve(solver, D);
2571
2572 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
2573 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
2574 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
2575
2576 MoFEMFunctionReturn(0);
2577}
2578
2579MoFEMErrorCode ContactPlasticity::postProcess() { return 0; }
2580
2581MoFEMErrorCode ContactPlasticity::checkResults() { return 0; }
2582
2583Range ContactPlasticity::getEntsOnMeshSkin() {
2584 Range body_ents;
2585 auto meshset = mField.getInterface<Simple>()->getMeshSet();
2586 CHKERR mField.get_moab().get_entities_by_dimension(meshset, 3, body_ents);
2587 // cout << "body_ents for skin: " << body_ents.size() << endl;
2588 Skinner skin(&mField.get_moab());
2589 Range skin_tris;
2590 CHKERR skin.find_skin(0, body_ents, false, skin_tris);
2591 Range boundary_ents;
2592 ParallelComm *pcomm =
2593 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
2594 if (pcomm == NULL) {
2595 pcomm = new ParallelComm(&mField.get_moab(), mField.get_comm());
2596 }
2597 CHKERR pcomm->filter_pstatus(skin_tris, PSTATUS_SHARED | PSTATUS_MULTISHARED,
2598 PSTATUS_NOT, -1, &boundary_ents);
2599
2600 return boundary_ents;
2601};
2602
2603static char help[] = "...\n\n";
2604
2605int main(int argc, char *argv[]) {
2606
2607#include <DefaultParams.hpp>
2608
2609 MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
2610
2611 try {
2612
2613 DMType dm_name = "DMMOFEM";
2614 CHKERR DMRegister_MoFEM(dm_name);
2615
2616 moab::Core mb_instance; ///< mesh database
2617 moab::Interface &moab = mb_instance; ///< mesh database interface
2618
2619 MoFEM::Core core(moab); ///< finite element database
2620 MoFEM::Interface &m_field = core; ///< finite element database insterface
2621
2622 Simple *simple = m_field.getInterface<Simple>();
2623 simple->getProblemName() = "Multifield plasticity";
2624 CHKERR simple->getOptions();
2625
2626 ContactPlasticity ex(m_field);
2627 CHKERR ContactPlasticity::getAnalysisParameters();
2628
2629 PetscBool is_partitioned = PETSC_FALSE;
2630 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_partitioned",
2631 &is_partitioned, PETSC_NULL);
2632 if (is_partitioned) {
2633 CHKERR simple->loadFile();
2634 } else {
2635 if (m_field.get_comm_size() > 1 && true) {
2636 CHKERR simple->loadFile("", "");
2637 CHKERR ContactPlasticity::myMeshPartition(m_field);
2638 } else
2639 CHKERR simple->loadFile("");
2640 }
2641
2642 // Add meshsets with material and boundary conditions
2643 CHKERR m_field.getInterface<MeshsetsManager>()->setMeshsetFromFile();
2644
2645 CHKERR ex.loadMeshBlockSets();
2646 CHKERR cache->getOptionsFromCommandLine();
2647 if (ex.data.scale_params)
2648 for (auto &mit : mat_blocks)
2649 CHKERR mit.second.scaleParameters();
2650
2651 CHKERR ex.runProblem();
2652 }
2653 CATCH_ERRORS;
2654
2655 CHKERR MoFEM::Core::Finalize();
2656}
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
param_file
std::string param_file
Definition: DefaultParams.hpp:113
Field_multiIndex
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.
Definition: FieldMultiIndices.hpp:467
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
main
int main()
Definition: adol-c_atom.cpp:46
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
ROW
@ ROW
Definition: definitions.h:123
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MAX_DOFS_ON_ENTITY
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
LASTBASE
@ LASTBASE
Definition: definitions.h:69
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
USER_BASE
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:82
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
H1
@ H1
continuous field
Definition: definitions.h:85
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
n
FTensor::Index< 'n', SPACE_DIM > n
Definition: fluid_structure_eigenproblem.cpp:66
m
FTensor::Index< 'm', SPACE_DIM > m
Definition: fluid_structure_eigenproblem.cpp:65
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
MoFEM::DMMoFEMTSSetIFunction
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMoFEM.cpp:786
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:412
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1183
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:572
MoFEM::DMMoFEMTSSetIJacobian
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMoFEM.cpp:839
MoFEM::CoreInterface::get_fields
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
MoFEM::PipelineManager::getOpDomainLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
Definition: PipelineManager.hpp:749
MoFEM::PipelineManager::getOpBoundaryLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
Definition: PipelineManager.hpp:797
MoFEM::PipelineManager::getOpBoundaryRhsPipeline
boost::ptr_deque< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
Definition: PipelineManager.hpp:821
MoFEM::PipelineManager::getOpSkeletonRhsPipeline
boost::ptr_deque< UserDataOperator > & getOpSkeletonRhsPipeline()
Get the Op Skeleton Rhs Pipeline object.
Definition: PipelineManager.hpp:869
MoFEM::PipelineManager::getOpDomainRhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Definition: PipelineManager.hpp:773
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::CoreInterface::get_field_structure
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
MoFEM::USER_ASSEMBLE
@ USER_ASSEMBLE
Definition: FormsIntegrators.hpp:104
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
D
double D
Definition: initial_diffusion.cpp:44
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
mat_blocks
std::map< int, BlockParamData > mat_blocks
Definition: multifield_plasticity.cpp:88
OpMiXLambdaGradURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixTensorTimesGradU< 3 > OpMiXLambdaGradURhs
Definition: multifield_plasticity.cpp:59
OpSpringLhs
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpSpringLhs
Definition: multifield_plasticity.cpp:66
help
static char help[]
Definition: multifield_plasticity.cpp:2603
OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
Definition: multifield_plasticity.cpp:63
OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
Definition: multifield_plasticity.cpp:53
OpMixUTimesDivLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
Definition: multifield_plasticity.cpp:61
OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
Definition: multifield_plasticity.cpp:73
OpCentrifugalForce2
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpCentrifugalForce2
Definition: multifield_plasticity.cpp:76
OpSpringRhs
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpSpringRhs
Definition: multifield_plasticity.cpp:70
OpMixDivURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
Definition: multifield_plasticity.cpp:57
cache
BlockParamData * cache
Definition: multifield_plasticity.cpp:87
OpSpringLhsNoFS
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpSpringLhsNoFS
Definition: multifield_plasticity.cpp:68
OpLogStrainMatrixLhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, 3, 3, 1 > OpLogStrainMatrixLhs
Definition: multifield_plasticity.cpp:46
space_test
FieldSpace space_test
Definition: multifield_plasticity.cpp:105
TEST_H1_SPACE
constexpr bool TEST_H1_SPACE
Definition: multifield_plasticity.cpp:103
zero_type
EntityType zero_type
Definition: multifield_plasticity.cpp:104
OpStiffnessMatrixLhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, 3, 3, 0 > OpStiffnessMatrixLhs
Definition: multifield_plasticity.cpp:48
OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
Definition: multifield_plasticity.cpp:50
OpLambdaGraULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
Definition: multifield_plasticity.cpp:55
FTensor
Tensors class implemented by Walter Landry.
Definition: FTensor.hpp:51
FTensor::dd
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
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::DMMoFEMTSSetMonitor
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1042
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
MoFEM::set_contact_dm
MoFEMErrorCode set_contact_dm(MoFEM::Interface &m_field, SmartPetscObj< DM > &dm, SmartPetscObj< DM > &dm_contact)
Definition: MultifieldPlasticity.cpp:102
MoFEM::smartCreateDMVector
DEPRECATED auto smartCreateDMVector(DM dm)
Definition: DMMoFEM.hpp:1011
MoFEM::addHOOpsFace3D
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition: HODataOperators.hpp:779
MoFEM::block_mat_container_ptr
boost::weak_ptr< BlockMatContainer > block_mat_container_ptr
Definition: multifield_plasticity.cpp:89
MoFEM::createTS
auto createTS(MPI_Comm comm)
Definition: PetscSmartObj.hpp:245
MoFEM::createSmartGhostVector
DEPRECATED auto createSmartGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Definition: PetscSmartObj.hpp:184
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
MoFEM::set_plastic_tau_dm
MoFEMErrorCode set_plastic_tau_dm(MoFEM::Interface &m_field, SmartPetscObj< DM > &dm, SmartPetscObj< DM > &dm_plastic_tau)
Definition: MultifieldPlasticity.cpp:76
MoFEM::set_plastic_ep_dm
MoFEMErrorCode set_plastic_ep_dm(MoFEM::Interface &m_field, SmartPetscObj< DM > &dm, SmartPetscObj< DM > &dm_plastic)
Definition: MultifieldPlasticity.cpp:32
OpElasticTools::tensor4_to_ddg
Ddg< double, 3, 3 > tensor4_to_ddg(Tensor4< T, 3, 3, 3, 3 > &t)
Definition: ElasticOperators.hpp:128
OpElasticTools::OpSchurPreconditionMassContactBound
struct OpSchurPreconditionMassContact< OpMassPrecHTensorBound > OpSchurPreconditionMassContactBound
Definition: ElasticOperators.hpp:1194
OpElasticTools::OpSchurPreconditionMassContactVol
struct OpSchurPreconditionMassContact< OpMassPrecHTensorVol > OpSchurPreconditionMassContactVol
Definition: ElasticOperators.hpp:1195
print_rollers
Definition: print_rollers.py:1
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
approx_order
static constexpr int approx_order
Definition: prism_polynomial_approximation.cpp:14
BlockParamData::getOptionsFromCommandLine
MoFEMErrorCode getOptionsFromCommandLine()
Definition: BasicFeTools.hpp:105
CacheParams::young_modulus
double young_modulus
Definition: BasicFeTools.hpp:20
ContactPlasticity::MMonitor::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: multifield_plasticity.cpp:315
ContactPlasticity::MMonitor::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: multifield_plasticity.cpp:294
ContactPlasticity::MMonitor::volSideFe
boost::shared_ptr< DomainSideEle > volSideFe
Definition: multifield_plasticity.cpp:579
ContactPlasticity::MMonitor::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: multifield_plasticity.cpp:585
ContactPlasticity::MMonitor::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: multifield_plasticity.cpp:584
ContactPlasticity::MMonitor::MMonitor
MMonitor(SmartPetscObj< DM > &dm, MoFEM::Interface &m_field, boost::shared_ptr< OpContactTools::CommonData > common_data_ptr, boost::shared_ptr< PostProcVolumeOnRefinedMesh > post_proc_fe, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > ux_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uy_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uz_scatter)
Definition: multifield_plasticity.cpp:248
ContactPlasticity::MMonitor::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: multifield_plasticity.cpp:317
ContactPlasticity::MMonitor::commonDataPtr
boost::shared_ptr< OpContactTools::CommonData > commonDataPtr
Definition: multifield_plasticity.cpp:576
ContactPlasticity::MMonitor::postProcFe
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProcFe
Definition: multifield_plasticity.cpp:577
ContactPlasticity::MMonitor::postProcFeSkin
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProcFeSkin
Definition: multifield_plasticity.cpp:581
ContactPlasticity::MMonitor::postProcFeSkeleton
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProcFeSkeleton
Definition: multifield_plasticity.cpp:580
ContactPlasticity::MMonitor::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: multifield_plasticity.cpp:583
ContactPlasticity::ProblemData::calc_reactions
boost::shared_ptr< DomainEle > calc_reactions
Definition: multifield_plasticity.cpp:195
ContactPlasticity::ProblemData::contact_pipeline
boost::shared_ptr< DomainEle > contact_pipeline
Definition: multifield_plasticity.cpp:196
ContactPlasticity::ProblemData::mapBlockTets
std::map< EntityHandle, int > mapBlockTets
Definition: multifield_plasticity.cpp:204
ContactPlasticity::ProblemData::update_roller
boost::shared_ptr< BoundaryEle > update_roller
Definition: multifield_plasticity.cpp:194
ContactPlasticity::boundaryMarker
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: multifield_plasticity.cpp:139
ContactPlasticity::base
FieldApproximationBase base
Definition: multifield_plasticity.cpp:157
ContactPlasticity::uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: multifield_plasticity.cpp:131
ContactPlasticity::uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: multifield_plasticity.cpp:129
ContactPlasticity::postProcess
MoFEMErrorCode postProcess()
Definition: multifield_plasticity.cpp:2579
ContactPlasticity::bcData
boost::ptr_vector< DataFromMove > bcData
Definition: multifield_plasticity.cpp:141
ContactPlasticity::loadMeshBlockSets
MoFEMErrorCode loadMeshBlockSets()
Definition: multifield_plasticity.cpp:793
ContactPlasticity::edge_forces
boost::ptr_map< std::string, EdgeForce > edge_forces
Definition: multifield_plasticity.cpp:136
ContactPlasticity::commonDataPtr
boost::shared_ptr< OpContactTools::CommonData > commonDataPtr
Definition: multifield_plasticity.cpp:127
ContactPlasticity::mField
MoFEM::Interface & mField
Definition: multifield_plasticity.cpp:117
ContactPlasticity::uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: multifield_plasticity.cpp:130
ContactPlasticity::createCommonData
MoFEMErrorCode createCommonData()
Definition: multifield_plasticity.cpp:1001
ContactPlasticity::neumann_forces
boost::ptr_map< std::string, NeumannForcesSurface > neumann_forces
Definition: multifield_plasticity.cpp:134
ContactPlasticity::getAnalysisParameters
static MoFEMErrorCode getAnalysisParameters()
Definition: multifield_plasticity.cpp:602
ContactPlasticity::postProcFe
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProcFe
Definition: multifield_plasticity.cpp:128
ContactPlasticity::volSideFe
boost::shared_ptr< DomainSideEle > volSideFe
Definition: multifield_plasticity.cpp:142
ContactPlasticity::dmContact
SmartPetscObj< DM > dmContact
Definition: multifield_plasticity.cpp:154
ContactPlasticity::myMeshPartition
static MoFEMErrorCode myMeshPartition(MoFEM::Interface &m_field)
Definition: multifield_plasticity.cpp:704
ContactPlasticity::dirichlet_bc_ptr
boost::shared_ptr< DirichletDisplacementRemoveDofsBc > dirichlet_bc_ptr
Definition: multifield_plasticity.cpp:137
ContactPlasticity::data
static ProblemData data
Definition: multifield_plasticity.cpp:244
ContactPlasticity::ContactPlasticity
ContactPlasticity(MoFEM::Interface &m_field)
Definition: multifield_plasticity.cpp:109
ContactPlasticity::nodal_forces
boost::ptr_map< std::string, NodalForce > nodal_forces
Definition: multifield_plasticity.cpp:135
ContactPlasticity::checkResults
MoFEMErrorCode checkResults()
Definition: multifield_plasticity.cpp:2581
ContactPlasticity::runProblem
MoFEMErrorCode runProblem()
Definition: multifield_plasticity.cpp:591
ContactPlasticity::dmEP
SmartPetscObj< DM > dmEP
Definition: multifield_plasticity.cpp:152
ContactPlasticity::feLambdaLhs
boost::shared_ptr< DomainEle > feLambdaLhs
Definition: multifield_plasticity.cpp:147
ContactPlasticity::dmTAU
SmartPetscObj< DM > dmTAU
Definition: multifield_plasticity.cpp:153
ContactPlasticity::feLambdaRhs
boost::shared_ptr< DomainEle > feLambdaRhs
Definition: multifield_plasticity.cpp:146
MeasureTime::stop
chrono::high_resolution_clock::time_point stop
Definition: multifield_plasticity.cpp:33
MeasureTime::start
chrono::high_resolution_clock::time_point start
Definition: multifield_plasticity.cpp:32
MetaEdgeForces::setOperators
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name, std::string mesh_node_positions="MESH_NODE_POSITIONS")
Set integration point operators.
Definition: EdgeForce.hpp:97
MetaEdgeForces::addElement
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: EdgeForce.hpp:62
MetaNeumannForces::addNeumannBCElements
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
Definition: SurfacePressure.cpp:1974
MetaNeumannForces::setMomentumFluxOperators
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
Definition: SurfacePressure.cpp:2069
MetaNodalForces::setOperators
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NodalForce > &nodal_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: NodalForce.hpp:128
MetaNodalForces::addElement
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: NodalForce.hpp:92
MethodForForceScaling::applyScale
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Definition: MethodForForceScaling.hpp:21
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
MoFEM::BcManager::getBlockIS
SmartPetscObj< IS > getBlockIS(const std::string block_prefix, const std::string block_name, const std::string field_name, const std::string problem_name, int lo, int hi, SmartPetscObj< IS > is_expand=SmartPetscObj< IS >())
Get block IS.
Definition: BcManager.cpp:242
MoFEM::BcManager::getBcMapByBlockName
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:200
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::remove_ents_from_finite_element
virtual MoFEMErrorCode remove_ents_from_finite_element(const std::string name, const EntityHandle meshset, const EntityType type, int verb=DEFAULT_VERBOSITY)=0
remove entities from given refinement level to finite element database
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MoFEM::Field
Provide data structure for (tensor) field approximation.
Definition: FieldMultiIndices.hpp:51
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::OpCalculateHOCoords
Calculate HO coordinates at gauss points.
Definition: HODataOperators.hpp:41
MoFEM::OpCalculateHVecTensorDivergence
Calculate divergence of tonsorial field using vectorial base.
Definition: UserDataOperators.hpp:2642
MoFEM::OpCalculateHVecTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2522
MoFEM::OpCalculateHVecTensorTrace
Calculate trace of vector (Hdiv/Hcurl) space.
Definition: UserDataOperators.hpp:2709
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1298
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:83
MoFEM::OpCalculateTensor2SymmetricFieldGradient
Get field gradients at integration pts for symmetric tensorial field rank 2.
Definition: UserDataOperators.hpp:1735
MoFEM::OpCalculateTensor2SymmetricFieldValuesDot
Calculate symmetric tensor field rates ant integratio pts.
Definition: UserDataOperators.hpp:1075
MoFEM::OpCalculateTensor2SymmetricFieldValues
Calculate symmetric tensor field values at integration pts.
Definition: UserDataOperators.hpp:981
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1520
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:450
MoFEM::OpGetHONormalsOnFace
Calculate normals at Gauss points of triangle element.
Definition: HODataOperators.hpp:282
MoFEM::OpHOSetContravariantPiolaTransformOnFace3D
transform Hdiv base fluxes from reference element to physical triangle
Definition: HODataOperators.hpp:299
MoFEM::OpSetBc
Set indices on entities on finite element.
Definition: FormsIntegrators.hpp:38
MoFEM::OpSetHOContravariantPiolaTransform
Apply contravariant (Piola) transfer to Hdiv space for HO geometry.
Definition: HODataOperators.hpp:181
MoFEM::OpSetHOInvJacToScalarBases
Set inverse jacobian to base functions.
Definition: HODataOperators.hpp:73
MoFEM::OpSetHOInvJacVectorBase
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Definition: HODataOperators.hpp:91
MoFEM::OpSetHOWeights
Set inverse jacobian to base functions.
Definition: HODataOperators.hpp:159
MoFEM::OpSetHOWeightsOnFace
Modify integration weights on face to take in account higher-order geometry.
Definition: HODataOperators.hpp:123
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM::PipelineManager::getDomainRhsFE
boost::shared_ptr< FEMethod > & getDomainRhsFE()
Definition: PipelineManager.hpp:405
MoFEM::PipelineManager::setDomainRhsIntegrationRule
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:530
MoFEM::PipelineManager::setBoundaryLhsIntegrationRule
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:557
MoFEM::PipelineManager::setBoundaryRhsIntegrationRule
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:584
MoFEM::PipelineManager::setDomainLhsIntegrationRule
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:503
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
MoFEM::Problem::subProblemData
boost::shared_ptr< SubProblemData > subProblemData
Definition: ProblemsMultiIndices.hpp:113
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:348
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition: PetscSmartObj.hpp:79
MoFEM::TSMethod::ts
TS ts
time solver
Definition: LoopMethods.hpp:155
MoFEM::TSMethod::ts_t
PetscReal ts_t
time
Definition: LoopMethods.hpp:162
MoFEM::TSMethod::ts_step
PetscInt ts_step
time step number
Definition: LoopMethods.hpp:159
MoFEM::TSMethod::ts_u
Vec & ts_u
state vector
Definition: LoopMethods.hpp:165
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
OpContactTools::OpConstrainBoundaryLhs_dTraction
Definition: ContactOperators.hpp:63
OpContactTools::OpConstrainBoundaryLhs_dU
Definition: ContactOperators.hpp:52
OpContactTools::OpGetClosestRigidBodyData
Definition: ContactOperators.hpp:74
OpContactTools::OpPostProcContact
Definition: ContactOperators.hpp:584
OpContactTools::OpPostProcVertex
[Common data]
Definition: ContactOperators.hpp:32
OpElasticTools::OpEssentialRHS_U
Definition: ElasticOperators.hpp:949
OpElasticTools::OpPostProcElastic
Definition: ElasticOperators.hpp:899
OpPlasticTools::OpCalculateConstraintLhs_dEP
Definition: PlasticOperators.hpp:185
OpPlasticTools::OpCalculateConstraintLhs_dTAU
Definition: PlasticOperators.hpp:195
OpPlasticTools::OpCalculateConstraintLhs_dU
Definition: PlasticOperators.hpp:175
OpPlasticTools::OpCalculatePlasticFlowLhs_dEP
Definition: PlasticOperators.hpp:144
OpPlasticTools::OpCalculatePlasticFlowLhs_dTAU
Definition: PlasticOperators.hpp:154
OpPlasticTools::OpCalculatePlasticFlowLhs_dU
Definition: PlasticOperators.hpp:134
OpPlasticTools::OpCalculatePlasticInternalForceLhs_dEP
Definition: PlasticOperators.hpp:121
OpPlasticTools::OpPostProcPlastic
Definition: PlasticOperators.hpp:205
OpRotatingFrameTools::OpPostProcPlasticSkeleton
Definition: RotatingFrameOperators.hpp:119
TimeForceScale
Force scale operator for reading two columns.
Definition: TimeForceScale.hpp:18
Generated by Doxygen 1.9.5 and hosted at University of Glasgow